Alessio Figalli

Marilena Ligabò

Thierry Paul

Semiclassical limit for mixed states with singular and rough potentials

We consider the semiclassical limit for the Heisenberg-von Neumann equation with a potential which consists of the sum of a repulsive Coulomb potential, plus a Lipschitz potential whose gradient belongs to BV ; this assumption on the potential guarantees the well posedness of the Liouville equation in the space of bounded integrable solutions. We find sufficient conditions on the initial data to ensure that the quantum dynamics converges to the classical one. More precisely, we consider the Husimi functions of the solution of the Heisenberg-von Neumann equation, and under suitable assumptions on the initial data we prove that they converge, as ε → 0, to the unique bounded solution of the Liouville equation (locally uniformly in time).

Introduction

The aim of this paper is to study the semiclassical limit for the Heisenberg-von Neumann (quantum Liouville) equation:

   iε∂ t ρε t = [H ε , ρt ε ], ρε 0 = ρ0,ε , (1.1)
{ρ 0,ε } ε>0 being a family of uniformly bounded (with respect to ε), positive, trace class operators, and with H ε = -ε 2 2 ∆ + U . When ρ0,ε is the orthogonal projector onto ψ 0,ε ∈ L 2 (R n ), (1.2) is equivalent (up to a global phase) to the Schrödinger equation

   iε∂ t ψ ε t = -ε 2 2 ∆ψ ε t + U ψ ε t = H ε ψ ε t , ψ ε 0 = ψ 0,ε ∈ L 2 (R n ), (1.2)
We recall that the Wigner transform W ε ψ of a function ψ ∈ L 2 (R n ) is defined as

W ε ψ(x, p) := 1 (2π) n R n ψ(x + ε 2 y)ψ(x - ε 2
y)e -ipy dy, 0 AF is partially supported by the NSF grant DMS-0969962.

and the one of a density matrix ρ is defined as

W ε ρ(x, p) := 1 (2π) n R n ρ(x + ε 2 y, x - ε 2
y)e -ipy dy,

where ρ(x, x ′ ) denotes the integral kernel associated to the operator ρ.

The weak limit of the Wigner function of the solution of (1.2) or (1.1) has been studied in many articles (e.g. [START_REF] Lions | Sur les mesures de Wigner[END_REF][START_REF] Gérard | Mesures semi-classiques et ondes de Bloch. Seminaire sur les Équations aux Dérivées Partielles[END_REF][START_REF] Gérard | Homogenization limits and Wigner transforms[END_REF], and more recently in strong topology in [START_REF] Athanassoulis | Strong phase-space semiclassical asymptotics[END_REF][START_REF] Athanassoulis | Strong and weak semiclassical limits for some rough Hamiltonians[END_REF]). More precisely, it is well-known that the limit dynamics of the Schrödinger equation is related to the Liouville equation

∂ t µ + p • ∇ x µ -∇U (x) • ∇ p µ = 0, (1.4)
and, roughly speaking, the above results state that: (A) If U is of class C 2 and there exists a sequence ε k → 0 such that W ε k ρ 0,ε k converges in the sense of distribution to some (nonnegative) measure µ 0 , then W ε k ρ ε k t → (Φ t ) # µ 0 (the convergence is again in the sense of distribution), where Φ t is the (unique) flow map associated to the Hamiltonian system ẋ = p, ṗ = -∇U (x) (1.5) so that µ t := (Φ t ) # µ 0 is the unique solution to (1.4) (here and in the sequel, # denotes the push-forward, so that µ t (A) = µ 0 (Φ -1 t (A)) for all A ⊂ R 2n Borel). (B) If U is of class C 1 and there exists a sequence ε k → 0 such that the curve t → W ε k ρ ε k t converges in the sense of distribution to some curve of (nonnegative) measure t → µ t , then µ t solves (1.4).

In the present paper we want to use some recent results proved in [START_REF] Ambrosio | Almost everywhere well-posedness of continuity equations with measure initial data[END_REF][START_REF] Ambrosio | Semiclassical limit of quantum dynamics with rough potentials and well posedness of transport equations with measure initial data[END_REF] to improve the literature in two directions:

(i) By lowering the regularity assumptions of (A) on the potential in order get convergence results for a more general class of potentials, as described below.

(ii) Get rid of the "after an extraction of a subsequence" argument, due to compactness, used in most of the available proofs where one is unable to uniquely identify the limit. More precisely, in (B) above one needs to take a subsequence along which the whole curve t → W ε k ρ ε k t converges for all t in order to obtain a solution to (1.4). Moreover, the limiting solution may depend on the particular subsequence. In our case we will be able to show that, for a class potential much larger than C 2 , once one assumes that the Wigner functions at time t = 0 have a limit, then the limit at any other time will converge to a "uniquely identified" solution of (1.4).

The price to pay for the lack of regularity of the potential will be to have some size condition on the initial datum which forbids the possibility of considering pure states. Even more, the Wigner function of the initial datum cannot concentrate at a point, a possibility which might actually enter in conflict with the fact that the underlying flow is not uniquely defined everywhere. Let us mention however that, with extra assumptions on the potential (but still allowing the possibility of not having uniqueness of a classical flow), it is possible to consider concentrating initial Wigner functions, giving rise to atomic measures whose evolution follows the "multicharacteristics" of the flow (see [START_REF] Athanassoulis | Strong and weak semiclassical limits for some rough Hamiltonians[END_REF]).

As described below, we will nevertheless show that, for general bounded and globally Lipschitz potential associated to locally BV vector fields (in addition to some Coulomb part), the Wigner measure of the solution at any time is the push-forward of the initial one by the Ambrosio-DiPerna-Lions flow [START_REF] Diperna | Ordinary differential equations, transport theory and Sobolev spaces[END_REF][START_REF] Ambrosio | Transport equation and Cauchy problem for BV vector fields[END_REF].

Our method will use extensively the Husimi transforms ψ → Wε ψ and ρ → Wε ρ, which we recall are defined in terms of convolution of the Wigner transform with the 2n-dimensional Gaussian kernel with variance ε/2:

Wε ψ := (W ε ψ) * G (2n) ε , Wε ρ := (W ε ρ) * G (2n) ε , G (2n) ε (x, p) := e -(|x| 2 +|p| 2 )/ε (πε) n = G (n) ε (x)G (n) ε (p).
(1.6) Of course, the asymptotic behaviour of the Wigner and Husimi transform is the same in the limit ε → 0. However, one of the main advantages of the Husimi transform is that it is nonnegative (see Appendix).

Let us observe that, thanks to (A.8), the L ∞ -norm of Wε ψ can be estimated using the Cauchy-Schwarz inequality:

Wε ψ(x, p) ≤ 1 ε n ψ 2 L 2 φ ε x,p 2
L 2 = ψ 2 L 2 ε n .
However, this estimate blows up as ε → 0. On the other hand we will prove that, by averaging the initial condition with respect to translations, we can get a uniform estimate as ε → 0 (Section 3.2). This gives us, for instance, an important family of initial data to which our result and the ones in [START_REF] Ambrosio | Semiclassical limit of quantum dynamics with rough potentials and well posedness of transport equations with measure initial data[END_REF] apply (see also the other examples in Section 3).

2 The main results

Setting

We are concerned with the derivation of classical mechanics from quantum mechanics, corresponding to the study of the asymptotic behaviour of solutions ρε t to the Heisenberg-von Neumann equation

   iε∂ t ρε t = [H ε , ρε t ] ρε 0 = ρ0,ε , (2.1)
as ε → 0, where

H ε = -ε 2 2 ∆ + U , and U : R n → R is of the form U b + U s on R n , where U s is a repulsive Coulomb potential U s (x) = 1≤i<j≤M Z i Z j |x i -x j | , M ≤ n/3, x = (x 1 , . . . x M , x) ∈ (R 3 ) M × R n-3M , Z i > 0, U b is globally bounded, locally Lipschitz, ∇U b ∈ BV loc (R n ; R n ), and ess sup x∈R n |∇U b (x)| 1 + |x| < +∞.
The formal solution of (1.2) is ρε t , where

ρε t := e -itHε/ε ρ0,ε e -itHε/ε
and its kernel is ρ ε t . Moreover, as shown for instance in [START_REF] Lions | Sur les mesures de Wigner[END_REF], W ε ρ ε t solves in the sense of distributions the equation

∂ t W ε ρ ε t + p • ∇ x W ε ρ ε t = E ε (U, ρ ε t ), (2.2)
where E ε (U, ρ) is given by

E ε (U, ρ)(x, p) := - i (2π) n R n U (x + ε 2 y) -U (x -ε 2 y) ε ρ(x + ε 2 y, x - ε 2
y)e -ipy dy.

(2.3)

Adding and subtracting ∇U (x) • y in the term in square brackets and using ye -ip•y = i∇ p e -ip•y , an integration by parts gives

E ε (U, ρ) = ∇U (x) • ∇ p W ε ρ + E ′ ε (U, ρ), where E ′ ε (U, ρ) is given by E ′ ε (U, ρ)(x, p) := - i (2π) n R n U (x + ε 2 y) -U (x -ε 2 y) ε -∇U (x) • y ρ(x + ε 2 y, x - ε 2
y)e -ipy dy.

(2.4) Let b : R 2n → R 2n be the autonomous divergence-free vector field b(x, p) := p, -∇U (x) . Then, by the discussion above, W ε ρ ε t solves the Liouville equation associated to b with an error term:

∂ t W ε ρ ε t + b • ∇W ε ρ ε t = E ′ ε (U, ρ ε t ). (2.5)
On the other hand, thanks to (2.2), it is not difficult to prove that Wε ρ ε t solves in the sense of distributions the equation

∂ t Wε ρ ε t + p • ∇ x Wε ρ ε t = E ε (U, ρ ε t ) * G (2n) ε - √ ε∇ x • [W ε ρ ε t * Ḡ(2n) ε ], (2.6)
where Ḡ(2n)

ε (y, q) := q √ ε G (2n) ε (y, q). (2.7)
Since W ε ρ ε t and Wε ρ ε t have the same limit points as ε → 0, the heuristic idea is that in the limit ε → 0 all error terms should disappear, and we should be left with the Liouville equation (which describes the classical dynamics)

∂ t ω t + b • ∇ω t = 0 on R 2n .

Preliminary results on the Liouville equations

Under the above assumptions on U one cannot hope for a general uniqueness result in the space of measures for the Liouville equation, as this would be equivalent to uniqueness for the ODE with vector field b (see for instance [START_REF] Ambrosio | Lecture Notes in Mathematics "Calculus of Variations and Non-Linear Partial Differential Equations[END_REF]). On the other hand, as shown in [1, Theorem 6.1], the equation

   ∂ t ω t + b • ∇ω t = 0 ω 0 = ω ∈ L 1 (R 2n ) ∩ L ∞ (R 2n ) and nonnegative, (2.8)
has existence and uniqueness in the space

L ∞ + ([0, T ]; L 1 (R 2n )∩L ∞ (R 2n
)). This means that there exist a unique W :

[0, T ] → L 1 (R 2n )∩L ∞ (R 2n ), nonnegative and such that ess sup t∈[0,T ] W t L 1 (R 2n ) + W t L ∞ (R 2n ) < +∞, that solves (2.8) in the sense of distributions on [0, T ] × R 2n .
One may wonder whether, in this general setting, solutions to the transport equation can still be described using the theory of characteristics. Even if in this case one cannot solve uniquely the ODE, one can still prove that there exists a unique flow map in the "Ambrosio-DiPerna-Lions sense". Let us recall the definition of Regular Lagrangian Flow (in short RLF) in the sense of Ambrosio-DiPerna-Lions:

We say that a (continuous) family of maps Φ t : R 2n → R 2n , t ≥ 0, is a RLF associated to (1.5) if:

-Φ 0 is the identity map.

-For L 2n -a.e. (x, p), t → Φ t (x, p) is an absolutely continuous curve solving (1.5).

-For every T > 0 there exists a constant C T such that (Φ t ) # L 2n ≤ C T L 2n for all t ∈ [0, T ],

where L 2n denotes the Lebesgue measure on R 2n . Observe that, since ∇U is not Lipschitz, a priori the ODE (1.5) could have more than one solution for some initial condition. However, the approach via RLFs allows to get rid of this problem by looking at solutions to (1.5) as a whole, and under suitable assumptions on U the RLF associated to (1.5) exists, and it is unique in the following sense: assume that Φ 1 and Φ 2 are two RLFs. Then, for L 2n -a.e. (x, p), Φ 1 t (x, p) = Φ 2 t (x, p) for all t ∈ [0, +∞). In particular, as shown in [START_REF] Ambrosio | Semiclassical limit of quantum dynamics with rough potentials and well posedness of transport equations with measure initial data[END_REF]Section 6], the unique solution to (2.8) is given by

ω t L 2n = (Φ t ) # ωL 2n .
(2.9)

Hence, the idea is that, if we can ensure that any limit point of the Husimi transforms Wε ρ ε t give rives to a curve of measure belonging to

L ∞ + ([0, T ]; L 1 (R 2n ) ∩ L ∞ (R 2n
)), by the aforementioned result we would deduce that the limit is unique (once the limit initial datum is fixed), and moreover it is transported by the unique RLF. In order to get such a result we need to make some assumptions on the initial data.

Assumptions on the initial data and main theorem

Let {ρ 0,ε } ε∈(0,1) be a family of initial data which satisfy ρ0,ε = ρ * 0,ε , ρ0,ε ≥ 0 and tr(ρ

0,ε ) = 1 ∀ ε ∈ (0, 1). Let ρ0,ε = j∈N µ (ε) j φ (ε) j , • φ (ε) j
be the spectral decomposition of ρ0,ε , and denote by ρ 0,ε its integral kernel. We assume: sup

ε∈(0,1) j∈N µ (ε) j H ε φ (ε) j 2 < +∞, (2.10) 1 ε n ρ0,ε ≤ C Id, (2.11) lim R→+∞ sup ε∈(0,1) R\B (n) R ρ 0,ε (x, x) dx = 0, (2.12)
and lim R→+∞ sup ε∈(0,1)

1 (2πε) n R\B (n) R Fρ 0,ε p ε , p ε dp = 0, (2.13) where B (n)
R is the ball of radius R in R n and F is the Fourier transform on R 2n , see (A.5). Conditions (2.12) and (2.13) are equivalent to asking that the family of probability measure { Wε ρ 0,ε } ε∈(0,1) is tight (see Appendix). By Prokhorov's Theorem, this is equivalent to the compactness of { Wε ρ 0,ε } ε∈(0,1) with respect to the weak topology of probability measures (i.e., in the duality with C b (R 2n ), the space of bounded continuous functions). Hence, up to extracting a subsequence, assumptions (2.12) together with (2.13) is equivalent to the existence of a probability density ω such that

w -lim ε→0 Wε ρ 0,ε L 2n = ωL 2n ∈ P R 2n , (2.14)
where P R 2n denotes the space of probability measure on R n . In order to avoid a tedious notation which would result by working with a subsequence ε k , we will assume that (2.14) holds along the whole sequence ε → 0, keeping in mind that all the arguments could be repeated with an arbitrary subsequence.

Let us observe that condition (2.10) is slightly weaker than sup ε∈(0,1) tr(H 2 ε ρ0,ε ) < +∞, as in order to give a sense to the latter we need the operator H 2 ε ρt ε to make sense (at least on a core). Concerning assumption (2.14), let us observe that the hypothesis tr(ρ 0,ε ) = 1 implies that Wε ρ 0,ε ∈ P R 2n (see Appendix).

To express in a better and cleaner way the fact that the convergence is uniform in time, we denote by d P any bounded distance inducing the weak topology in P R 2n . Recall also that Φ t denotes the unique RLF associated to b(x, p) = (p, -∇U (x)), so that (Φ t ) # ωL 2n is the unique nonnegative solution of (2.8)

in L ∞ + ([0, T ]; L 1 (R 2n ) ∩ L ∞ (R 2n )).
Theorem 2.1. Let U be as in Section 2.1. Under the assumptions (2.10), (2.11) and (2.14)

lim ε→0 sup [0,T ] d P ( Wε ρ ε t L 2n , (Φ t ) # ωL 2n ) = 0. (2.15) Moreover, if we define W t L 2n = (Φ t ) # ωL 2n , for every smooth function ϕ ∈ C ∞ c (R 2n
) the map t → R 2n ϕW t dx dp is continuously differentiable, and

d dt R 2n ϕW t dx dp = R 2n b • ∇ϕW t dx dp.
The rest of the paper will be concerned with the proof of Theorem 2.1. However, before proceeding with the proof, we first provide some example and sufficient conditions for our result to apply.

Examples

We will give three types of examples of density matrices satisfying the assumptions of the preceding section, so that Theorem 2.1 applies.

Average of an orthonormal basis

For simplicity, we set up our first example in the one-dimensional case. In particular, there is no Coulomb interaction (that is, U = U b ), since by assumption Coulomb interactions are three-dimensional. We leave to the interested reader the extension to arbitrary dimension (the only difference in the case U s = 0 appears when checking assumption (2.10)).

Let us consider the orthonormal basis of L 2 (R) given by the (semiclassical) Hermite functions

ψ (ε) j (x) = e -x 2 /2ε 2 j j!(πε) 1/4 H j x √ ε , j ∈ N,
where H j 's are the Hermite polynomials, i.e.

H j (x) = (-1) j e x 2 d j dx j e -x 2 .
The following holds:

Proposition 3.1. Let {µ (ε)
j } j∈N be a sequence of positive numbers, and define the density matrix ρ ε given by

ρ ε = j∈N µ (ε) j ψ (ε) j , • ψ (ε) j .
Assume that (2.11), and (2.14) hold.

• 0 ≤ µ (ε) j ≤ Cε, j∈N µ (ε) j = 1; • ε 2 j∈N µ (ε) j j 2 ≤ C < +∞; • w -lim ε→0 j∈N µ (ε) j δ(x 2 + p 2 -jε) = ωL 2 ∈ P R 2 . Then (2.10),
Proof. The first assumption is equivalent to (2.11) and the trace-one condition.

Concerning (2.10), using the well-know fact that

ε d dx ψ (ε) j = ε 2 jψ (ε) j-1 -j + 1ψ (ε) j+1 ,
by a simple calculation it follows that

H ε ψ (ε) j = - ε 2 2 d 2 dx 2 ψ (ε) j + U b ψ (ε) j = - ε 4 j(j -1)ψ (ε) j-2 -(2j + 1)ψ (ε) j + (j + 1)(j + 2)ψ (ε) j+2 + U b ψ (ε) j .
Hence

H ε ψ (ε) j 2 ≤ ε 4 j(j -1) + (2j + 1) + (j + 1)(j + 2) + U b ∞ 2 ≤ C(1 + ε 2 j 2 ),
and (2.10) follows from the first two assumptions. Finally, the third assumption implies (2.14) by noticing that

w - lim ε→0,j→∞,jε→a Wε ψ (ε) j = δ(x 2 + p 2 -a)
∀ a ≥ 0 (see, for instance, [15, Exemple III.6]).

Töplitz case

Let φ ∈ H 2 (R n ; C) with R n |φ(x)| 2 dx = 1.
Given ǫ, ς > 0, for any w, q ∈ R n let ψ ε w,q be defined by

ψ ε w,q (x) := 1 ς n/2 φ x -q ς e i w•x ε .
Then, using Plancherel theorem, one can easily check that the identity

1 ε n R 2n |ψ ε w,q ψ ε w,q | dw dq = (2π) n Id (3.1)
holds, where |ψ ψ| is the Dirac notation for the orthogonal projection onto a normalized vector ψ ∈ L 2 (R n ). Thanks to (3.1) and the fact that orthogonal projectors are nonnegative operators, we immediately obtain the following important estimate: for every nonnegative bounded function

χ ε : R 2n → R, it holds 1 ε n R 2n χ ε (w, q)|ψ ε w,q ψ ε w,q | dw dq ≤ χ ε ∞ (2π) n Id. (3.2) Set now ρ0,ε := R 2n χ ε (w, q)|ψ ε w,q ψ ε w,q | dw dq, ε ∈ (0, 1),
where {χ ε } ε∈(0,1) is a family of nonnegative bounded functions such that R 2n χ ε (w, q) dw dq = 1, and let S be the singular set of U s as defined in (4.27) below.

Proposition 3.2. Let ς = ς(ε) = ε α with α ∈ (0, 1), and assume that

• sup ε∈(0,1) χ ε ∞ < +∞. • w -lim ε→0 χ ε L 2n = ωL 2n ∈ P R 2n • R 2n χ ε (w, q) |w| 4 + 1 dist(q,S) 2 dw dq ≤ C < +∞.
Then (2.10), (2.11), and (2.14) hold for the family of initial data {ρ 0,ε } ε∈(0,1) .

Proof. (2.11) follows from the first assumption and (3.2). Since ς = ε α with α ∈ (0, 1) we have that for all (w, q) ∈ R 2n

wlim ε→0 Wε ψ ε w,q L 2n = δ (w,q) , see [15, Exemple III.3], and so (2.14) follows from our second assumption.

To show that the third assumption implies (2.10), we notice that in this case (2.10) can be written as follows

R 2n χ ε (w, q) H ε ψ ε w,q , H ε ψ ε w,q dw dq < +∞. (3.3)
Since α < 1, and φ ∈ H 2 (R n ; C), by a simple computation we get

H ε ψ ε w,q , H ε ψ ε w,q ≤ ε 4 2 ∆ x ψ ε w,q , ∆ x ψ ε w,q + 2 U ψ ε w,q , U ψ ε w,q ≤ C 1 + |w| 4 + C R n U (x) 2 1 ς n φ 2 x -q ς dx.
Since U b is bounded, |U s (q)| ≤ C/dist(q, S), and R n |φ(x)| 2 dx = 1, a simple estimate analogous to the one in Section 4.4 shows that (2.10) holds. We leave the details to the interested reader.

.

Conditions on the Wigner function

Here we consider a general family of density matrices {ρ 0,ε } ε∈(0,1) which satisfies the tightness conditions (2.12) and (2.13) (so that (2.14) is satisfied up to the extraction of a subsequence).

In the next proposition we show some simple sufficient conditions on the Wigner functions {W ε ρ 0,ε } ε∈(0,1) in order to ensure the validity of assumptions (2.10) and (2.11).

Proposition 3.3. Assume that

• max |α|,|β|≤[ n 2 ]+1 ∂ α x ∂ β p W ε ρ 0,ε ∞ ≤ C < +∞, • R 2n |p| 4 4 + U 2 (x) + |p| 2 U (x) -nε 2 2 ∆U (x) W ε ρ 0,ε (x, p) dx dp ≤ C < +∞.
Then (2.10) and (2.11) hold.

Proof. Let us recall first that the Weyl symbol of an operator ρ of integral kernel ρ(x, y) is, by definition, given by

σ ε (ρ)(x, p) := R n ρ(x + y 2 , x - y 2 )e -iy•p/ε dy,
that is equal to (2πε) n W ε ρ. Moreover, using (A.3) and (A.4), it holds

tr(ρ) = R 2n W ε ρ(x, p) dx dp (3.4)
Now, we remark that the first assumption gives (2.11) using Calderón-Vaillancourt Theorem [START_REF] Boulkhemair | L 2 estimates for Weyl quantization[END_REF].

Concerning (2.10), we will prove that sup ε∈(0,1)

tr(H 2 ε ρ0,ε ) < +∞ (as observed in Section 2.3, this condition is slightly stronger than (2.10)). To this aim, we first note that

H 2 ε = ε 4 4 ∆ 2 + U 2 - ε 2 2 ∆U - ε 2 2 U ∆. (3.5)
Moreover, let us observe that if ρ1 and ρ2 have kernels ρ 1 and ρ 2 respectively, then the kernel associated to the operator ρ1 ρ2 is given by ρ 1 (•, z)ρ 2 (z, •) dz. By this fact and (3.4), a simple computation shows that the identity

tr(Aρ ε ) = R 2n σ ε (A)(x, p)W ε ρ 0,ε (x, p) dx dp
holds for any "suitable" operator A (here σ ε (A) is the Weyl symbol of A). Hence, in our case,

tr(H 2 ε ρε ) = R 2n σ ε (H 2 ε )(x, p)W ε ρ 0,ε (x, p) dx dp.
We claim that the Weyl symbol of

H 2 ε is σ ε (H 2 ε )(x, p) = |p| 4 4 + U 2 (x) + |p| 2 U (x) - nε 2 2 ∆U (x).
Indeed, let f (x, p) := |p| 2 = σ ε (-ε 2 ∆)(x, p) and g(x, p) := U (x) = σ ε (U )(x, p). Then, using Moyal expansion,

σ ε (H 2 ε )(x, p) = σ ε ε 4 4 ∆ 2 + U 2 - ε 2 2 ∆U - ε 2 2 U ∆ (x, p) = f (x, p) 2 4 + g(x, p) 2 + 1 2 f ♯g(x, p) + 1 2 g♯f (x, p),
where by definition

h 1 ♯h 2 (x, p) := e i ε 2 (∂x∂ p ′ -∂p∂ x ′ ) h 1 (x, p)h 2 (x ′ , p ′ ) x ′ =x,p ′ =p .
In our case, in the expansion of the exponential

e i ε 2 (∂x∂ p ′ -∂p∂ x ′ ) = j∈N 1 j! i ε 2 (∂ x ∂ p ′ -∂ p ∂ x ′ ) j
we can stop at the second order term, since f (x, p) = |p| 2 . Therefore

f ♯g(x, p) = |p| 2 U (x) -iεp • ∇U (x) - nε 2 2 ∆U (x),
and

g♯f (x, p) = |p| 2 U (x) + iεp • ∇U (x) - nε 2 2 ∆U (x).
This proves the claim and conclude the proof of the proposition.

Proof of Theorem 2.1

The proof of the theorem is split into several steps: first we show some basic estimates on the solutions, and we prove that the family Wε ρ ε t is tight in space and uniformly weakly continuous in time (this is the compactness part). Then we show that Wε ρ ε t solves the Liouville equation (away from the singular set of the Coulomb potential) with an error term which converges to zero as ε → 0. Combining this fact with some uniform decay estimate for Wε ρ ε t away from the singularity, we finally prove that any limit point is bounded and solves the Liouville equation. By the uniqueness of solution to the Liouville equation in the function space

L ∞ + ([0, T ]; L 1 (R 2n ) ∩ L ∞ (R 2n
)), we conclude the desired result.

Let us observe that some of our estimates can be found [START_REF] Ambrosio | Passage from quantum to classical molecular dynamics in the presence of Coulomb interactions[END_REF] and [START_REF] Ambrosio | Semiclassical limit of quantum dynamics with rough potentials and well posedness of transport equations with measure initial data[END_REF]. However, the setting and the notation there are slightly different, and in some cases one would have to recheck the details of the proofs in [START_REF] Ambrosio | Passage from quantum to classical molecular dynamics in the presence of Coulomb interactions[END_REF][START_REF] Ambrosio | Semiclassical limit of quantum dynamics with rough potentials and well posedness of transport equations with measure initial data[END_REF] to verify that everything works also in our case. Hence, for sake of completeness and in order to make this paper more accessible, we have decided to include all the details.

Basic estimates

Conserved quantities

The spectral decomposition of ρε t is

ρε t = j∈N µ (ε) j φ (ε) j,t , • φ (ε) j,t ,
where φ (ε)

j,t = e -itHε/ε φ (ε) j solves (1.2)
. By standard results on the unitary propagator e -itHε/ε follows that

j∈N µ (ε) j φ (ε) j,t , H ε φ (ε) j,t = j∈N µ (ε) j φ (ε) j , H ε φ (ε) j (4.1) and j∈N µ (ε) j H ε φ (ε) j,t 2 = j∈N µ (ε) j H ε φ (ε) j 2 (4.2)
for all t ∈ R and ε ∈ (0, 1). Therefore, using (2.10) we have sup ε∈(0,1)

sup t∈R j∈N µ (ε) j φ (ε) j,t , H ε φ (ε) j,t < +∞, (4.3)
sup ε∈(0,1)

sup t∈R j∈N µ (ε) j H ε φ (ε) j,t 2 < +∞. (4.4)

A priori estimates

From (4.1), (4.2) and from the fact that U s > 0 and U b ∈ L ∞ (R n ), follows that for all ε ∈ (0, 1)

sup t∈R R n U 2 s (x)ρ ε t (x, x) dx ≤ j∈N µ (ε) j H ε φ (ε) j 2 + 2 U b ∞   j∈N µ (ε) j φ (ε) j , H ε φ (ε) j + U b ∞   (4.5) and sup t∈R 1 2 j∈N µ (ε) j R n |ε∇φ (ε) j,y (x)| 2 dx ≤ j∈N µ (ε) j φ (ε) j H ε φ (ε) j + U b ∞ . (4.6)
Hence, by (4.3) and (4.4) we obtain sup ε∈(0,1)

sup t∈R R n U 2 s (x)ρ ε t (x, x) dx ≤ C 1 (4.7)
and sup ε∈(0,1)

sup t∈R j∈N µ (ε) j R n |ε∇φ (ε) j,t (x)| 2 dx ≤ C 2 . (4.8)

Propagation of (2.11) and consequences

Observe that, by unitarity of e itHε/ε , we have, for all t ∈ R,

1 ε n ρε t ≤ C Id. (4.9) Hence, since Wε ρ ε t (y, p) = 1 (2π) n φ ε y,p , ρε t φ ε y,p , (4.10)
(see Appendix), using (4.9) we have sup ε∈(0,1)

sup t∈R Wε ρ ε t ∞ ≤ Cε n (2π) n φ ε y,p 2 = C (2π) n (4.11) (because φ ε y,p = ε -n/2
). Now, define for all x, y ∈ R n and ε, λ > 0

g ε,λ,y (x) = ( √ 2ε) n/2 (πλ) n/4 G (n) λε 2 (x -y).
Observe that

1 ε n g ε,λ,y , ρε t g ε,λ,y = 1 ε n j∈N µ (ε) j | g ε,λ,y , φ (ε)) j,t | 2 = 1 ε n j∈N µ (ε) j |( √ 2ε) n/2 (πλ) n/4 φ (ε)) j,t * G (n) λε 2 (y)| 2 = 2 n/2 (πλ) n/2 j∈N µ (ε) j |φ (ε)) j,t * G (n) λε 2 (y)| 2 ,
therefore, since g ε,λ,y = 1, by (4.9) we have that

2 n/2 (πλ) n/2 j∈N µ (ε) j |φ (ε)) j,t * G (n) λε 2 (y)| 2 ≤ C. ( 4
sup t∈[0,T ] R 2n \C (2n) R Wε ρ t ε (x, p) dx dp = 0. (4.14)
Observe that for all R > 0 sup ε∈(0,1)

sup t∈[0,T ] R 2n \C (2n) R Wε ρ t ε (x, p) dx dp ≤ sup ε∈(0,1) sup t∈[0,T ] 1 2 (R n \C (n) R )×R n Wε ρ t ε (x, p) dx dp + R n ×(R n \C (n) R )
Wε ρ t ε (x, p) dx dp , so we can check the tightness property separately for the first and the second marginals of Wε ρ t ε . From (2.14) follows immediately that the family { Wε ρ ε,0 L 2n } ε∈(0,1) is tight (because, by Prokhorov's Theorem, a family of nonnegative finite measures on R 2n is tight if and only if it is relatively compact in the duality with C b (R 2n )). Therefore lim

R→+∞ (R n \C (n) R )×R n Wε ρ ε,0 (x, p) dx dp = 0. (4.15) Let χ ∈ C(R n ), 0 ≤ χ ≤ 1 such that χ(x) = 0 if |x| < 1/2 and χ(x) = 1 if |x| > 1, and define χ R (x) := χ(x/R). Observe that ∇χ R ∞ ≤ C ′ /R and ∆χ R ∞ ≤ C ′ /R 2 .
We define the following operator:

A (ε) R ψ(x) = χ R * G (n) ε (x)ψ(x), ψ ∈ L 2 (R n ).
Observe that d dt tr(A

(ε) R ρt ε ) = - i ε tr([A (ε) R , H ε ]ρ t ε )
and that [A

(ε) R , H ε ] = ε 2 ∆(χ R * G (n) ε )/2 + ∇(χ R * G (n) ε ) • ∇ . So, using (4.8), d dt tr(A (ε) R ρt ε ) = d dt R 2n χ R (x) Wε ρ t ε (x, p) dx dp ≤ C ′ ε R 2 + C ′ √ C 2 R ≤ C ′ R 2 + C ′ √ C 2 R , which gives (R n \C (n) 2R )×R n Wε ρ t ε (x, p) dx dp ≤ R 2n χ R (x) Wε ρ t ε (x, p) dx dp ≤ R 2n χ R (x) Wε ρ 0,ε (x, p) dx dp + C ′ R 2 + C ′ √ C 2 R T ≤ (R n \C (n) R )×R n Wε ρ 0,ε (x, p) dx dp + C ′ R 2 + C ′ √ C 2 R T.
Therefore, using (4.15), we get lim R→+∞ sup ε∈(0,1)

sup t∈[0,T ] (R n \C (n) 2R )×R n Wε ρ t ε (x, p) dx dp = 0, (4.16)
as desired. For the second marginal we observe first that

R 2n |p| 2 Wε ρ t ε (x, p) dx dp = R 2n |p| 2 W ε ρ t ε (x, p) dx dp + nε 2 (4.17) and R 2n |p| 2 W ε ρ t ε (x, p) dx dp = j∈N µ (ε) j R n 1 (2πε) ε/2 φ(ε) j,t p ε 2 |p| 2 dp = j∈N µ (ε) j R n |ε∇φ (ε) j,t (
sup t∈[0,T ] R n ×(R n \C (n) R ) Wε ρ t ε (x, p) dx dp ≤ 1 R 2 C 2 + n 2 → 0 as R → +∞.

Weak Lipschitz continuity in time

Here we prove that for all φ ∈ C ∞ c (R 2n ) the map

t ∈ R → f ε,φ (t) := R 2n φ(x, p) Wε ρ ε t (x, p) dx dp is differentiable and sup ε∈(0,1) sup t∈R d dt f ε,φ (t) ≤ C φ , (4.19)
where C φ is a constant depending only on φ. First observe that

f ε,φ (t) = R 2n W ε ρ ε t (x, p)φ ε (x, p) dx dp, (4.20)
where

φ ε := φ * G (2n) ε
. Therefore, using (2.2), we have

d dt f ε,φ (t) = R 2n E ε (U b , ρ ε t )(x, p)φ ε (x, p) dx dp + R 2n E ε (U s , ρ ε t )(x, p)φ ε (x, p) dx dp + R 2n (p • ∇ x φ ε (x, p))W ε ρ ε t (x, p) dx dp. (4.21)
For the first term it is easy to check that

R 2n E ε (U b , ρ ε t )(x, p)φ ε (x, p) dx dp ≤ ∇U b ∞ (2π) n R n |y| sup x∈R n |F p φ ε |(x, y) dy. (4.22)
In the case of the Coulomb potential we follow a specific argument borrowed from [5, proof of Theorem 1.1(ii)]), based on the inequality

1 |z + w/2| - 1 |z -w/2| ≤ |w| |z + w/2||z -w/2| (4.23) with z = (x i -x j ) ∈ R 3 , w = ε(y i -y j ) ∈ R 3
. By estimating the difference quotients of U s as in (4.23), using (4.7) we obtain

R 2n E ε (U s , ρ ε t )(x, p)φ ε (x, p) dx dp ≤ C * R n |y| sup x ′ ∈R n |F p φ ε (x ′ , y)| dy R n U 2 s (x)ρ ε t (x, x) dx ≤ C * C 1 R n |y| sup x ′ ∈R n |F p φ ε (x ′ , y)| dy, (4.24)
with C * depending only on the numbers Z 1 , . . . , Z M , and C 1 is the constant defined in (4.7).

For the last term it is easy to see that

R 2n (p • ∇ x φ ε (x, p))W ε ρ ε t (x, p) dx dp ≤ 1 (2π) n R n sup x ′ ∈R n |F p φε |(x ′ , y) dy, (4.25) where φε (x, p) = p • ∇ x φ ε (x, p).
Therefore we have only to bound

R n |y| sup x∈R n |F p φ ε (x, y)| dy and R n sup x ′ ∈R n |F p φε |(x ′ , y) dy
with a constant depending only on φ.

For the first term

R n |y| sup x∈R n |F p φ ε (x, y)| dy = R n |y| sup x∈R n R n G (n) ε (x -x ′ )F p φ(x ′ , y) dx ′ |e -y 2 ε/4 | dy ≤ R n |y| sup z∈R n |F p φ(z, y)| dy ≤ C (1)
φ .

For the second term

R n sup x ′ ∈R n |F p φε (x ′ , y)| dy = R n sup x∈R n R 3n dp dx ′ dp ′ e -ip•y φ(x ′ , p ′ )G (n) ε (p -p ′ ) p • ∇ x G (n) ε (x -x ′ ) dy. (4.26) Now observe that R 3n dp dx ′ dp ′ e -ip•y φ(x ′ , p ′ )G (n) ε (p -p ′ ) p • ∇ x G (n) ε (x -x ′ ) = n k=1 R 2n dx ′ dp ′ ∂ x k G (n) ε (x -x ′ )G (n) ε (p ′ ) R n dp p k φ(x ′ , p -p ′ )e -ip•y = e -εy 2 /4 R n dx ′ (∇ x • F p g(x -x ′ , y))G (n) ε (x ′ ) + iε 2 R 2n dx ′ (y • ∇ x F p φ(x -x ′ , y))G (n) ε (x ′ ) ,
where g(x, p) = pφ(x, p). Now, since ε ∈ (0, 1)

R n sup x ′ ∈R n |F p φε (x ′ , y)| dy ≤ R n sup x∈R n R n dx ′ (∇ • F p g(x -x ′ , y))G (n) ε (x ′ ) + ε 2 R n sup x∈R n R 2n dx ′ (y • ∇F p φ(x -x ′ , y))G (n) ε (x ′ ) ≤ R n dy sup z∈R n |∇ • F p g(z, y)| + ε 2 R n dy |y| sup z∈R n |∇ z F p φ(z, y)| ≤ C (2)
φ .

Therefore sup ε∈(0,1)

sup t∈R d dt f ε,φ (t) ≤ ∇U b ∞ (2π) n C (1) φ + C * C 1 C (1) φ + C (2) φ (2π) n .

Uniform decay away from the singularity

The singular set of U s is given by We start with the second term:

S = 1≤i<j≤M S ij , S i,j = x = (x 1 , . . . , x M , x) ∈ (R 3 ) M × R n-3M : x i = x j for some i = j , ( 4
R 2n dx dp 1 dist(x, S) 2 Wε ρ t ε (x, p) = B (n) R ×R n dx dx ′ ρ t ε (x ′ , x ′ )G (n) ε (x -x ′ ) dist(x, S) 2 ≤ R n dx ′ ρ t ε (x ′ , x ′ ) dist(x ′ , S) 2 ≤ 1 c R n dx ′ U s (x ′ ) 2 ρ t ε (x ′ , x ′ ) ≤ C 1 c ,
where c is defined in (4.28), C 1 is defined in (4.7), and we used (4.28).

To prove the second estimate we observe that

R 2n |p| 4 Wε ρ t ε (x, p) dx dp ≤ R 2n |p| 4 W ε ρ t ε (x, p) dx dp + nε 2 R 2n |p| 2 W ε ρ t ε (x, p) dx dp + n(n + 2)ε 2 4 .
Thanks to (4.18), it suffices to control the first integral in the right hand side:

R 2n |p| 4 W ε ρ t ε (x, p) dx dp = j∈N µ (ε) j R n 1 (2πε) ε/2 φ(ε) j,t p ε 2 |p| 4 dp = j∈N µ (ε) j R n |ε 2 ∆φ (ε) j,t (x)| 2 dx ≤ 2 j∈N µ (ε) j R n |H ε φ (ε) j,t (x)| 2 + U 2 (x)|φ (ε) j,t (x)| 2 dx,
and the last term is uniformly bounded thanks to (4.4), (4.7), and the boundedness of U b .

Limit continuity equation away from the singularities

We want to prove that

lim ε→0 T 0 ϕ ′ (t) R 2n φ(x, p) Wε ρ t ε (x, p) dx dp + ϕ(t) R 2n b(x, p) • ∇φ(x, p) Wε ρ t ε (x, p) dx dp dt = 0 (4.30) for all φ ∈ C ∞ c (R 2n \ (S × R n )) and ϕ ∈ C ∞ c (0, T ).
Hence, recalling (2.6), we have to show that

lim ε→0 sup t∈[0,T ] R 2n dx dp E ε (U, ρ ε t ) * G (2n) ε (x, p)φ(x, p) + R 2n dx dp ∇U (x) • ∇ p φ(x, p) Wε ρ t ε (x, p) = 0, (4.31) and lim
ε→0 T 0 dt ϕ(t) R 2n dx dp √ ε∇ x • [W ε ρ ε t * Ḡ(2n) ε ]φ(x, p) = 0, (4.32) for all φ ∈ C ∞ c (R 2n \ (S × R n )) and ϕ ∈ C ∞ c (0, T ).

Verification of (4.31)

We can consider separately the contributions of U b and U s . We start with the contribution of U s . We have to prove that

lim ε→0 sup t∈[0,T ] R 2n dx dp E ε (U s , ρ ε t ) * G (2n) ε (x, p)φ(x, p)+ R 2n dx dp ∇U s (x)•∇ p φ(x, p) Wε ρ t ε (x, p) = 0 (4.33) for all φ ∈ C ∞ c (R 2n \ (S × R n )). We know that lim ε→0 sup t∈[0,T ] R 2n ϕ(x, p)W ε ρ t ε (x, p) dx dp - R 2n ϕ(x, p) Wε ρ t ε (x, p) dx dp = 0 (4.34) for all ϕ ∈ C ∞ c (R 2n
). First of all, we see that we can apply (4.34) with ϕ(x, p) = ∇U s (x) • ∇ p φ(x, p) to replace the integrals

R 2n ∇U s (x) • ∇ p φ(x, p) Wε ψ ε dx dp with R 2n ∇U s (x) • ∇ p φ(x, p)W ε ψ ε dx dp
in the verification of (4.33). Analogously, using (4.7) and (4.24) we see that we can replace

R 2n E ε (U s , ρ t ε ) * G (2n) ε (x, p)φ(x, p) dx dp with R 2n E ε (U s , ρ t ε )(x, p)φ(x, p) dx dp.
Thus, we are led to show the convergence

lim ε→0 sup t∈[0,T ] R 2n E ε (U s , ρ t ε )φ dx dp + R 2n ∇U s (x) • ∇ p φ(x, p)W ε ρ t ε (x, p) dx dp = 0 (4.35) for all φ ∈ C ∞ c ((R n \ S) × R n ) . Since R 2n E ε (U s , ρ t ε )φ dx dp = R 2n U s (x + ε 2 y) -U s (x -ε 2 y) ε ρ t ε x + εy 2 , x - εy 2 F p φ(x, y) dxdy
we can split the region of integration in two parts, where √ ε|y| > 1 and where √ ε|y| ≤ 1. The contribution of the first region can be estimated as in (4.24), with

C * { √ ε|y|>1} |y| sup x ′ |F p φ(x ′ , y)| dy R n U 2 s (x)ρ t ε (x, x) dx,
which is infinitesimal, using (4.7) again, as ε → 0. Since

U s (x + ε 2 y) -U s (x -ε 2 y) ε → ∇U s (x) • y
uniformly as √ ε|y| ≤ 1 and x belongs to a compact subset of R n \ S, the contribution of the second part is the same of

R 2n (∇U s (x) • y)ρ t ε x + εy 2 , x - εy 2 F p φ(x, y) dxdy which coincides with R 2n ∇U s (x) • ∇ p φ(x, p)W ε ρ t ε (x, p) dx dp.
Now we consider the contribution of U b . We have to prove that lim

ε→0 sup t∈[0,T ] R 2n E ε (U b , ρ t ε )(x, p)φ ε (x, p) dx dp + R 2n ∇U b (x) • ∇ p φ(x, p) Wε ρ ε t dx dp = 0 (4.36) for all φ ∈ C ∞ c (R 2n ), where φ ε = φ * G (2n) ε
. The proof of (4.36) is divided in two parts: first we prove that

lim ε→0 sup t∈[0,T ] R 2n E ε (U b , ρ t ε )(x, p)φ(x, p) dx dp + R 2n ∇U b (x) • ∇ p φ(x, p) Wε ρ ε t dx dp = 0 (4.37)
for all φ ∈ C ∞ c (R 2n ), and then, using the following estimate

R 2n E ε (U b , ρ ε t )(x, p)ϕ(x, p) dx dp ≤ ∇U b ∞ (2π) n R n |y| sup x∈R n |F p ϕ|(x, y) dy. (4.38)
for all ϕ ∈ C ∞ c (R 2n ), we can replace φ by φ ε in the first summand of (4.37), obtaining (4.36). The proof of (4.37) is achieved by a density argument. The first remark is that linear combinations of tensor functions φ(x, p) = φ 1 (x)φ 2 (p), with φ i ∈ C ∞ c (R n ), are dense for the norm considered in (4.38). In this way, we are led to prove convergence in the case when φ(x, p) = φ 1 (x)φ 2 (p). The second remark is that convergence surely holds if U b is of class C 2 (by the arguments in [START_REF] Lions | Sur les mesures de Wigner[END_REF], [START_REF] Ambrosio | Passage from quantum to classical molecular dynamics in the presence of Coulomb interactions[END_REF]). Hence, combining the two remarks and using the linearity of the error term with respect to the potential, we can prove convergence by a density argument, by approximating U b uniformly and in W 1,2 topology on the support of φ 1 by potentials V k ∈ C 2 (R n ) with uniformly Lipschitz constants; then, setting

A k = (U b -V k )φ
sup t∈[0,T ] R 2n E ε (U b -V k , ρ ε t )(x, p)φ 1 (x)φ 2 (p) dx dp = 0.
As for the term in (4.36) involving the Husimi transforms, we can use (4.11) to obtain that lim sup k→∞ sup ε∈(0,1)

sup t∈[0,T ] R 2n Wε ρ ε t ∇(U b (x) -V k (x)) • ∇φ 2 (p)φ 1 (x) dx dp ≤ C (2π) n lim sup k→∞ R n |φ 1 (x)||∇U b (x) -∇V k (x)| dx R n |∇φ 2 (p)| dp = 0.
So we need only to prove the following lemma:

Lemma 4.1 (A priori estimate). For all λ > 0, we have that sup ε∈(0,1)

sup t∈[0,T ] R 2n E ε (U b , ρ ε t )(x, p)φ 1 (x)φ 2 (p) dx dp (4.39) ≤ φ 1 1 ∇U b ∞ sup y∈R n |y|| φ2 (y) -φ2 * G (n) λ (y)| + √ λ ∇A ∞ φ2 1 R n |u|G (n) 1 (u) du (4.40) + √ C ∇A 2 (2πλ) n/4 R n |z|| φ2 |(z) dz + U b ∞ ∇φ 1 ∞ R n |y|| φ2 * G (n) λ |(
sup t∈[0,T ] R 2n E ε (U b , ρ ε t )φ 1 (x)φ 2 (p) dx dp - R 2n E ε (U b , ρ ε t )φ 1 (x)φ 2 (p)e -|p| 2 λ dx dp ≤ φ 1 1 ∇U b ∞ sup y∈R n |y|| φ2 (y) -φ2 * G (n) λ (y)|
we recognize the first error term in (4.39). So we have only to estimate sup ε∈(0,1)

sup t∈[0,T ] R 2n E ε (U b , ρ ε t )φ 1 (x)φ 2 (p)e -|p| 2 λ dx dp . (4.42) Observe that R 2n E ε (U b , ρ ε t )φ 1 (x)φ 2 (p)e -|p| 2 λ dx dp = I ε,t + II ε,t -III ε,t ,
where

I ε,t := R 2n A(x + ε 2 y) -A(x -ε 2 y) ε φ2 * G (n) λ (y)ρ ε t (x + ε 2 y, x - ε 2 y) dxdyd, (4.43)
II ε,t := R 2n U b (x + ε 2 y) φ 1 (x) -φ 1 (x + ε 2 y) ε φ2 * G (n) λ (y)ρ ε t (x + ε 2 y, x - ε 2 y) dxdy, (4.44)
III ε,t := - R 2n U b (x - ε 2 y) φ 1 (x) -φ 1 (x -ε 2 y) ε φ2 * G (n) λ (y)ρ ε t (x + ε 2 y, x - ε 2
sup t∈[0,T ] |II ε,t | + |III ε,t | ≤ U b ∞ ∇φ 1 ∞ R n |y|| φ2 * G (n) λ |(y) dy.
The estimate of I ε,t is more delicate: we first perform some manipulations of this expression,then we estimate the resulting terms with the help of (4.13).

We expand the convolution product and make the change of variables

u = x + ε 2 y v = x - ε 2 y
to get

I ε,t = 1 (πλ) n/2 ε n R 3n dudvdz A(u) -A(v) ǫ e -|εz-(u-v)| 2 ε 2 λ ρ ε t (u, v) φ2 (z) = 1 ε j∈N µ (ε) j R 2n (Aφ (ε) j,t ) * G (n) λε 2 (v + εz)φ (ε) j,t (v) φ2 (z) dvdz - 1 ε j∈N µ (ε) j R 2n A(v)(φ (ε) j,t * G (n) λε 2 )(v + εz)φ (ε) j,t (v) φ2 (z) dvdz = 1 ε j∈N µ (ε) j R 2n (Aφ (ε) j,t ) * G (n) λε 2 (v + εz) -A(v + εz)(φ (ε) j,t * G (n) λε 2 )(v + εz) φ (ε) j,t (v) φ2 (z) dvdz + 1 ε j∈N µ (ε) j R 2n [A(v + εz) -A(v)] (φ (ε) j,t * G (n) λε 2 )(v + εz)φ (ε) j,t (v) φ2 (z) dvdz. ( 4
µ (ε) j R 2n (Aφ (ε) j,t ) * G (n) λε 2 (v + εz) -A(v + εz)(φ (ε) j,t * G (n) λε 2 )(v + εz) φ (ε) j,t (v) φ2 (z) dvdz = R n dz φ2 (z) R 2n dudv A(v + εz -u) -A(v + εz) ε G (n) λε 2 (u)φ (ε) j,t (v)φ (ε) j,t (v + εz -u) ≤ ∇A ∞ R n dz | φ2 (z)| R 2n dudv |u| ε G (n) λε 2 (u)|φ (ε) j,t (v)||φ (ε) j,t (v + εz -u)| ≤ √ λ ∇A ∞ φ2 1 R 2n |u|G (n) 1 (u) du.
For the second summand in (4.46), using (4.13), we have

1 ε j∈N µ (ε) j R 2n [A(v + εz) -A(v)] (φ (ε) j,t * G (n) λε 2 )(v + εz)φ (ε) j,t (v) φ2 (z) dvdz ≤ j∈N µ (ε) j R 2n A(v + εz) -A(v) ε |(φ (ε) j,t * G (n) λε 2 )(v + εz)||φ (ε) j,t (v)|| φ2 (z)| dvdz ≤ C (2πλ) n/2 ∇A 2 R n |z|| φ2 (z)| dz.
This completes the estimate of the term in (4.43) and the proof.

Verification of (4.32)

This is easy, taking into account the fact that

R 2n W ε ρ t ε * Ḡ(2n) ε (x, p) • ∇ x φ(x, p) dx dp = R 2n W ε ρ t ε ∇ x • [φ * Ḡ(2n) ε ] dx
d P (W (ε k ) t , W t ) = 0. ( 4

.47)

We now prove the following assertions:

(i) W : [0, T ] → P R 2n is weakly continuous and, for all t ∈ [0, T ],

W t = Wt L 2n for some function Wt ∈ L 1 (R 2n ) ∩ L ∞ (R 2n ). Moreover Wt ≥ 0 and sup t∈[0,T ] Wt L 1 (R 2n ) + Wt L ∞ (R 2n ) ≤ C. In particular, W ∈ L ∞ + ([0, T ]; L 1 (R 2n ) ∩ L ∞ (R 2n )). (ii) b ∈ L 1 loc (0, T ) × R 2n
; dt dW t , so the continuity equation (2.8) with ω t = Wt makes sense;

(iii) W solves (2.8) in the sense of distributions on [0, T ] × R 2n ; (iv) For any φ ∈ C ∞ c (R 2n ), t → R 2n φ dW t belongs to C 1 ([0, T ]).
Proof of (i): Observe that (4.47) implies that W : [0, T ] → P R 2n is weakly continuous because it is uniform limit of the weakly continuous maps W (ε k ) . The second part of the proposition follows immediately from (4.11). Indeed, for all φ ∈ L 1 (R 2n ), sup ε∈(0,1) Proof of (iii): First we prove that W solves (2.8) in R 2n \ (S × R n ), where S is the singular set of U s defined in (4.27). Unfortunately this does not follow immediately by (4.30) because we have no information about the singular set Σ of ∇U b , so we cannot control the limit k → ∞ of A Notations and some notions about density operators Observe that the trace condition on ρ can be expressed as follows in terms of its kernel tr(ρ) =

sup t∈[0,T ] R 2n φ(x, p) Wε ρ ε t (x, p) dx dp ≤ C (2π) n R 2n φ(x,
T 0 dt ϕ(t) R 2n dx dp ∇U b (x) • ∇ p φ(x, p) Wε k ρ t ε k (x,

R

  = {y = (y 1 , . . . , y k ) ∈ R k : |y j | ≤ R, j = 1, . . . , k}. We want to prove that lim

  1 and choosing a sequence λ k in Lemma 4.1 converging slowly to 0 for k → +∞, in such a way that ∇A k 2 = o(λ n/4 k ) for k → +∞. In this way we obtain lim k→∞ sup ε∈(0,1)

  (4.50) Observe that (4.50) implies that W t (S × R n ) = 0 for every t ∈ (0, T ). The proof of the global validity of the continuity equation uses the classical argument of removing the singularity by multiplying any test function φ ∈ C ∞ c (R 2n ) by χ k , where χ k (x) = χ(kdist(x, S)) and χ is a smooth cut-off function equal to 0 on [0, 1] and equal to 1 on [2, +∞), with 0 ≤ χ ′ ≤ 2. If we use φχ k as a test function, since χ k depends on x only, we can use the particular structure of b, namely b(x, p) = (p, -∇U (x)), to write the term depending on the derivatives of χ k as k R 2n φχ ′ (kdist(x, S))p • ∇dist(x, S) dW t (x, p)dt. If K is the support of φ, the integral above can be bounded by 2 sup K |pφ| {x∈K:kdist(x,S)≤2}

A density operator on L 2 (

 2 R n ) is a positive, self-adjoint, trace-class operator, namely ρ = ρ * , ρ ≥ 0 and tr(ρ) = 1, where the trace is defined as follows: tr(ρ) := j∈N ϕ j , ρϕ j (A.1)with {ϕ j } j∈N is any orthonormal basis of L 2 (R n ). It can be shown that each density operator ρ is a compact operator, so it can be decomposed as followsρ = j∈N λ j ψ j , • ψ j (A.2)where 0 ≤ λ 1 ≤ λ 2 ≤ • • • ≤ 1, and {ψ j } j∈N is a orthonormal basis of eigenvectors of ρ. Therefore ρ is an integral operator and its kernel is ρ(x, y) = j∈N λ j ψ j (x)ψ j (y), so that ρψ(x) = R n ρ(x, y)ψ(y) dy.

  R n ρ(x, x) dx = 1. (A.3)

  y) dy (4.41) where A := U b φ 1 and C is the constant in (2.11).

	Proof. Set φ2 = F p φ 2 . Observe that since (4.38) gives that
	sup
	ε∈(0,1)

  The estimate b ∈ L 1 loc (0, T ) × R 2n ; dW t dt follows easily from (4.29) and (4.18).

					p) dx dp	(4.48)
	and so	sup t∈[0,T ] R 2n	φ(x, p)dW t (x, p) ≤	C (2π) n R 2n	φ(x, p) dx dp.	(4.49)
	Proof of (ii):

  p) in (2.8) with (4.47). But we can proceed by a density argument because, using the regularity conditions (4.48) and (4.49), we can approximate ∇U b in L 1 on suppφ by bounded continuous functions.In order to prove that W solves (2.8) in [0, T ] × R 2n we use (4.29) to obtain that sup t∈[0,T ] R 2n 1 dist(x, S) 2 dW t (x, p) dt < +∞.

  k dW t (x, p) dt ≤ 8 max K |pφ| k Since the distributional derivative of t → R 2n φW t dx dp is given by R 2n b • ∇φ dW t , we have to show that the map is continuous. Observing that the map t → W t is weakly continuous andW t = Wt L 2n with W ∈ L ∞ + ([0, T ]; L 1 (R 2n ) ∩ L ∞ (R 2n)), the only delicate term is Define the nonnegative Hamiltonian function H = |p| 2 /2+U + U b ∞ . Taking the limit in (4.29) as ε → 0 we easily deduce that supt∈[0,T ] R 2n H 2 dW t ≤ C sup t∈[0,T ] R 2n 1 + |p| 4 + U 2 s (x) dW t < +∞.Since the Hamiltonian is preserved by the Liouville dynamics (under our assumptions on the potential, this fact is contained in the proof of [1, Theorem 6.1]), the above bound implies Hence, if we define the sets A N := {U s ≤ N }, the functionst → f N (t) := A N ∇U s (x) • ∇ p φ dW tare continuous and converge uniformly to R 2n ∇U s (x) • ∇ p φ dW t as N → ∞. This proves (iv).To conclude the proof of the theorem, recalling that W denote the unique distributional solution of(2.8) in L ∞ + ([0, T ]; L 1 (R 2n ) ∩ L ∞ (R2n)) starting from ωL 2n (see[START_REF] Ambrosio | Semiclassical limit of quantum dynamics with rough potentials and well posedness of transport equations with measure initial data[END_REF] Theorem 6.1]), we have proved W = W , and solim P ( Wε k ρ t ε k L 2n , W t L 2n ) = 0. (4.51)Since the limit W t L 2n is independent of the chosen subsequence, this implies the convergence of the whole family, namely

	k→∞	sup t∈[0,T ]
	lim ε→0	sup t∈[0,T ]	d P ( Wε ρ t ε L 2n , W t L 2n ) = 0,	(4.52)
	as desired.
				t →	R 2n	b • ∇φ dW t
	R 2n ∇U sup t∈[0,T ] {H≥N } H 2 dW t =	{H≥N }	H 2 dW 0 → 0	as N → ∞.
	As U s ≤ H, this implies
	sup t∈[0,T ] {Us≥N }	U 2 s dW t ≤

K 1 dist 2 (x, S) dW t (x, p),

and the right hand side is infinitesimal (uniformly in t) as k → ∞.

Proof of (iv): s (x) • ∇ p φ(x, p) dW t .

{H≥N } H 2 dW 0 → 0 as N → ∞. d

The Wigner transform of ρ is defined as

y)e -ipy dy, and the Husimi transform of ρ as

It is easy to check that the marginals of W ε ρ are

where

Similarly the marginals of Wε ρ are

Moreover, the Husimi transform is nonnegative: indeed (see for instance [START_REF] Lions | Sur les mesures de Wigner[END_REF]),

where •, • is the scalar product on L 2 (R n ) and

Hence Wε ψ ≥ 0, and using the spectral decomposition (A.2) one obtains the non-negativity of Wε ρ for any trace-class operator ρ. Moreover, combining (A.3) and (A.6), it follows that Wε ρ is a probability measure.


Stéphane Caro
email: stephane.caro@irccyn.ec-nantes.fr

Guillaume Moroz
email: guillaume.moroz@irccyn.ec-nantes.fr

Thibault Gayral
email: thibault.gayral@ens-cachan.org

Damien Chablat
email: damien.chablat@irccyn.ec-nantes.fr

Chao Chen
email: chao.chen@monash.edu

Singularity Analysis of a Six-dof Parallel Manipulator using Grassmann-Cayley Algebra and Gröbner Bases

published or not. The documents may come

Introduction

The challenge in Virtual Reality (VR) today is to develop "transparent" haptic interface. VR aims at creating high-fidelity reproductions of a real environment, where a user perceives the illusion of being in this environment and interacts with it to perform specific tasks in an effective and comfortable way [START_REF] Alonso | Stepping into Virtual Reality[END_REF]. Haptics is critical in VR applications, where the user actively manipulates the simulated world, such as in surgical simulators [START_REF] Burdea | Force and Touch Feedback for Virtual Reality[END_REF]. There are two main classes of haptic devices: admittance devices and impedance devices [START_REF] Siciliano | Handbook of Robotics[END_REF]. Admittance devices sense the force applied by the operator and then control the operator's position. They are often more expensive because of the cost of the required force sensors, and more complex because of the weight of the force sensors. Impedance devices, also termed force feedback devices, are simpler, cheaper and more common. The fidelity in VR requires light and stiff interfaces. In this vein, Chen et al. developed a new six-dof three-legged parallel manipulator for force-feedback interface [START_REF] Chen | A new 6-dof 3-legged parallel mechanism for forcefeedback interface[END_REF]. The advantage of the proposed device is that all motors are mounted on the ground. This feature significantly reduces the overall mass and inertia of the moving parts, which will eventually yield a more transparent haptic device.

Here, we focus on the singularity analysis of the foregoing manipulator by using Grassmann-Cayley algebra (GCA) and Gröbner bases (GB). Singularities can be found using either numeric, symbolic or geometric methods. For Gough-Stewart parallel manipulators, the rows of the backward Jacobian matrix are the Plücker coordinate vectors of six finite lines that are the lines of action of six actuation forces applied by the actuators to the moving platform. The parallel singularities of such manipulators that Merlet [START_REF] Merlet | Singular Configurations of Parallel Manipulators and Grassmann Geometry[END_REF] analyzed using Grassmann line geometry occur when those lines become linearly dependent. Ben-Horin and Shoham analyzed the parallel singularities of six-dof parallel manipulators using GCA [START_REF] Ben-Horin | Singularity Condition of Six-Degree-Of-Freedom Three-Legged Parallel Robots Based on Grassmann-Cayley Algebra[END_REF][START_REF] Ben-Horin | Application of Grassmann-Cayley Algebra to Geometrical Interpretation of Parallel Robot Singularities[END_REF]. Kanaan et al. [START_REF] Kanaan | Singularity Analysis of Lower-Mobility Manipulators Using Grassmann-Cayley Algebra[END_REF] enlarged the application of GCA to lower-mobility manipulators, in which the rows of the backward Jacobian are either finite lines (zero pitch wrenches) or infinite lines (infinite pitch wrenches). These wrenches, also known as governing lines, are actuation and constraint wrenches applied to the moving platform.

This paper is organized as follows. Section 2 describes the manipulator under study. Sections 3 and 4 are devoted to its singularity analysis with Grassmann-Cayley algebra and Gröbner bases, respectively. Finally, the parallel singularities of the manipulator are plotted in its orientation workspace.

Manipulator Description

Figure 1 illustrates the parallel manipulator under study, which is a simplified kinematic version of the manipulator proposed in [START_REF] Chen | A new 6-dof 3-legged parallel mechanism for forcefeedback interface[END_REF]. It is composed of an equilateral moving platform connected to the base by means of three identical legs. Each leg is composed of three orthogonal prismatic joints and one spherical joint, the first two prismatic joints being actuated. P stands for a prismatic joint whereas S stands for a spherical joint. An underline letter denotes an actuated joint. As a consequence, the manipulator is named 3-PPPS manipulator and provides six-degree-of-freedom motions, i.e., three translations and three rotations.

Parameterization

Let C 1 , C 2 and C 3 be the corners of the moving platform (MP) of side length r. Let F p (C p , X p ,Y p , Z p ) be the frame attached to the moving platform, its origin C p being the centroid of the MP. Y p is parallel to line (C 2 C 3 ) and Z p is normal to the MP. Accordingly,

c 1p =   2r √ 3/6 0 0   , c 2p =   -r √ 3/6 r/2 0   , c 3p =   -r √ 3/6 -r/2 0   (1) O X Y Y Y Z u v w A 1 A 2 A 3 C 1 C 1 C 2 C 2 C 3 C 3 ε1 ∞1 ε1 ∞2 ε1 ∞3
a 1b =   x 1 y 1 0   , a 2b =   0 y 2 z 2   , a 3b =   x 3 0 z 3   (2)
be the Cartesian coordinate vectors of points A 1 , A 2 and A 3 .

Orientation Space

The orientation space can be fully represented with the variables

(Q 2 , Q 3 , Q 4 )
, a subset of the quaternions coordinates. Indeed, the quaternions represent the rotations of the platform with a rotation axis u and an angle θ . The relation between the quaternions and the axis and angle representation can be found in [START_REF] Khalil | Modeling, identification and control of robots[END_REF]:

Q 1 = cos(θ /2), Q 2 = u x sin(θ /2), Q 3 = u y sin(θ /2), Q 4 = u z sin(θ /2) (3)
where u 2 x + u 2 y + u 2 z = 1 and 0 ≤ θ ≤ π. Thus, each rotation can be mapped onto a point of the unit ball in the space defined by the variables (Q 2 , Q 3 , Q 4 ) with the following bijection:

Φ : S ×]0, π] → B \ {   0 0 0   }   u x u y u z   , θ →   Q 2 := u x sin(θ /2) Q 3 := u y sin(θ /2) Q 4 := u z sin(θ /2)   ( 4
)
where S is the unit sphere in a 3-dimension space, and B is the closed unit ball. When θ is equal to zero, the corresponding rotation matrix is the identity that does not depend on the rotation axis u. It also maps to the center of B in the quaternions representation.

Geometric Model

Let c b = c x c y c z T be the Cartesian coordinate vector of point C, the centroid of the MP, expressed in F b . The following equations characterize the geometric model of the 3-PPPS manipulator:

c x -1/3 √ 3Q 1 2 -1/3 √ 3Q 2 2 + 1/6 √ 3 -Q 2 Q 3 + Q 1 Q 4 -x 1 = 0 (5a) c y -1/3 √ 3Q 2 Q 3 -1/3 √ 3Q 1 Q 4 -Q 1 2 -Q 3 2 + 1/2 -y 1 = 0 (5b) c y -1/3 √ 3Q 2 Q 3 -1/3 √ 3Q 1 Q 4 + Q 1 2 + Q 3 2 -1/2 -y 2 = 0 (5c) c z -1/3 √ 3Q 2 Q 4 + 1/3 √ 3Q 1 Q 3 + Q 3 Q 4 + Q 1 Q 2 -z 2 = 0 (5d) c x + 2/3 √ 3Q 1 2 + 2/3 √ 3Q 2 2 -1/3 √ 3 -x 3 = 0 (5e) c z + 2/3 √ 3Q 2 Q 4 -2/3 √ 3Q 1 Q 3 -z 3 = 0 (5f) Q 1 2 + Q 2 2 + Q 3 2 + Q 4 2 -1 = 0 (5g)
3 Singularity Analysis with Grassmann-Cayley Algebra

Grassmann-Cayley Algebra

The Grassmann-Cayley Algebra (GCA), also known as exterior algebra, was developed by H. Grassmann as a calculus for linear varieties operating on extensors with the join and meet operators. The latter are associated with the span and intersec-tion of vector spaces of extensors. Extensors are symbolically denoted by Plücker coordinates of lines and characterized by their step. In the four-dimensional vector space V associated with the three-dimensional projective space P 3 , extensors of step 1, 2 and 3 represent points, lines and planes, respectively. They are also associated with subspaces of V , of dimension 1, 2 and 3, respectively. Points are represented with their homogeneous coordinates, while lines and planes are represented with their Plücker coordinates. The notion of extensor makes it possible to work at the symbolic level and therefore, to produce coordinate-free algebraic expressions for the geometric singularity conditions of spatial parallel manipulators (PMs). For further details GCA, the reader is referred to [START_REF] Ben-Horin | Singularity Condition of Six-Degree-Of-Freedom Three-Legged Parallel Robots Based on Grassmann-Cayley Algebra[END_REF][START_REF] White | The Bracket Ring of a Combinatorial Geometry I[END_REF][START_REF] White | The bracket of 2-extensors[END_REF][START_REF] White | Grassmann-Cayley Algebra and Robotics Applications[END_REF][START_REF] Amine | Singularity Analysis of Lower-Mobility Parallel Robots with an Articulated Nacelle[END_REF].

Twist System of the 3-PPPS Manipulator

A unit screw is given by $ = [s, (s 0 × s + λ s)] T where s is a unit vector along the screw axis, s 0 is the position vector of a point on the screw axis, with respect to a reference frame and λ is the pitch of the screw. A screw of intensify ρ is written as: $ = ρ $. A zero pitch screw $0 = (s, r × s) T (λ = 0) corresponds to the Plücker coordinate vector of a finite line in P 3 . An infinite pitch screw $∞ = (0, s) T (λ → ∞) corresponds to the Plücker coordinate vector of an infinite line in the projective space P 3 .

A twist is a screw representing the instantaneous motion of a rigid body, a wrench is a screw representing a system of forces and moments acting on a rigid body. Let ε denote a twist and τ denote a wrench.

(a) An infinite pitch twist ε ∞ represents a pure translation; (b) A zero pitch twist ε 0 represents a pure rotation; (c) A pure force constrains the translation along its line of action and is represented by a zero pitch wrench τ 0 ; (d) A pure moment constrains the rotation about its direction and is represented by an infinite pitch wrench τ ∞ .

Each leg of the 3-PPPS manipulator provides three independent translations and three independent rotations, represented with six independent twists εi ∞1 , εi ∞2 , εi ∞3 , εi 01 , εi 02 , εi 03 , i = 1, 2, 3, that span its twist sytem T i . As shown in Fig. 1, the independent twists can be defined as:

εi 01 = u c i × u , εi 02 = v c i × v , εi 03 = w c i × w (6)
and

ε1 ∞2 = ε2 ∞3 = ε3 ∞1 = [0, u] T (7a) ε1 ∞1 = ε2 ∞2 = ε3 ∞3 = [0, v] T (7b) ε1 ∞3 = ε2 ∞1 = ε3 ∞2 = [0, w] T (7c)
where c 1 , c 2 and c 3 are the Cartesian coordinate vectors of points C 1 , C 2 and C 3 , respectively; u, v and w are the unit vectors along the X, Y and Z axes of F b .

T = 3 i=1 T i (8)

Wrench System of the 3-PPPS Manipulator

The actuated joints of the 3-PPPS manipulator are the first two prismatic joints of each leg. The actuation wrench τi 01 corresponding to the first prismatic joint of the ith leg is reciprocal to all the twists in T i , but to εi ∞1 . Likewise, the actuation wrench τi 02 corresponding to the second prismatic joint of the ith leg is reciprocal to all the twists in T i , but to εi ∞2 . As a result,

τ1 01 = v c 1 × v , τ2 01 = w c 2 × w , τ3 01 = u c 3 × u (9)
and

τ1 02 = u c 1 × u , τ2 02 = v c 2 × v , τ3 02 = w c 3 × w (10)
In a non-singular configuration, the six actuation wrenches τ1 01 , τ1 02 , τ2 01 , τ2 02 , τ3 01 and τ3 02 span the actuation wrench system of the 3-PPPS manipulator. As the 3-PPPS manipulator does not have any constraint wrench, its global wrench system amounts to its actuation wrench, namely,

W 3-PPPS = span( τ1 01 , τ1 02 , τ2 01 , τ2 02 , τ3 01 , τ3 02 ) (11)
The legs of the 3-PPPS manipulator apply six actuation forces to its movingplatform. Its global wrench system is a six-system. A parallel singularity occurs when the wrenches in the six-system become linearly dependent and span a k-system with k < 6.

Wrench Diagram of the 3-PPPS Manipulator in P 3

The six actuation forces τ1 01 , τ1 02 , τ2 01 , τ2 02 , τ3 01 and τ3 02 form a basis of the global wrench system W 3-PPPS . Those wrenches are represented by six finite lines in P 3 . To obtain the six extensors of the 3-PPPS superbracket, we have to select twelve projective points on the six projective lines, i.e., two points on each line. The extensor of a finite line can be represented by either two distinct finite points or one finite point and one infinite point.

Let c 1 , c 2 , c 3 be the intersection points of τ1 01 and τ1 02 , τ2 01 and τ2 02 , τ3 01 and τ3 02 , respectively. From Fig. 1,τ1 02 and τ3 01 are parallel and intersect at the infinite plane Π ∞ at point x = (x, 0) T , which corresponds to the X direction. Likewise, τ1

01 and τ2 02 are parallel and intersect at the infinite plane Π ∞ at point y = (y, 0) T , which corresponds to the Y direction. Similarly, τ2 01 and τ3 02 are parallel and intersect at the infinite plane Π ∞ at point z = (z, 0) T , which corresponds to the Z direction.

Therefore, the six points have been selected to express the 3-PPPS superbracket τ1 01 ≡ c 1 y, τ1 02 ≡ c 1 x, τ2 01 ≡ c 2 z, τ2 02 ≡ c 2 y, τ3 01 ≡ c 3 x and τ3 02 ≡ c 3 z. The wrench diagram of the 3-PPPS manipulator is shown in Fig. 2.

Superbracket of the 3-PPPS Manipulator

The rows of the backward Jacobian matrix of a parallel manipulator are the Plücker coordinates of six lines in P 3 . The superjoin of these six vectors in P 5 corresponds to the determinant of their six Plücker coordinate vectors up to a scalar multiple, which is the superbracket in GCA Λ (V (2) ) [START_REF] White | Grassmann-Cayley Algebra and Robotics Applications[END_REF]. Thus, a singularity occurs when these six Plücker coordinate vectors are dependent, which is equivalent to a superbracket equal to zero. The expression of the 3-PPPS superbracket is [c 1 x c 1 y c 2 y c 2 z c 3 x c 3 z], which corresponds to six points selected in the robot wrench diagram. This expression can be developed into a linear combination of 24 bracket monomials [START_REF] Ben-Horin | Singularity Condition of Six-Degree-Of-Freedom Three-Legged Parallel Robots Based on Grassmann-Cayley Algebra[END_REF][START_REF] Mcmillan | White The Dotted Straightening Algorithm Journal of Symbolic Computation[END_REF], each one being the product of three brackets of four projective points. The 3-PPPS superbracket was simplified by means of the user interface developed in [START_REF] Ben-Horin | User Interface for the Singularity Analysis of Parallel Robots Based on GrassmannCayley Algebra[END_REF]:

[c 1 x c 1 y c 2 y c 2 z c 3 x c 3 z] = [c 3 x z y] [c 3 c 2 z c 1 ][c 2 x c 1 y] - [c 3 x z c 2 ][c 3 c 2 y c 1 ] (12)

Geometric Condition for the 3-PPPS Manipulator Singularities

Let Π 1 be the plane passing through point C 1 and normal to vector w. Let Π 2 be the plane passing through point C 2 and normal to vector u. Let Π 3 be the plane passing through point C 3 and normal to vector v. Let Π 4 be the plane passing through points C 1 , C 2 and C 3 . From Eq. ( 12) and the user interface developed in [START_REF] Ben-Horin | User Interface for the Singularity Analysis of Parallel Robots Based on GrassmannCayley Algebra[END_REF], it turns out that the 3-PPPS manipulator reaches a singular configuration if and only if planes Π 1 , Π 2 , Π 3 and Π 4 intersect at least at one point.

Singularity Analysis with Gröbner Bases

In this section, we focus on the computation of the parallel singularities of the 3-PPPS manipulator using the Jacobian and Gröbner bases. We derive the algebraic relations of the singularities satisfied by the orientation variables.

Jacobian Formulation

The formula we use to define the parallel singularities is the determinant of a Jacobian matrix. This criterion was introduced in [START_REF] Gosselin | Singularity analysis of closed-loop kinematic chains[END_REF], where parallel singularities were referred to singularities of the first type. Equations (5a)-(g) depend on six joint variables T = (x 1 , y 1 , y 2 , z 2 , x 3 , z 3 ), six pose variables (c x , c y , c z , Q 2 , Q 3 , Q 4 ) and one passive variable (Q 1 ). We denote by X the union of the pose and the passive variables. Let A be the Jacobian matrix of these seven equations with respect to X, i.e.,

A = ∂ F i ∂ X =                         1 0 0 - 2 3 √ 3Q 2 -Q 3 -Q 2 Q 1 - 2 3 √ 3Q 1 + Q 4 0 1 0 - 1 3 √ 3Q 3 - 1 3 √ 3Q 2 -2 Q 3 - 1 3 √ 3Q 1 - 1 3 √ 3Q 4 -2 Q 1 0 1 0 - 1 3 √ 3Q 3 - 1 3 √ 3Q 2 + 2 Q 3 - 1 3 √ 3Q 1 - 1 3 √ 3Q 4 + 2 Q 1 0 0 1 - 1 3 √ 3Q 4 + Q 1 1 3 √ 3Q 1 + Q 4 - 1 3 √ 3Q 2 + Q 3 1 3 √ 3Q 3 + Q 2 1 0 0 4 3 √ 3Q 2 0 0 4 3 √ 3Q 1 0 0 1 2 3 √ 3Q 4 - 2 3 √ 3Q 1 2 3 √ 3Q 2 - 2 3 √ 3Q 3 0 0 0 2 Q 2 2 Q 3 2 Q 4 2 Q 1                         (13)
Moreover, let B be the Jacobian matrix of Eqs. (5a)-(g) with respect to the joint variables. It appears that B is the negative identity matrix. Denoting by F(X, T) the vector of seven polynomials on the left-hand side of Eqs. (5a)-(g), we have:

F(X, T) = 0 (14)
Differentiating Eq. ( 14) with respect to time we obtain:

A Ẋ -Ṫ = 0 (15)
In particular, as in [START_REF] Gosselin | Singularity analysis of closed-loop kinematic chains[END_REF], we can infer from Eq. ( 15) that parallel singularities occur when the determinant of A vanishes:

det(A) = -8 Q 4 √ 3Q 3 3 + 48 Q 2 2 Q 3 Q 1 -48 Q 2 Q 3 2 Q 4 -24 √ 3Q 2 2 Q 3 Q 4 + 48 Q 4 Q 2 Q 1 2 + 24 √ 3Q 2 Q 1 Q 4 2 + 24 Q 3 √ 3Q 1 2 Q 4 -48 Q 3 Q 1 Q 4 2 -24 Q 2 Q 1 √ 3Q 3 2 + 8 √ 3Q 4 3 Q 3 -8 √ 3Q 2 3 Q 1 + 8 Q 2 √ 3Q 1 3 = 0 (16)
Besides, it turns out that the 3-PPPS manipulator does not have any serial singularity as matrix B is always invertible.

Singularities in the Workspace

The singularities of this mechanism can be represented in terms of the pose variables. To this end, we need to eliminate the joint and passive variables from Eq. ( 16) and Eqs. (5a)-(g). This can be achieved with methods based on Gröbner basis theory.

We can notice that these equations depend only on the orientation variables (Q 2 , Q 3 , Q 4 ). This means that the parallel singularities do not depend on the position of the centroid of the moving platform. As a matter of fact, the parallel singularities of the 3-PPPS manipulator can be represented in its orientation workspace only, the latter being characterized with variables (Q 2 , Q 3 , Q 4 ) as shown in Fig. 3.

Q 2 Q 3 Q 4 Q 2 Q 3 Q 4
Fig. 3 The parallel singularities of the 3-PPPS manipulator in the orientation workspace: given a point M, vector --→ OM defines the orientation axis and its Euclidean norm --→ OM is the sine of the half-angle of rotation Finally, note that Eq. ( 18) was also obtained from the geometric condition for the 3-PPPS manipulator singularities given in Sec. 3.6.

Conclusions

The subject of this paper was the singularity analysis of a six-dof three-legged parallel mechanism to be used as a force-feedback interface. A geometric condition for the manipulator singularities was obtained by means of Grassmann-Cayley algebra. The parallel singularities of the manipulator were computed using the Jacobian and Gröbner bases. Consequently, the algebraic relations of the singularities satisfied by the orientation variables were derived. Finally, the parallel singularities of the manipulator were plotted in its orientation workspace.

2 Π 3 Fig. 1

 231 Fig. 1 The 3-PPPS manipulator

Fig. 2

 2 Fig. 2 Wrench diagram in P 3 of the 3-PPPS manipulator

Acknowledgments

This work has been supported by the "Agence Nationale de la Recherche (ANR SIROPA)".

Gröbner basis for elimination

Let P be a set of polynomials in the variables (x 1 , y 1 , y 2 , z 2 , x 3 , z 3 ) and (c x , c y , c z , Q 2 ,Q 3 , Q 4 ). Moreover, let V be the set of common roots of the polynomials in P, and let W be the projection of V on the workspace. It might not be possible to represent W by polynomial equations. Let W be the smallest set defined by polynomial equations that contains W . Our goal is to compute the polynomial equations defining W .

These polynomial equations are computed with Gröbner-basis theory. A Gröbner basis of a polynomial system is a polynomial system equivalent to the first one, and satisfying some additional specific properties. The Gröbner basis of a system depends on an ordering on the monomials. In our case, if we choose an elimination ordering eliminating X, then the Gröbner basis of P will contain exactly the polynomials defining W . This theory is not the subject of this article and shall not be further detailed here. An introduction on elimination with Gröbner basis is available in [START_REF] Cox | Ideals, Varieties, and Algorithms[END_REF]Chapter 3].

Equations of the parallel singularities in the workspace

We can now use the elimination of the previous paragraph for our problem to obtain the polynomial equations defining implicitly the parallel singularities in the workspace. Let us consider the polynomial set:

We compute a Gröbner basis of system [START_REF] Cox | Ideals, Varieties, and Algorithms[END_REF] with respect to elimination ordering eliminating X. This computation yields directly the relation satisfied by the parallel singularities in the orientation workspace, namely,