Calculus of the coarse Ricci curvatures κ and κ on an example of degenerated diffusion

Laurent Veysseire

This work is about studying an example of a degenerated diffusion on which the Ollivier's generalization of coarse Ricci curvature is positive, thus allowing to get a Poincaré inequality.

In section 1, we give the formulas to compute the two different curvatures, κ and κ, introduced in [4]. The curvature κ is the classical Ollivier's coarse Ricci curvature, whereas the curvature κ is a modified version of it, which only exists in special cases and is always smaller than κ. In the example presented in this paper, the curvature κ is less usefull than the curvature κ to get a good constant in the Poincaré inequality, because, due to the symmetries, the curvatures will be the same everywhere on the manifold, so the harmonic mean of κ will always be smaller than the infimum of κ.

In section 2, the example of degenerated diffusion is introduced by defining its generator, along with a lot of notations which will be used in the further sections.

Section 3 describes the eigenvalues and eigenvectors of the generator, and in particular, this gives the true optimal Poincaré constant for the diffusion we are studying here. This section is mostly independent from the other sections and is quite long, so it can be skipped by the reader if he is more interested in the curvature computations than in the spectral analysis of the generator.

Section 4 is where the curvatures κ and κ are actually computed thanks to the formulas given in section 1.

In the section 5 is a discussion about the regularity of the optimal couplings used to define the curvatures κ and κ in this particular example of diffusion, and how it behaves differently than is the non-degenerate case.

The final section 6 shows that the more classical Bakry-Emery curvature approach fails in the example we study here, proving the utility of the Ollivier's coarse Ricci curvature.

This work was supported by the Ecole Normale Superieure de Lyon (France), the mathematics laboratory LAMA (UMR 8050) at the university 1 of Marne-la-Vallée (France) and the TECHNION (Haifa, Israel).

1 Formulas for κ and κ

The coarse Ricci curvatures κ and κ we are computing here are the ones introduced in [4] for diffusions. We first recall the formulas for computing them.

Let (M, g) be a smooth Riemannian manifold, and let

L = 1 2 A ij ∇ i ∇ j + F i ∇ i
be the generator of a diffusion on M, with A a smooth symmetric nonnegative tensor field, F a smooth vector field, and ∇ the Levi-Civita connexion on M.

We take x and y two different points on M , such that (x, y) is not on the cut-locus of M. Then there exist u(x, y) ∈ T x M and u(y, x) ∈ T y M two unit norm vectors such that y = exp x (d(x, y)u(x, y)) and x = exp y (d(x, y)u(y, x)), and they are unique.

M

x y geodesic from x to y

T x M T y M u(x, y) u(y, x)
The distance admits the following Taylor expansion in the neighborhood of (x, y): d(exp x (εv), exp y (εw)) = d(x, y) -ε[g(x)(u(x, y), v) + g(y)(u(y, x), w)] + ε 2 2d(x, y)

[q [START_REF] Bakry | Hypercontractivité de semi-groupes de diffusion[END_REF] (x, y)(v, v) + q (2) (x, y)(w, w) + 2q (12) (x, y)(v, w)] + O(ε 3 ) where q (1) (x, y), q [START_REF] Bakry | Diffusions hypercontractives[END_REF] (x, y) and q (12) (x, y) are quadratic forms, and can be computed by integrating the Jacobi field along the geodesic joining x to y.

Then the coarse Ricci curvature κ between x and y takes the following form: κ(x, y) = 1 d(x, y) (g(x)(u(x, y), F (x)) + g(y)(u(y, x), F (y)))

-1 2 (tr(A(x)q [START_REF] Bakry | Hypercontractivité de semi-groupes de diffusion[END_REF] (x, y)) + tr(A(y)q [START_REF] Bakry | Diffusions hypercontractives[END_REF] (x, y)))

+ tr A(x)q (12) (x, y)A(y)q (12)T (x, y)

The curvature κ is only defined when the tensor field A satisfies the following condition: (H) :

∇ i A jk + ∇ j A ki + ∇ k A ij = 0, with A ij = g ik A kl g lj .
In this case, we take ρ(x, y) := (g(x)A(x)g(x)) (u(x, y), u(x, y)) = (g(y)A(y)g(y)) (u(y, x), u(y, x)).

If ρ(x, y) = 0 we define Ã(x) = A(x) -A(x)g(x)u(x, y)u T (x, y)g(x)A(x) ρ(x, y) and Ã(y) = A(y) -A(y)g(y)u(y, x)u T (y, x)g(y)A(y) ρ(x, y) .

The tensor Ã(x) corresponds to the orthogonal projection of A(x) on Im(A(x))∩ (g(x)u(x, y)) ⊥ , with respect to the scalar product induced by A(x) on Im(A(x)), and the same holds by replacing Ã(x) with Ã(y), A(x) with A(y) and g(x)u(x, y) with g(y)u(y, x). Then the curvature κ can by computed using the following formula:

κ(x, y) = 1 d(x, y) (g(x)(u(x, y), F (x)) + g(y)(u(y, x), F (y))) - 1 2
(tr(A(x)q (1) (x, y))+tr(A(y)q (2) (x, y)))+ q (12) (A(x)g(x)u(x, y), A(y)g(y)u(y, x)) ρ(x, y) + tr Ã(x)q (12) (x, y) Ã(y)q (12)T (x, y) .

(Sorry, I wrote a wrong formula in the first version of [4].)

In the case when ρ(x, y) = 0, the variance term of the distance is 0 for every coupling, so the maximal decay rate of the distance under the constraint that this term cancels (which was the definition of κ) is given by κ, so if ρ(x, y) = 0, we have κ(x, y) = κ(x, y).

Our example of degenerated diffusion

In our example, the manifold M will be the sphere of odd dimension S 2n-1 , and the diffusion will be a Brownian motion restricted to the directions orthogonal to the fibers of the Hopf fibration on S 2n-1 . Due to the symmetries, this diffusion is reversible with respect to the usual normalised volume measure. We will denote by σ 2n-1 this measure.

So we will assume that n ≥ 2, because the resulting diffusion in the case n = 1 is a bit too degenerated.

We will work on C n equipped with an hermitian product y, x . Here, we take the convention where this product is semilinear with respect to y and linear with respect to x. We identify S 2n-1 with the unit sphere of C n for the hermitian norm. The Hopf fiber of x ∈ S 2n-1 is {e iθ .x, θ ∈ [0, 2π)}. The quotient space is naturally the complex projective space of dimension n -1.

The hermitian scalar product on C n inducts an Euclidean one, which is its real part, and a symplectic form, which is its imaginary part. We denote by g 0 the symmetric tensor corresponding to the Euclidean norm, and by ξ the antisymmetric one which corresponds to the symplectic form. For any orthonormal C-basis (x 1 , . . . , x n ), we have:

g 0 = n j=1 x * j ⊗ x * j + (ix j ) * ⊗ (ix j ) * and ξ = n j=1 x * j ⊗ (ix j ) * -(ix j ) * ⊗ x * j ,
where x * = g 0 .x is the R-linear form on C n defined by x * (y) = Re( x, y ). Beware, the tensor products are here over R, so for any two non-zero vectors x and y, the four tensors x ⊗ y, x ⊗ iy, ix ⊗ y and ix ⊗ iy form a linearly free system, even if x = y! The metric we put on S 2n-1 is the one induced by the Euclidean norm, and the tensors of the diffusion are:

A(x) = g -1 0 -x ⊗ x -(ix) ⊗ (ix)
and F (x) = 0. In the particular case where n = 2, we can identify S 3 to the set of unit norm quaternions, and we may approximate our diffusion by Markov chains as follows. Let ε > 0, we discretize the time by taking intervals of length ε 2 . At each step N , we randomly choose p N = ±1 and q N = ±1, independently and taking the value +1 with probability [START_REF] Bakry | Hypercontractivité de semi-groupes de diffusion[END_REF] 2 and the value -1 with probability [START_REF] Bakry | Hypercontractivité de semi-groupes de diffusion[END_REF] 2 , and we jump to

X (N +1)ε 2 = e ε(p N j+q N k) X N ε 2 ,
where the exponential function and the product are the ones in the set of quaternions. Our diffusion is the limit of this random walk when ε tends to 0. Then, with this interpretation, we see a simple coupling between pathes: we take one path Z t of the diffusion starting at 1, and for any two starting points x and y, we set X t = Z t .x and Y t = Z t .y. This is an isometric coupling between the pathes starting at x and those ones starting at y. The simple fact that such a coupling exists implies that A satisfies the (H) condition, and that κ(x, y) ≥ 0 (and thus that κ(x, y) ≥ 0, since κ ≥ κ). This remark also holds more generally for diffusions unvariant by right multiplication (which means that P t f (y.x) = P t f x (y), with f x (z) = f (z.x)) on finite dimentional Lie groups equipped with a metric which is unvariant by left multiplication (i.e. d(x, y) = d(z.x, z.y)).

For n > 2, we saw in Example 22 in [4] that the A which satisfy (H) on S 2n-1 could be defined thanks to a four dimensional tensor in T ∈ R 2n C n having the symmetry of a Riemann tensor in the following way:

A(x)(g(x)u, g(x)u) = T (x, u, x, u)
for every u ∈ T x S 2n-1 , that is, for every u g 0 -orthogonal to x. Now let us consider the four dimensional tensor defined by: T (x, y, z, t) = g 0 (x, z)g 0 (y, t)-g 0 (x, t)g 0 (y, z)-1 3 (2ξ(x, y)ξ(z, t)+ξ(x, z)ξ(y, t)+ξ(x, t)ξ(z, y)).

It has the symmetry of a Riemann tensor, i.e.

T (x, y, z, t) = -T (x, y, t, z) = T (z, t, x, y) T (x, y, z, t) + T (x, z, t, y) + T (x, t, y, z) = 0.

Checking this equalities is easy and is left to the reader. Now, we have, for any u g 0 -orthogonal to x, T (x, u, x, u) = g 0 (u, u)-(ξ(x, u)) 2 , and

A(x)(u, u) = g 0 (u, u) -((ix) * .u) 2 , but we have (ix) * .u = Re( ix, u ) = -Im( x, u ) = -ξ(x, u)
so the tensor T defines the tensor field A in the way stated above, thus A satisfies (H).

Eigenvalues and eigenvectors of the generator

A positive lower bound of the curvature κ or κ implies a positive lower bound for the spectral gap of the generator, as seen in [4].

In this section, we give the whole spectrum and a nice basis of eigenvectors of the generator of our diffusion. So we can get its spectral gap, and we will compare it later to the bounds obtained by computing κ and κ.

A simple description of the eigensubspaces

We first remark we can write our generator as a difference or two more simple operators:

L = 1 2 ∆ S - 1 2 ∆ H
where ∆ S is the classical Laplacian operator on the 2n -1-dimensional sphere, and ∆ H is the Laplacian operator on the Hopf fiber (more explicitely,

∆ H f (x) = d 2 dθ 2 θ=0
f (e iθ x)). As we will see below, the operators ∆ S and ∆ H commute, the nontrivial inersections of an eigensubspace of ∆ S and an eigensubspace of ∆ H form an orthogonal decomposition of the space L 2 (σ 2n-1 ) on subspaces of eigenvectors of L.

Let us remind the eigensubspaces and the eigenvalues of the classical Laplacian ∆ on the sphere S N , which can be found in [3], for example. Definition 1 Let N ∈ N, N ≥ 1. We denote by Π d the space of restrictions to S N of polynomials of degree d on N + 1 variables x 0 , . . . , x N . We also set

E 0 = Π 0 and E d = Π d ∩ (Π d-1 ) ⊥ for d ≥ 1
, where the orthogonal of Π d-1 is taken with respect to the scalar product of L 2 (σ N ).

The following proposition gives us the eigenvalues and eigensubspaces of ∆.

Proposition 2

The subspaces E d are the eigensubspaces of ∆, and the corresponding eigenvalues are -d(d + N -1).

Proof : Let f be a smooth function on R N +1 , and x ∈ S N . We want to compute ∆ f (x), where f is the restriction of f to S N . Let e 1 , . . . , e N be an orthonormal basis of T x S N . Then we have:

∆ f (x) = N i=1 d 2 dt 2 | t=0 f (γ i (t))
where γ i (t) = exp x (t.e i ) = cos(t)x + sin(t)e i is the geodesic starting at x with speed e i . We have

d 2 dt 2 f (cos(t)x + sin(t)e i ) = -cos(t)∇ x f (cos(t)x + sin(t)e i ) + sin(t)∇ e i f (cos(t)x + sin(t)e i ) + sin(t) 2 ∇ 2 xx f (cos(t)x + sin(t)e i ) -2 sin(t) cos(t)∇ 2 xe i f (cos(t)x + sin(t)e i ) + cos(t) 2 ∇ 2 e i e i f (cos(t)x + sin(t)e i ),
where ∇ and ∇ 2 are the usual gradient and hessian in R N +1 . So we get:

∆ f (x) = -N ∇ x f (x) + N i=1 ∇ 2 e i e i f (x) = ∆ 0 f (x) -N ∇ x f (x) -∇ 2 xx f (x),
where ∆ 0 is the usual Laplacian on R N +1 , and the last equality is due to the simple fact that ∆

0 f = ∇ 2 xx f + N i=1 ∇ 2 e i e i f , because x, e 1 , . . . , e N is an orthonormal basis of R N +1 .
Let us apply the preceding formula for f being a monomial of degree d:

f (x) = N i=0 x α i i , with N i=0 α i = d. We have ∆ f (x) = N i=0   α i (α i -1)x α i -2 i j∈[|0,N |],j =i x α j j   -N N i=0 x i   α i x α i -1 i j∈[|0,N |],j =i x α j j   - N i=0 N j=0 x i x j   1 i=j α i (α i -1)x α i -2 i k∈[|0,N |],k =i x α k k +1 i =j α i α j x α i -1 i x α j -1 j k∈[|0,N |],k =i,k =j x α k k   =g(x) + f (x)(-N N i=0 α i - N i=0 N j=0 (α i α j -1 i=j α i )) =g(x) -d(d + N -1)f (x) where g = ∆ 0 f is an homogenous polynomial of degree d -2 (or g = 0 if d < 2).
By linearity, we get that Π d is stabilized by ∆, and more precisely, for any f ∈ Π d , we have ∆f

+ d(d + N -1)f ∈ Π d-1 . Let f ∈ E d , we set g = ∆f +d(d+N -1)f ∈ Π d-1 . We have S N ḡ∆f dσ N = S N ḡ(-d(d+N -1)f + g)dσ N = S N |g| 2 dσ N because f ∈ E d , so f is orthogonal to g. But ∆ is a symmetric operator of L 2 (σ N ), so we have S N ḡ∆f dσ N = S N f ∆(ḡ)dσ N = 0 because ∆g ∈ Π d-1 , so it is orthogonal to f . Thus S N |g| 2 dσ N = 0, hence g = 0 and ∆f = -d(d + N -1)f , that is to say, provided f is not 0, it is an eigenvector of ∆ with the eigenvalue -d(d + N -1).
The linear span of the E d is the set of the restrictions of polynomials to S N , which is dense in the space of continuous functions on S N , and thus in L 2 (σ N ), so we have found all the eigenvalues.

One can also compute the dimensions of the spaces E d to get the multiplicity of the eigenvalues.

As Π d-1 ⊂ Π d , we have dim(E d ) = dim(Π d ) - dim(Π d-1
). Two distinct polynomials can have the same restriction on S N , in which case their difference is a polynomial that cancels on S N , so it is a multiple of

P 0 = N i=0 x 2 i -1.
So the kernel of the linear map which associates to a polynomial of degree d or less its restriction to S N is exactly the set of the products of P 0 and a polynomial of degree less or equal to d -2.

The dimension of the space of polynomials of degree d or less in N + 1 variables is the binomial coefficient

N +d+1 N +1 . So the dimension of Π d is N +d+1 N +1 -N +d-1 N +1 = N +d N + N +d-1 N . And thus the dimension of E d is N +d N -N +d-2 N = N +d-1 N -1 + N +d-2 N -1 .
Classical Fourier analysis gives us the eigenvalues and eigenvectors of ∆ H .

Definition 3 Let k ∈ Z, we set F k the set of L 2 functions on S 2n-1 satis- fying f (e iθ x) = e kiθ f (x) for all x ∈ S 2n-1 and θ ∈ [0, 2π). Proposition 4 Any non 0 element of F k is an eigenvector of ∆ H for the eigenvalue -k 2 . The eigensubspaces of ∆ H are F 0 , and the F k ⊕ F -k for k ≥ -1. Proof : The fact that the vectors of F k are eigenvectors of ∆ H is trivial. Let f ∈ L 2 (σ 2n-1
), then, by Fubini, the restriction of f on alomst every Hopf fibers (with respect to the quotient measure of σ 2n-1 on the set of Hopf fibers P C n ) is L 2 . One can use the classical Fourier decomposition on each Hopf fiber for which the restriction of f is L 2 , so using Fubini again, we have f = k∈Z f k , where f k (z) = 2π 0 e -ikθ f (e iθ z)dθ, which is well defined for almost all z, and f k ∈ F k .

Polynomials on the real and imaginary parts x j and y j of the complex coordinates z j (j ∈ [|1, n|]) can be written as polynomials on the z j 's and their complex conjugates z j 's.

Remark 5 For any monomial n j=1 z α j j z j β j , its restriction to S 2n-1 belongs to F k , with k = n j=1 α j -β j .

Definition 6 Let M d,k be the set of monomials n j=1 z

α j j z j β j of degree d = n j=1 α j + β j whose restriction to S 2n-1 belongs to F k .
Well note that M d,k is empty if d and k do not have the same parity or if |k| > d. More precisely, we have

#M d,k = d+k 2 + n -1 n -1 d-k 2 + n -1 n -1 1 d+k∈2Z,|k|≤d ,
since the α j 's and β j 's are nonnegative integers satisfying n j=1 α j = d+k 2 and n j=1

β j = d-k 2 .
Definition 7 We denote by Π d,k the vector space of the restrictions to S 2n-1 of the polynomials which are linear combinations of the monomials in M d,k .

The dimension of Π d,k is #M d,k , because the restriction to S 2n-1 of homogeneous polynomials of degree d is injective. We remark that Π d,k ⊂ Π d+2,k , since when one multiplies a monomial of M d,k by the polynomial n j=1 z j z j , whose restriction to S 2n-1 is 1, every monomial of the product belongs to M d+2,k .

Definition 8 We set E d,k = Π d,k ∩ Π ⊥ d-2,k , where Π d-2,k = {0} if d < 2. The dimension of E d,k is #M d,k -#M d-2,k , which is non-zero if n > 1, k + d ∈ 2Z and |k| ≤ d or n = 1 and |k| = d.
Lemma 9 Let E d be as in Definition 1, with N = 2n -1 Then, we have

E d,k = E d ∩ F k . Proof : We first prove E d,k ⊂ E d ∩ F k . The inclusion E d,k ⊂ F k is trivial, let us prove that E d,k ⊂ E d . Let f ∈ E d,k
, and g ∈ Π d-1 . We trivially have f ∈ Π d , it only remains to check S 2n-1 f g = 0. The function g is the restriction of some polynomial P (z 1 , z 1 , . . . , z n , z n ) of degree smaller than d -1. The restriction of each monomial of P belongs to some F k , with |k | ≤ d -1. So we can write P = d-1 k =-d+1 P k , where each monomial of P k belongs to some M d ,k , with d ≤ d -1. The restriction g k of P k to S 2n-1 belongs to F k (actually, g k is the orthogonal projection of g on F k ). Thus f is orthogonal to all the g k with k = k (remember f ∈ F k ). So it suffices to show that f is orthogonal to

g k . If d + k / ∈ 2Z, then f = 0, and this is trivial. If d + k ∈ 2Z, every monomial of P k lies in some M d ,k , with d ≤ d -1 and d + k ∈ 2Z. Thus g k ∈ d 2 j=1 Π d-2j,k = Π d-2,k , so f is orthogonal to g k , since f ∈ E d,k . So we have just proved the inclusion E d,k ⊂ E d ∩ F k . It remains to prove E d ∩ F k ⊂ E d,k . Let f ∈ E d ∩ F k . Since f ∈ E d ,
there exists a polynomial P of degree at most d such that f is the restriction of P to S 2n-1 . One can use the same decomposition as above P = d k =-d P k . As f ∈ F k , the restriction of P k to S 2n-1 must be 0 for all k = k, so we can assume P = P k . If P has no term of degree d, then f ∈ Π d-1 , but since f ∈ E d , f must be orthogonal to f , and then f = 0, so f ∈ E d,k . If P = P k has at least one term of degree d, then k + d ∈ 2Z, and every monomial of

P k lies in some M d ,k with d + k ∈ 2Z and d ≤ d, thus f ∈ Π d,k . Let g ∈ Π d-2,k ⊂ Π d-2 ⊂ Π d-1 . Since f ∈ E d , f is orthogonal to g, so f ∈ Π d,k ∩ Π ⊥ d-2,k = E d,k . Hence E d ∩ F k ⊂ E d,k
. So the lemma is proved. Now we are able to prove the following theorem, which gives the eigenvalues and eigensubspaces of the operator L:

Theorem 10 For n ≥ 2, the eigenvalues of L are non-positive integers, and more precisely those ones which can be written

λ d,k := - d(d + 2n -2) -k 2 2 ,
where d is a non-negative integer and k an integer such that k + d ∈ 2Z and |k| ≤ d. For a such eigenvalue λ, the corresponding eigensubspace is Proof :If n = 1, L = 0, so it is trivial. If n ≥ 2, we have λ 0,0 = 0 and E 0,0 is the set of constant functions, and so has dimension 1. We have λ 1,1 = λ 1,-1 = -(n -1), and for any d ≥ 2, we have for any

(d,k)|λ d,k =λ E d,
k such that k + d ∈ 2Z and |k| ≤ d λ d,k ≤ -d(n -1) ≤ -2(n -1) < λ 1,1 , so λ 1,1
is the greatest non-zero eigenvalue.

A nice orthonormal eigenvector basis

In this subsection, we give a nice orthonormal basis of E d,k . To check that a given set of vectors of E d,k is an orthonormal basis, it suffices to show they are orthogonal with each other, have norm 1, and that the cardinal of the set equals the (finite) dimension of E d,k .

We introduce some systems of coordinates to carry out the computations. We set

r j := |z j | 2 and θ j := arg(z j ) for j ∈ [|1, n|].
The angular coordinate θ j is not uniquely defined if z j = 0, but this only happens on a set of measure 0 with respect to σ 2n-1 . The coordinates r j are nonnegative and satisfy the relation

n j=1 r j = 1.
So we parametrize S 2n-1 with the 2n-1 coordinates θ 1 , . . . , θ n , r 1 , . . . , r n-1 , where the θ j 's belong to [0, 2π) and the n -1-uple (r 1 , . . . , r n-1 ) belongs to the domain D n of R 2n-1 defined by r j ≥ 0 and n-1 j=1 r j ≤ 1. This domain D n is the closed simplex whose vertices are 0 and the n -1 vectors of the canonical basis of R n-1 .

If we denote by ϕ the parametrization:

ϕ : [0, 2π) n × D n → S 2n-1 (θ 1 , . . . , θ n , r 1 , . . . , r n-1 ) → √ r 1 e iθ 1 , . . . , √ r n-1 e iθ n-1 , 1 -n-1 j=1 r j e iθn
we have the following.

Proposition 12 The pulled back measure ϕ (σ 2n-1 ) of σ 2n-1 by ϕ is the normalized Lebesgue measure on [0, 2π) n × D n .

Proof : We compute the matrix G of the pulled back metric, since the density of the pulled back Riemannian volume measure with respect to the Lebesgue measure is det(G). We have for 1 ≤ j ≤ n ∂ϕ ∂θ j = (0, . . . , 0, i √ r j e iθ j , 0, . . . , 0)

with r n = 1 -n-1 j=1 r j , and for 1 ≤ j ≤ n -1 ∂ϕ ∂r j = (0, . . . , 0, e iθ j 2 √ r j , 0, . . . , 0, -e iθn 2 √ r n ).

We get then

G(θ 1 , . . . , θ n , r 1 , . . . , r n-1 ) = Diag(r 1 , . . . , r n ) 0 0 Diag( 1 r 1 , . . . , 1 r n-1 ) + 1 rn 1 n-1
,

where Diag(a 1 , . . . , a k ) is the diagonal matrix whose diagonal coefficients are a 1 , . . . , a k , and 1 k is the k × k square matrix whose coefficients all equal to 1.

The determinant of G is thus

det(G) = r 1 . . . r n 1 r 1 + 1 rn 1 rn • • • 1 rn 1 rn 1 r 2 + 1 rn . . . . . . . . . . . . . . . 1 rn 1 rn • • • 1 rn 1 r n-1 + 1 rn = r 1 . . . r n 1 r 1 + 1 rn 1 rn • • • 1 rn -1 r 1 . . . -1 r 1 Diag( 1 r 2 , . . . , 1 r n-1 )
,

where we have substracted the first line of the matrix to the other ones to get the last equality. Now we use the formula det(M ) = σ∈S n-1 (σ) n-1 j=1 m j,σ(j) for the determinant of a (n -1) × (n -1) square matrix. The last matrix has a lot of zeros, and to get a non-zero term, the permutation σ must satisfy σ(j) = 1 or j for every 2 ≤ j ≤ n -1, which is only possible if σ is the identity or a transposition (1, j). Thus we get:

det(G) = r 1 . . . r n   ( 1 r 1 + 1 r n ) n-1 j=2 1 r j - n-1 j=2 (- k =j 1 r k )   = r n +r 1 + n-1 j=2 r j = 1.
Hence the pulled back Riemannian volume measure on S 2n-1 by ϕ is the Lebesgue measure on [0, 2π) n × D n . We can remark that the Lebesgue measure of [0, 2π) n × D n is (2π) n (n-1)! , which is indeed the Riemannian volume of S 2n-1 .

We only have to renormalize the measures to get the proposition. The domain D n can be parametrized on a simple way. For 1

≤ j ≤ n-1, we set X j := n k=j+1 r k n k=j r k ∈ [0, 1]
. The X j 's are all well defined as soon as r n-1 and r n are not both 0, which occurs on a subset of Lebesgue measure 0. We have

r j = (1 -X j ) j-1 k=1 X k for 1 ≤ j ≤ n -1. The application ψ : [0, 1] n-1 → D n (X 1 , . . . , X n-1 ) → (1 -X 1 , X 1 (1 -X 2 ), . . . , X 1 . . . X n-2 (1 -X n-1 ))
is a diffeomorphism from (0, 1) n-1 to the interior of D n . The matrix of its differential is lower triangular, and its diagonal coefficients are -1, -X 1 , -X 1 X 2 , . . . , -X 1 . . . X n-2 . So its determinant is (-1) n n-1 j=1 X n-1-j j , and thus we have:

dψ (1 Dn Leb n-1 ) = n-1 j=1 X n-1-j j 1 [0,1] (X j )dX j . (1)
To give the expression of our eigenvectors, we define the following polynomials:

Definition 13 Let a and b be two real numbers greater than -1. For any k ∈ N, we set P a,b,d the polynomial of degree d obtained by performing the Gram-Schmidt orthogonalisation algorithm on the canonical basis of polynomials for the following scalar product

f, g a,b = 1 0 f (x)g(x)x a (1 -x) b dx.
Remark 14 The polynomials P a,b,d have the following expressions:

P a,b,d (x) = (-1) d (a + b + 2d + 1)Γ(a + b + d + 1) d!Γ(a + d + 1)Γ(b + d + 1) 1 x a (1 -x) b d d dx d (x a+d (1 -x) b+d ) =(-1) d (a+b+2d+1)Γ(a+b+d+1) d!Γ(a+d+1)Γ(b+d+1) d k=0 (-1) k d!Γ(a+d+1)Γ(b+d+1) k!(d-k)!Γ(a+k+1)Γ(b+d-k+1) x k (1-x) d-k =(-1) d (a+b+2d+1)Γ(a+b+d+1) d!Γ(a+d+1)Γ(b+d+1) d k=0 (-1) k d!Γ(a+d+1)Γ(a+b+d+k+1) k!(d-k)!Γ(a+k+1)Γ(a+b+d+1) x k = (-1) d P b,a,d (1 -x),
where Γ(x) = +∞ 0 t x-1 e -t dt is the Euler's gamma function.

Proof: Let a and b be two real numbers greater than -1. We first show that the functions Q a,b,d (x) :=

1 x a (1-x) b d d dx d (x a+d (1 -x) b+d
) are polynomials of degree d and are orthogonals to each other for ., . a,b . After that, we will compute

Q a,b,d 2 a,b to renormalize Q a,b,d . We have Q a,b,d (x) = 1 x a (1 -x) b d d dx d (x a+d (1 -x) b+d ) = d k=0 d k 1 x a d d-k dx d-k (x a+d ) 1 (1 -x) b d k dx k ((1 -x) b+d ) = d k=0 (-1) k d k Γ(a + d + 1) Γ(a + k + 1) x k Γ(b + d + 1) Γ(b + d -k + 1) (1 -x) d-k , ( 2
)
so Q a,b,d is a polynomial of degree at most d.

Let 0 ≤ d 1 ≤ d 2 .
We have:

Q a,b,d 1 , Q a,b,d 2 a,b = 1 0 x a (1 -x) b Q a,b,d 1 (x)Q a,b,d 2 (x)dx = 1 0 Q a,b,d 1 (x) d d 2 dx d 2 (x a+d 2 (1 -x) b+d 2 )dx = d 2 k=1 (-1) k-1 Q (k-1)
a,b,d 1 (x) According to equation 2, we have

d d 2 -k dx d 2 -k (x a+d 2 (1 -x) b+d 2 ) 1 x=0 + (-1) d 2 1 0 Q (d 2 ) a,b,d 1 (x)x a+d 2 (1 -x) b+d 2 dx (n -2 integrations by parts) = d 2 k=1 (-1) k-1 Q (k-1) a,b,d 1 (x)x a+k (1 -x) b+k Q a+k,b+k,d 2 -k (x) 1 x=0 + (-1) d 2 1 0 Q (d 2 ) a,b,d 1 (x)x a+d 2 (1 -x) b+d 2 dx = (-1) d 2 1 0 Q (d 2 ) a,b,d 1 (x)x a+d 2 (1 -x) b+d 2 dx, (3)
Q a,b,d (x) = d k=0 (-1) k d k Γ(a + d + 1) Γ(a + k + 1) x k Γ(b + d + 1) Γ(b + d -k + 1) (1 -x) d-k = d k=0 d-k l=0 (-1) k+l d k d -k l Γ(a + d + 1)Γ(b + d + 1) Γ(a + k + 1)Γ(b + d -k + 1) x k+l = d k =0 (-x) k d k Γ(a + d + 1)Γ(b + d + 1) Γ(a + b + d + 1) k k=0 k k Γ(a + b + d + 1) Γ(a + k + 1)Γ(b + d -k + 1) (change of variables k = k + l) = d k =0 (-x) k d k Γ(a + d + 1)Γ(b + d + 1) Γ(a + b + d + 1) Γ(a + b + d + k + 1) Γ(b + d + 1)Γ(a + k + 1) , (4)
where the last equality occurs because of the identity

k k=0 k k Γ(z 1 + 1) Γ(z 2 -k + 1)Γ(z 1 -z 2 + k + 1) = Γ(z 1 + k + 1) Γ(z 2 + 1)Γ(z 1 -z 2 + k + 1)
, which holds for any k ∈ N, and for any z 1 and z 2 in C such that z 1 is not a negative integer. This identity can be easily shown by induction on k , using the fact that

k k = k -1 k-1 + k -1
k , but it also has a combinatoric interpretation when z 1 is a nonnegative integer and z 2 is an integer not less than k . In this case, the identity can be written as

k k=0 k k z 1 z 2 -k = z 1 + k z 2 .
The right hand term is the number of ways to choose z 2 elements among a set of z 1 + k elements, and the generic term of the sum of the left hand term is the number of ways to choose z 2 elements among a set of z 1 + k elements with the constraint that k of them belong to a fixed subset of cardinal k . So according to equation 4, we have lc(Q a,b,d ) = (-1) d Γ(a+b+2d+1) Γ(a+b+d+1) . Thus we have

Q a,b,d 2 a,b = d!Γ(a + b + 1)Γ(b + d + 1) (a + b + 2d + 1)Γ(a + b + d + 1)
.

Hence P a,b,d = ± (a+b+2d+1)Γ(a+b+d+1) d!Γ(a+d+1)Γ(b+d+1) Q a,b,d . Because of its definition, the leading coefficient of P a,b,d is positive, so finally we have

P a,b,d = (-1) d (a + b + 2d + 1)Γ(a + b + d + 1) d!Γ(a + d + 1)Γ(b + d + 1) Q a,b,d .
Replacing Q a,b,d with its definition in this last equality gives us the first equality of Remark 14. For the second one, one uses equation 2, and for the third one, one uses equation 4. The last equality of Remark 14 is a straightforward consequence of the first one or the second one. Now our eigenvectors will be defined on the following way:

Definition 15 Let α 1 , . . . , α n be n integers, and d 1 , . . . , d n-1 be n -1 nonnegative integers. We set

F α 1 ,...,αn,d 1 ,...,d n-1 (θ 1 ,...,θn,X 1 ,...,X n-1 )= 1 √ (n-1)! n j=1 ((1-Xj) j-1 k=1 X k ) |α j | 2 e iα j θ j n-1 j=1 P a j ,|α j |,d j (X j ) j-1 k=1 X d j k ,
where we take X n = 0 in the last term of the first product, and for 1 ≤ j ≤ n -1,

a j = n -1 -j + |α n | + n-1 k=j+1 (|α k | + 2d k ).
The function f α 1 ,...,αn,d 1 ,...,d n-1 will be defined on almost all S 2n-1 by

f α 1 ,...,αn,d 1 ,...,d n-1 (z 1 , . . . , z n ) = F α 1 ,...,αn,d 1 ,...,d n-1 (θ 1 , . . . , θ n , X 1 , . . . , X n-1 ),
with θ j = arg(z j ) and X j = Proof : The fact that f α 1 ,...,αn,d 1 ,...,d n-1 belongs to F α 1 +...+αn is rather trivial, since when one multiplies z by e iθ , one just adds θ (modulo 2π) to each θ j , and so the value of the function is multiplied by e i(α 1 +...+αn)θ .

To prove that f α 1 ,...,αn,d 1 ,...,d n-1 is the restriction to S 2n-1 of a polynomial function, it suffices to check that the j-th term of the first product in the definition of F α 1 ,...,αn,d 1 ,...,d n-1 is nothing but z α j j if α j is nonnegative, and z j -α j else, and that the j-th term of the second product is an homogenous polynomial of degree d j of the |z l | 2 's (and thus an homogenous polynomial of degree 2d j of the z l 's and z l 's). More precisely, according to Remark 14, it is (-1) d j (a j +|α j |+2d j +1)(a j +|α j |+d j )! d j !(a j +d j )!(|α j |+d j )! d j h=0

(-1) h d j !(a j +d j )!(|α j |+d j )! h!(d j -h)!(a j +h)!(|α j |+d j -h)! (|z j+1 | 2 +...+|zn| 2 ) h (|z j | 2 ) d j -h .
Thus the product of these guys is an homogenous polynomial of degree

n j=1 |α j | + 2 n-1 j=1 d j .
Let us check the orthonormality of these functions. According to Proposition 12 and equation 1, we have:

S 2n-1 f α 1 ,...,αn,d 1 ,...,d n-1 f α 1 ,...,α n ,d 1 ,...,d n-1 dσ 2n-1 = (n-1)! F α 1 ,...,αn,d 1 ,...,d n-1 (θ 1 ,...,θn,X 1 ,...X n-1 ) F α 1 ,...,α n ,d 1 ,...,d n-1 (θ 1 ,...,θn,X 1 ,...X n-1 ) X n-2 1 ...X n-2 dX 1 ...dX n-1 dθ 1 2π ... dθn 2π = n j=1 2π 0 e i(α j -α j )θ j dθ j 2π n-1 j=1 1 0 (1-X j ) |α j |+|α j | 2 X n-1-j+ n k=j+1 |α k |+|α k | 2 + n-1 k=j+1 d k +d k j P a j ,|α j |,d j (X j )P a j ,|α j |,d j (X j )dX j .
If the α j = α j for some j, the j-th integral of the first product is 0. Now assume that α j = α j for all 1 ≤ j ≤ n. If for some j, we have d j = d j , let's assume that j is the largest index satisfying this inequality. Then the j-th integral of the second product is

1 0 (1 -X j ) |α j | X a j j P a j ,|α j |,d j (X j )P a j ,|α j |,d j (X j ),
and this integral is 0 because of the orthogonality of P a j ,|α j |,d j and P a j ,|α j |,d j for the scalar product ., . a j ,|α j | . Finally, if α j = α j and d j = d j , all the integrals are 1, so our functions have norm 1.

It remains to prove that f α 1 ,...,αn,d 1 ,...,d n-1 is orthogonal to every polynomial of degree less than d = n j=1 |α j | + 2 n-1 j=1 d j , and that there are For the second fact, we associate to each (2n-1)-uple (α 1 , . . . , α n , d 1 , . . . , d n-1 ) that satisfies n j=1 |α j | + 2 n-1 j=1 d j = d and n-1 j=1 α j = k the monomial n j=1 z α j j z j β j , where for each 1 ≤ j ≤ n we take α j = α j + d j and β j = d j if α j ≥ 0, and α j = d j and β j = -α j + d j if α j < 0, with d n = 0. This is a bijection between our set of (2n -1)-uples and the monomials in M d,k which satisfy α n = 0 or β n = 0. The monomials in M d,k which satisfy α n > 0 and β n > 0 are just z n z n times a monomial of M d-2,k . So the cardinal of our set of (2n -1)-uples is #M )) when the α j 's are not all 0.

dim(E d,k ) (2n -1)-uples (α 1 , . . . , α n , d 1 , . . . , d n-1 ) satisfying n j=1 |α j | + 2 n-1 j=1 d j = d and n-1 j=1 α j = k.

Calculus of κ and κ

Since our diffusion is degenerated, we can't use the formulas given in [4] for the close-points limit curvatures κ(x, u) and κ(x, u). We will compute κ(x, y) and κ(x, y) and take the limit when d(x, y) → 0. Now we assume that x and y are fixed on S 2n-1 , with y = ±x. Then there exists an unit norm vector z ∈ C n , C-orthogonal to x, such that u(x, y) = cos(α)ix + sin(α)z, with α the angle between ix and u(x, y) (z is unique if α / ∈ {0, π}). We choose e 3 , e 4 , . . . , e n an orthonormal C-basis of the C-orthogonal of Vect C (x, z). The orthonormal basis of T x S 2n-1 we choose for our computations will be B(x) = (ix, z, iz, e 3 , ie 3 , e 4 , ie 4 , . . . , e n , ie n ).

If we note d = d(x, y) ∈ (0, π), we have y = cos(d)x + sin(d)u(x, y) and u(y, x) = -cos(d)u(x, y) + sin(d)x. Then we get u(y, x) = -cos(α)iysin(α)z , with

z = cos(d)z -sin(d) cos(α)iz -sin(d) sin(α)x.
z is C-orthogonal to y, its norm is 1. Then Vect C (y, z ) = Vect C (x, z), because the first one is clearly included in the second one and they have both dimension 2. So the following basis of T y S 2n-1 is orthonormal: B(y) = (iy, z , iz , e 3 , ie 3 , . . . , e n , ie n ).

The matrix of A(x) in the basis B(x) is the same that the matrix of A(y) in B(y). They both are

A =     0 0 0 0 1 0 0 0 1 0 0 I 2n-4     .
It remains to compute the matrices of q [START_REF] Bakry | Hypercontractivité de semi-groupes de diffusion[END_REF] , q [START_REF] Bakry | Diffusions hypercontractives[END_REF] and q (12) . For v ∈ T x S 2n-1 and w ∈ T y S 2n-1 , we have

exp x (εv) = 1 - ε 2 v 2 2 x+εv+O(ε 3 ) and exp y (εw) = 1 - ε 2 w 2 2 y+εw+O(ε 3 ).
Then, we have We identify the second order terms after simplification, and get:

d(exp x (εv), exp y (εw)) = arccos g 0 1 - ε 2 v 2 2 x + εv, 1 - ε 2 w 2 2 y + εw + O(ε 3 ) = arccos(cos(d)+ε(g 0 (x, w)+g 0 (v, y))+ ε 2 2 (2g 0 (v, w)-( v 2 + w 2 ) cos(d))+O(ε 3 )) = d-ε g 0 (x, w) + g 0 (y, v) sin(d) + ε 2 2 ( v 2
q (1) (v, v) = cos(d) d sin(d) ( v 2 -g 0 (v, u(x, y)) 2 ), q (2) (w, w) = cos(d) d sin(d) ( w 2 -(g 0 (w, u(y, x))) 2 )
and q (12) (v, w) = 1 d sin(d) (-g 0 (v, w) -cos(d)g 0 (v, u(x, y))g 0 (w, u(y, x))).

To get the matrices of q [START_REF] Bakry | Hypercontractivité de semi-groupes de diffusion[END_REF] , q [START_REF] Bakry | Diffusions hypercontractives[END_REF] and q [START_REF] Bakry | Diffusions hypercontractives[END_REF] in the bases B(x) and B(y), we only have to apply these quadratic forms to the vectors of the bases and so we get the coefficients.

The matrix of q [START_REF] Bakry | Hypercontractivité de semi-groupes de diffusion[END_REF] in the basis B(x) is the same that the matrix of q [START_REF] Bakry | Diffusions hypercontractives[END_REF] in the basis B(y) and equals to cos(d) d sin(d) Q, where

Q =     sin 2 (α) -sin(α) cos(α) 0 -sin(α) cos(α) cos 2 (α) 0 0 0 1 0 0 I 2n-4     .
The matrix of q (12) in the bases B(x) and B(y

) is 1 d sin(d) Q 12 , with Q 12 =     -cos(d) sin 2 (α) cos(d) sin(α) cos(α) sin(d) sin(α) cos(d) sin(α) cos(α) -cos(d) cos 2 (α) -sin(d) cos(α) -sin(d) sin(α) sin(d) cos(α) -cos(d) 0 0 -I 2n-4     .
To compute κ, we also need the expressions of Ã(x) and Ã(y). Here, we have A(x)g(x)u(x, y) = sin(α)z, A(y)g(y)u(y, x) = -sin(α)z and thus ρ(x, y) = sin 2 (α). So we need to compute κ only if α / ∈ {0, π}, in which case, the matrices of Ã(x) and Ã(y) in the bases B(x) and B(y) both equal to the very nice

à =     0 0 0 0 0 0 0 0 1 0 0 I 2n-4     .
Then, according to our formulas, we have in the case when sin(α) = 0, i.e in the case when x and y are not on the same Hopf fiber. We have tr(AQ) = cos 2 (α) + 2n -3, and q (12) (z, z ) = -cos(d) cos 2 (α)

κ(x, y) = 0 - cos (
d sin(d)
. The eigenvalues of the matrix AQ 12 AQ T 12 are 0, r 1 , r 2 and 1 with multiplicity 2n -4, with Using a simple trick, we can simplify √ r 1 + √ r 2 : we have (

√ r 1 + √ r 2 ) 2 = r 1 + r 2 + 2 √ r 1 r 2 . Thus, we get √ r 1 + √ r 2 = r 1 + r 2 + 2 √ r 1 r 2 = cos 2 (d) sin 4 (α) + 2 cos 2 (α) + 2 cos 4 (α) = cos 2 (d) sin 4 (α) + 4 cos 2 (α).
The eigenvalues of ÃQ 12 ÃQ T 12 are 0 with multiplicity 2, cos 2 (d) and 1 with multiplicity 2n -4.

So finally we get the following expressions of κ and κ:

κ(x, y) = cos 2 (d) sin 4 (α) + 4 cos 2 (α) + 2n -4 -cos(d)(cos 2 (α) + 2n -3) d sin(d) κ(x, y) = (2n -4)(1 -cos(d)) + | cos(d)| -cos(d) d sin(d) .
To get κ(x, u) with u ∈ T x S 2n-1 , we just make d tend to 0. The denominator in the formulas is equivalent to d 2 , so we have to compute the second order Taylor expansion of the numerator. For κ, we get κ(x, u) = n -2 for every u not colinear to ix. For κ, we have:

cos 2 (d) sin 4 (α) + 4 cos 2 (α) = (1 -d 2 )(1 -cos 2 (α)) 2 + 4 cos 2 (α) + O(d 4 ) = (1 + cos 2 (α)) 2 -d 2 (1 -cos 2 (α)) 2 + O(d 4 ) = 1 + cos 2 (α) -d 2 (1 -cos 2 (α)) 2 2(1 + cos 2 (α)) + O(d 4 ). So we get κ(x, u) = 2n -3 2 + 3 cos 2 (α) -1 2(1 + cos 2 (α)) .
We see that κ(x, u) is not the restriction to the unit sphere of a quadratic form, but it is however the quotient of two quadratic forms. The maximal value of κ is n -1, when α = 0 or π, which is the half (remind the There are 2n -4 more dimensions not represented on this figure, which correspond to the e j 's and ie j 's (j = 3, . . . , n), but this vectors are orthogonal to those ones represented, and are left unvariant by the parallel transport along the geodesic joining x to y. While doing this parallel transport, we see the Hopf fiber, which is directed by ix, revolving around the space of codimension 2 generated by u(x, y) and the e j 's and the ie j 's, so the hyperplane Im(A(x)) in which the diffusion moves revolves on the same way.

The couplings corresponding to the extremal covariances are the isometries beween Im(A(x)) and Im(A(y)). Such a coupling will be optimal if he maximizes the mean of the scalar product of the orthogonal projections of v and w on u(x, y) ⊥ and u(y, x) ⊥ , these both spaces being identified via the parallel transport, with v and w being standard gaussian vectors in A(x) and A(y).

In the case α = π 2 , one optimal coupling consists in makig a rotation of angle d or π -d (depending on whether or not d < π

2 ), but the fact that the functionnal to minimize does not depend on the u component allows to take the coupling we can obtain from this last one by a reflection with respect to u ⊥ . If α = d = π 2 , then any coupling preserving the e j 's and the ie j 's is optimal.

For α = 0, π, the supplementary constraint that the covariance term of the distance cancels imposes that the coupling sends z to z . It remains to make the coupling between the remaining u ⊥ components of v and w. The optimal coupling in the sense of κ leaves the e j 's and ie j 's invariant, sends z to z , as previously said, and sends iz to iz if d ≤ π 2 (in which case it is a rotation), or to -iz if d ≥ π 2 (in which case it is a reflection). When α → 0 or π, with a fixed d = 0, π 2 , π, the optimal coupling in the sense of κ tends to a rotation of angle d or a reflection (in fact, that is not totally true, because the 2-plane of the rotation, or the hyperplane of the reflection depends on the direction in which y approaches the Hopf fiber except in the case when n = 2 and d < π 2 ), which is not the optimal coupling for α = 0 or π, this latter one being the parallel transport. This explains the discontinuity of κ.

The Bakry-Émery curvature of our diffusion

In this section, we prove that our diffusion does not satisfy any curvaturedimension inequality CD(ρ, N ) in the sense of Bakry-Émery (these inequalities were introduced in [1] and [2]).

In the non degenerated case, one usually take a metric adapted to the generator of the diffusion. Here, it would imply that the metric takes the value +∞ in the directions which are not orthogonal to the Hopf fibers. The distance we would obtain will be similar to the Carnot-Carathéodory distance on the Heisenberg group. So we prefer to keep the usual metric on S 2n-1 to do our computations.

We first recall the definition of the Bakry-Émery curvature of a diffusion process on a manifold. Let L = 1 2 A ij ∇ i ∇ j + F i ∇ i be the generator of a diffusion on a Riemannian manifold. The carré du champ operator is defined as

Γ(f, g) = 1 2 (L(f g) -f L(g) -gL(f )) = 1 2 A ij ∇ i f ∇ j g.
The iterated carré du champ opérator is defined as Γ 2 (f, g) = 1 2 (L(Γ(f, g)) -Γ(f, L(g)) -Γ(g, L(f ))).

The general formula of Γ 2 is a bit long and usually people do not write it, but here we need it, so it is given below

Γ 2 (f, f ) = 1 8 A kl ∇ k ∇ l A ij + 1 4 A kl A im R kml j + 1 4 F k ∇ k A ij - 1 2 A ik ∇ k F j ∇ i f ∇ j f + 1 2 A jl ∇ l A ik - 1 4 A il ∇ l A jk ∇ i f ∇ j ∇ k f + 1 4 A ik A jl ∇ i ∇ j f ∇ k ∇ l f.

2 E

 2 k , which has finite dimension. Proof : According to Lemma 9, any non-0 vector of E d,k belongs to E d , and then is an eigenvector of ∆ S for the eigenvalue -d(d + 2n -2) and also belongs to F k , and then is an eigenvector of ∆ H for the eigenvalue -k 2 . Thus such a vector is an eigenvector of L for the eigenvalue λ d,k . As we remarked just below the definition 8, if n ≥ 2, such a vector does exist as soon as d and k have the same parity and |k| ≤ d. The eigenvalue λ d,k is an integer because d(d + 2n -2) -k 2 is an even integer if k and d are both odd integers or both even integers. Doing as in the proof of Lemma 9, we can write each polynomial of degree less than d as a sum of polynomials of degree less than d whose restriction to S 2n-1 belongs to F k , so Π d = d ≤d,|k|≤d Π d ,k . Furthermore, we trivially have, by induction Π d,k = 0≤p≤ d-|k| d-2p,k , so the linear span of the E d,k is nothing but the set of the restrictions of polynomials to S 2n-1 , which is dense in L 2 , so ∆ H and ∆ S commute, and any eigenvalue of L is a λ d,k . It remains to check that the λ d,k are non-positive, and that for a given λ, the equation λ = λ d,k has a finite set of solutions among the set of pairs of integers (d, k) satisfying k + d ∈ 2Z and |k| ≤ d. For a fixed d, the maximal possible λ d,k is obtained when k = d or -d, so λ d,k ≤ λ d,d = -(n -1)d ≤ 0, so the λ d,k are non-positive. If we have λ = λ d,k , we get λ ≤ -(n -1)d, so d ≤ λ n-1 , thus we have a finite number of possible d, and for each one a finite number of possible k. Corollary 11 The spectral gap of L is n -1.

  where the last equality occurs because when k ≥ 1, a + k and b + k are positive, so each term of the sum is 0. If d 1 < d 2 , one gets 0 because Q a,b,d 1 is a polynomial of degree at most d 1 , so differentiating it d 2 times, one gets 0. So Q a,b,d 1 and Q a,b,d 2 are orthogonals for ., . a,b , as announced. Using equation 3 with d 1 = d 2 = d, one has Q a,b,d 2 a,b = (-1) d d!lc(Q a,b,d )Γ(a + d + 1)Γ(b + d + 1) Γ(a + b + 2d + 2) , where lc(Q a,b,d ) is the leading coefficient of the polynomial Q a,b,d .

n k=j+1 |z j | 2 n k=j |z j | 2 .

 22 Then we have Theorem 16 The functions f α 1 ,...,αn,d 1 ,...,d n-1 , with n j=1 α j = k and n j=1 |α j |+ 2 n-1 j=1 d j = d, form an orthonormal basis of E d,k .

  For the first fact, we use induction on d. Assume the theorem is proved for d < d. Then the set of the f α 1 ,...,α n ,d 1 ,...,d n-1 such that n j=1 |α j | + 2 n-1 j=1 d j < d is an orthonormal basis of Π d-1 , and each of this functions is orthogonal to f α 1 ,...,αn,d 1 ,...,d n-1 . So f α 1 ,...,αn,d 1 ,...,d n-1 ∈ E d , and thus from Lemma 9, we get f α 1 ,...,αn,d 1 ,...,d n-1 ∈ E d,k .

+ w 2 )

 2 cos(d)-2g 0 (v,w) sin(d) -cos(d)(g 0 (x,w)+g 0 (y,v)) 2 sin 3 (d) +O(ε 3 ). As we have x = cos(d)y + sin(d)u(y, x) and y = cos(d)x + sin(d)u(x, y), the linear term becomes g 0 (x, w) + g 0 (y, v) sin(d) = cos(d)g 0 (y, w) + sin(d)g 0 (u(y, x), w) + cos(d)g 0 (x, v) + sin(d)g 0 (u(x, y), v) sin(d) = g 0 (u(y, x), w) + g 0 (u(x, y), v) since v ∈ T x S 2n-1 and w ∈ T y S 2n-1 . That is well what was expected in the Taylor expansion of the distance.

r 1 = 2 and r 2 = cos 2

 1222 cos 2 (d) sin 4 (α) + 2 cos 2 (α) -| cos(d)| sin 2 (α) cos 2 (d) sin 4 (α) + 4 cos 2 (α) (d) sin 4 (α) + 2 cos 2 (α) + | cos(d)| sin 2 (α) cos 2 (d) sin 4 (α) + 4 cos 2 (α) 2.

  d,k -#M d-2,k , which is precisely the dimension of E d,k . So the theorem is proved. Remark 17 If one wants a real valued eigenvectors basis, we can replace the pairs of conjugated complex eigenvectors (f α 1 ,...,αn,d 1 ,...,d n-1 , f -α 1 ,...,-αn,d 1 ,...,d n-1 ) by the pairs ( √ 2Re(f α 1 ,...,αn,d 1 ,...,d n-1 ), √ 2Im(f α 1 ,...,αn,d 1 ,...,d n-1

in L) of the Ricci curvature of S 2n-1 . The minimal value of κ is n -2, when α = π 2 , so exactly the value of κ. Remark that κ is not a continuous function of u, because it jumps from the value n -1 when u is colinear to ix to the value n -2 everywhere else. This is due to the fact that our optimal coupling defining κ is not continuous.

unicity and regularity of the optimal couplings

We first recall what happends in the non degenerate case.

We saw in [4] that the unicity of the covariance of the optimal coupling is equivalent to rk(A(x)q (12) A(y)q (12)T ) = min(rk(A(x)), rk(A(y))), with rk(M ) the rank of the matrix M . In this case, the optimal coupling between the tangent spaces is unique and linear.

In the case of non degenerated diffusions, the optimal coupling defining κ was not unique. Outside the cut-locus of the manifold, the set of optimal covariances was a segment [C -, C + ]. An optimal coupling for κ realized a kind of parallel transport between the A(x)-orthogonal projection of v on the A(x)-orthogonal of u on one side, and the A(y)-orthogonal projection of w on the A(y)-orthogonal of u on the other side. The coupling between the remaining u components of v and w could be freely chosen.

In the case of a non degenerated diffusion which satisfied (H), the optimal coupling under the constraint that the variance terms vanishes was unique and linear (for x and y close enough again). It sent A(x)g(x)u on A(y)g(y)u and realized a kind of parallel transport between the g orthogonals of u in T x M and T y M (not the true one, because the gaussian measures are not exactly the same ones). This linear coupling varied smoothly in a neighborhood of each point such that rk(q (12) ) = dim(M) -1.

In our degenerated diffusion, we have r 1 r 2 = cos 4 (α), so rk(A(x)q (12) A(y)q (12)T ) = 2n-2 if α = π 2 , in which case the optimal coupling in the sense of κ is unique and linear. If α = π 2 , then r 1 = 0 and the optimal coupling is not unique, and if furthermore d = π 2 , then r 2 = 0 too and the set of optimal couplings is even greater.

With the constraint to cancel the variance term, the optimal coupling in the sense of κ is unique if and only if rk( ÃQ 12 ÃQ T 12 ) = 2n -3, which is the case if and only if d = π 2 . The following figure illustrates how do the Hopf fibers and thus the hyperplanes of the diffusion evolve along a geodesic:

The diffusion is said to satisfy a curvature-dimension inequality CD(ρ, N ), in the sense of Bakry-Émery, if for any regular enough function f , we have:

So the optimal ρ for a given N is the infimum of

Thanks to the expressions of Γ and Γ 2 , we see that

, where G 1 and G 2 are quadratic functions, and G 1 is nonnegative. At a given point x, for any u * ∈ T *

x M and any symmetric m ∈ (T * x M) ⊗2 , there exists a regular function f such that ∇f = u * and ∇ 2 f = m.

It seems then natural to try to minimize (Γ 2 (f, f ) -1 N (Lf ) 2 )(x) with a fixed ∇f . So for any u * ∈ T *

x M such that G 1 (u * ) = 0, we set

We have

and G

(2)

2 three quadratic forms, and we can assume that G

2 is not nonnegative, then the infimum is -∞, by taking m = λm 0 , with m 0 such that G (2) 2 (m 0 , m 0 ) < 0 and making λ tend to +∞. In the classical case where A ij = g ij , it was the case if and only if [START_REF] Bakry | Hypercontractivité de semi-groupes de diffusion[END_REF] N > 1 n , with n the dimension of the manifold.

If G

2 is nonnegative, then the infimum is

2 is invertible, the minimum can be written as

If no such m u * exist, the infimum is -∞ because the non existence of m u * is equivalent to the existence of a m u * such that G [START_REF] Bakry | Diffusions hypercontractives[END_REF] (m u * , m u * ) = 0 and G (12) 2 (u * , m u * ) = 0, and then we can take m = λm u * , and make λ tend to +∞ or -∞. Now let us see what we get for our degenerated diffusion. We have

If M and M are the matrices of m and m in the basis (ix, e 2 , ie 2 , . . . , e n , ie n ) of T x S 2n-1 , with the e j chosen such that (x, e 2 , . . . , e n ) is a C-orthonormal basis, then we have

with M and M obtained from M and M by removing the first line and the first column. So G

2 is nonnegative if and only if

2 , we have to compute ∇A, and thus ∇(ix), since we have

As previously, we take u ∈ T x M an unit norm vector and y = exp x (du). Then ∇ u (ix) = d dd τ yx (iy)| d=0 , where τ yx is the parallel transport along the geodesic joining y to x. This parallel transport corresponds to a rotation in the 2-plane generated by x and u. We have:

wich is the orthogonal projection of iu on T x S 2n-1 . Let us denote by J(x) the endomorphism of T x S 2n-1 which sends u to iu + g 0 (u, ix)x. Then g(x)J(x) is antisymmetric, J(x)J(x) = -I + (ix) ⊗ (ix) * and J(x)(ix) = 0. We have ∇ i (ix) j = J j i .

Thus we get

where we used the antisymmetry of J(x)g -1 (x) and the symmetry of m jk to cancel the second term and to remark that the three other ones were equal (up to the prefactor 1 2 or 1 4 ). Here we are in the case where m u * does not exist since for any u * not colinear to (ix) * (in which case G 1 (u * , u * ) = 0), we can take m u * = 1 2 ((J(x) T u * ) ⊗ (ix) * + (ix) * ⊗ (J(x) T u * ).

Thus ρ N (x, u * ) = -∞, and we have no Bakery-Émery curvature-dimension inequality. So in our example, it is useful to consider a metric not adapted to the diffusion and to use the curvatures κ and κ (to get Poincaré inequalities when n ≥ 3, for example).


Jiajie Wang

Patha S Roop
email: p.roop@auckland.ac.nz

Alain Girault
email: alain.girault@inria.fr

Energy and timing aware synchronous programming *

The synchronous paradigm is widely used for the design of safety critical systems. Such systems, especially in the medical devices domain, must meet strict timing requirements while also ensuring long battery life. As a consequence, they are subject to very strict constraint both regarding their WCRT (Worst-Case Reaction Time) and their WCEC (Worst-Case Energy Consumption, the equivalent constraint for the energy consumption). Many techniques exist to compute an upper bound on the WCRT, but few techniques exist that address both the WCRT and the WCEC. We propose here a static analysis framework where conventional WCRT analysis interacts with a DVFS (Dynamic Voltage Frequency Scaling) algorithm to minimize also the WCEC of the given synchronous program. Our algorithm is able to compute the Pareto front of non-dominated solutions in the (WCRT, WCEC) space. Experimental results reveal that the proposed approach is scalable in terms of analysis time while providing more non-dominated solutions compared to two existing approaches. To the best of our knowledge, the proposed approach is the first to produce energy and timing aware synchronous programs.

INTRODUCTION

The synchronous paradigm [START_REF] Benveniste | The synchronous languages 12 years later[END_REF] has proved to be very useful in designing safety-critical systems [START_REF] Dormoy | SCADE 6: a model based solution for safety critical software development[END_REF][START_REF] Huhn | Towards certifiable software for medical devices: The pacemaker case study revisited[END_REF]. The advantages are that compilers are grounded on sound semantic foundations based on causality and constructiveness, and hence generate correct-by-construction code that guarantees determinism and reactivity. Determinism ensures that, from any state, given any input vector, at most one reaction is enabled. Reactivity (also called "input receptiveness") ensures

Related Work

A processor capable of DVFS may operate at several discrete frequencies. To ensure correct processor operation, each frequency is accompanied by a voltage setting. Thus a processor has several operation modes, each one being a pair (voltage, frequency), denoted (V, f ). A point during program execution, where the processor changes its operation mode, is referred to as a DVFS control point. DVFS can be controlled by the software to optimize the performance and energy consumption, and find good trade-offs, by changing the operation mode dynamically. For a real-time system, the energy consumption can be reduced by consuming the slack (i.e., the time between task/program completion and the deadline) in exchange for a lower processor frequency.

DVFS algorithms need to make two decisions: the choice of the control point locations in the program and the frequency value associated with each control point. DVFS algorithms may be broadly classified as online, offline, or hybrid. Online [START_REF] Aydin | System-Level Energy Management for Periodic Real-Time Tasks[END_REF][START_REF] Aboughazaleh | Collaborative Operating System and Compiler Power Management for Real-time Applications[END_REF] algorithms decide both the control points and the frequency values during program execution over a Real-Time Operating System (RTOS). They are closely associated with the underlying scheduling policy and hence are scheduler dependent. Offline algorithms [START_REF] Xie | Compile-time Dynamic Voltage Scaling Settings: Opportunities and Limits[END_REF][START_REF] Huang | Leakage-aware intraprogram voltage scaling for embedded processors[END_REF], on the other hand, make both decision at compile time. Hybrid approaches [START_REF] Seo | Profile-based optimal intra-task voltage scheduling for hard real-time applications[END_REF][START_REF] Huang | Efficient and Scalable Compiler-directed Energy Optimization for Realtime Applications[END_REF], combine offline and online decisions: DVFS control points are inserted directly into the program statically, but the frequency values of the control points are decided on-line. Both the online and hybrid approaches can be classified as best effort techniques: they rely on extracting the run-time information, such as the actual execution times, in order to compute the next processor frequency. Thus, their effectiveness can only be evaluated using simulation, which may yield variable results depending on the simulation and benchmark settings [START_REF] Kim | Performance comparison of dynamic voltage scaling algorithms for hard real-time systems[END_REF]. Hence, these techniques are unsuitable for safety-critical medical devices, which require strict guarantees.

To ensure that bounds of the system can be computed, we have to consider offline DVFS schemes. These algorithms make compile time decisions regarding both the control point locations and the associated frequency values. The problem may be considered as a bicriteria optimization problem over two antagonistic goals: execution time and energy consumption. For instance, [START_REF] Xie | Compile-time Dynamic Voltage Scaling Settings: Opportunities and Limits[END_REF][START_REF] Huang | Leakage-aware intraprogram voltage scaling for embedded processors[END_REF] solve this problem by first abstracting the computation of timing. Instead of considering all possible execution paths, they only consider a single execution path derived using profiling. Then, they use Integer Linear Programming (ILP) to optimize the energy consumption, where the execution time of the program is expressed as an ILP constraint. The drawback of such linearized approaches is that they may produce underapproximations, which are not desirable in safety-critical systems.

Our approach

We propose, for the first time, a compile time DVFS scheme for synchronous programs. While increasing the frequency reduces the WCRT, such an increase will increase the Worst-Cases Energy Consumption (WCEC). Hence, these two criteria are antagonistic.

To compare tradeoffs in the (WCRT,WCEC) space, we rely on the notion of dominance and Pareto optima:

• The point (x, y) weakly dominates the point (x , y ) iff (x < x ∧ y = y ) ∨ (x = x ∧ y < y ).

• The point (x, y) strongly dominates the point (x , y ) iff (x < x ∧ y < y ).

• A point is a weak Pareto optimum iff there does not exist another point that strongly dominates it (but there might exist other points that weakly dominate it).

• A point is a strong Pareto optimum iff there does not exist another point that (weakly or strongly) dominates it.

• The Pareto front is the set of weak and strong Pareto optima.

To avoid potential under-approximations like existing techniques, we do not linearize the problem. Instead, we perform bicriteria minimization using an adaptation of the wellknown ε-constraint method [START_REF] Haimes | On a bicriterion formulation of the problems of integrated system identification and system optimization[END_REF]. This method is used while performing optimization for two antagonistic criteria such as the WCRT and the WCEC. The ε-constraint method builds the Pareto front incrementally by transforming one of the two criteria into a constraint and minimizing the other criterion under this constraint to obtain one tradeoff. Then, a new value is chosen for the constraint, therefore obtaining a new tradeoff, and so on. At the end, all the dominated tradeoffs are removed to obtain the Pareto front. In our case, we transform WCRT into a constraint, i.e. WCRT ≤ δ, where δ is a given deadline value. We then minimize the WCEC under this constraint. Each value of δ will potentially produce a different solution in the (WCRT, WCEC) plane. More points are obtained by varying the deadline value by a fixed amount. The developed algorithm is detailed in section 4.

Additionally, we use the four DVFS frequency points of the selected processor to create four points on the Pareto front. This approach is known as the fixed frequency approach and it has been established that using fixed frequencies will produce dominant points on the Pareto front [START_REF] Rizvandi | Some observations on optimal frequency selection in DVFS-based energy consumption minimization[END_REF]. Since these frequencies are far apart, the resulting fixed frequency points are also far apart in the (WCRT,WCEC) plane. The goal of our proposed method is precisely to find more tradeoffs between these fixed-frequency points. Our algorithm, while being sub-optimal, performs no backtracking and hence is very efficient and scalable. Due to suboptimality, the points obtained may not be guaranteed to be dominant. However, through benchmarking, we show that the algorithm produces dominant points in many benchmarks.

The main contributions of this paper is summarized below:

• We propose the first framework that allows DVFS to be used with the synchronous paradigm for the design of safety-critical systems [START_REF] Alemzadeh | Analysis of safety-critical computer failures in medical devices[END_REF].

• We developed an analysis technique that tackles the bicriteria optimization problem (time, energy) for synchronous programs, without using linearization.

• We have benchmarked our method with real-life synchronous programs. The Pareto optima that we produce dominate the tradeoffs produced by the linearized approach. Moreover, we produce a lot more Pareto optima than the four points produced by the fixed frequency approach, therefore giving the end user more tradeoffs to consider.

The remainder of the paper is organized as follows: We introduce synchronous languages and the target intermediate format in Section 2. Then we formalize the proposed compile time DVFS scheme in Section 3. In Section 4, we present the overview of the proposed framework. Subsequently, the interactive steps of timing analysis and the DVFS algorithm are detailed in Section 5 and Section 6 respectively. We present the evaluation of the proposed approach relative to two existing approaches in Section 7. Finally, Section 8 concludes with discussions on limitations and scope for future work.

BACKGROUND

Synchronous languages

The foundation of synchronous paradigm is the synchrony hypothesis, which assumes that the execution of the program is faster than the environment. Synchronous programs execute in discrete instants, call ticks. During each tick, the environment is sampled at the beginning, then the computations take place, and the outputs are emitted at the end. For the synchrony hypothesis to hold, the execution time of any tick must be shorter than the minimum arrival interval of events, and the analysis to compute the longest execution time of a tick is known as WCRT analysis [START_REF] Kuo | Efficient WCRT analysis of synchronous programs using reachability[END_REF][START_REF] Wang | ILPc: A novel approach for scalable timing analysis of synchronous programs[END_REF][START_REF] Raymond | Timing analysis enhancement for synchronous program[END_REF]. The technique we propose in this paper is based on the synchronous language PRET-C and its intermediate format Timed Concurrent Control Flow Graph (TCCFG) [START_REF] Andalam | A Predictable Framework for Safety-Critical Embedded Systems[END_REF]. TCCFG captures both the high-level control flow of PRET-C programs and the timing information of the underlying platform. As an example, a TCCFG is presented in Fig. 1, which demonstrates a subset of its features. A TCCFG is a directed graph {B, E} where B = {B0, B1, B2. . . } is the set of nodes and E = {E0, E1, E2. . . } is the set of edges. The control flow of a PRET-C program is captured using various notations in TCCFG. The start/end nodes (e.g., B0 and B10) mark the beginning and the end of the program. Blocks of sequential computation instructions are encapsulated inside computation nodes (e.g., B4, B6 and B8), and conditional branches ("if" and "else") are captured by condition nodes (e.g., B1). Concurrent threads are respectively forked and joined using fork and join nodes (e.g., B2 and B9). When the parent thread reaches a fork node, it spawns all its child threads and suspends itself, and it resumes when the last child thread terminates (i.e., reaches the join node). The "End of Tick" nodes (EOT) mark the tick/state boundaries of the threads. During each tick, all active threads execute until they have reached their respective EOT nodes, which is therefore a synchronization barrier. When all the active threads have reached their barrier, the program advances to the next tick, and the execution of threads resumes from the EOT nodes where they stopped at the end of the previous tick. It follows that concurrent threads execute in a lockstep manner. PRET-C threads are interleaved only at the tick granularity, i.e, thread switching can only occur after an EOT or a join node.

Timed Concurrent Control Flow Graph

A fragment of the TCCFG is outlined in Fig. 1 with computation nodes B4, B6, and B8, which will serve as a motivat-ing example to illustrate the proposed framework. Regarding the timing information, each node is annotated with the exact number of processor cycles required to execute it on the target processor. This number is extracted by analyzing the binary code [START_REF] Andalam | A Predictable Framework for Safety-Critical Embedded Systems[END_REF]. To have a clear presentation, only the execution costs of B4, B6, and B8 are shown in the outlined fragment, which are assumed to be 90, 30, and 110 processor cycles respectively.

Worst-Case Energy Consumption (WCEC)

We want to apply the principle of computing the WCRT as the maximal execution time of all ticks to the energy, with the Worst-Cases Energy Consumption (WCEC). The motivation is to derive, for instance, a lower bound on the life time of a battery powered embedded system. Each PRET-C program will therefore be characterized by its two WCRT and WCEC values. Then, in order to compute several tradeoffs between energy and time, we shall use DVFS, as explained in the next section.

DVFS SCHEMES FOR SYNCHRONOUS PROGRAMS

The core of the proposed framework is a compile time DVFS framework to produce tradeoffs in the (WCRT,WCEC) plane. Our DVFS framework is tightly coupled with the synchronous semantics: the program has a set of control points C, which are the boundaries at which the context switching between threads occurs, i.e., the start, EOT and join nodes. These are ideal places for altering the frequency as they are well marked in synchronous specifications. For example, for the TCCFG in Fig. 1, the set of DVFS control points is C = {B0, B3, B5, B7, B9}. In addition, we are given a set of pre-defined processor frequencies F. Let the cardinality of C be n and that of F be m. Then, there are m n possible alternatives to assign frequency values to the control points. Each such assignment is known as a DVFS scheme. The goal of the paper is stated below: Goal: Given a set C of DVFS control points in a program and a set F of available frequencies on the processor, a DVFS scheme is a function λ : C → F. Given the intermediate TCCFG representation of the program and a timing deadline T , the objective of this paper is to compute a suitable λ from the search space, such that the WCRT of the resulting program is less than T and its WCEC is as small as possible. Then, by varying the deadline T , it is possible to obtain a set of tradeoffs, from which all the dominated point (in the Pareto sense, see below) can be removed to obtain the Pareto front. This is the key challenge of the proposed framework. In Sec 4, we present our static analysis technique which solves this problem by combining an ILP-based WCRT analysis with a greedy heuristic.

A motivating example

In this section, we illustrate how we compute the execution time and energy consumption with a DVFS scheme using the running example outlined in Fig 1 . The program starts from the start node B0, and executes the condition node B1, which has two branches E2 and E9. The branch E2 leads to the fork node B2, which subsequently spawns two threads and eventually leads to the execution of the computation nodes B4 and B6. The alternative branch E9 leads to the execution of the computation node B8. In other words, the program either executes B8 or executes B4 and B6 concurrently. The computation nodes B4, B6 and B8 may execute at different frequencies, depending on the values of the DVFS control points at the EOT nodes B3, B5 and B7 respectively. Moreover, since B0 is also a DVFS control point, nodes B1 and B2 should execute at the frequency set in B0. Yet, since there is no computation node between B0 and B3/B5/B7, for the sake of simplicity we do not take this frequency into account in the computations presented in the paper, and similarly for the join node B9. But of course, our tool does take it into account. The execution time of a node varies with the processor frequency. In this paper, we let the execution time of a node executing at the fastest frequency to be equal to its execution cost in processor cycles, and scale the execution time linearly for the slower frequencies. For example, if the running example is executing on a processor that has four frequencies -0.25Mhz, 0.5Mhz, 0.75Mhz and 1Mhz -then the execution time of the computation nodes B4, B6 and B8 at different frequencies are as in Table 1. The energy consumed by the execution also varies according to the frequency (Table 2). The energy consumption of a node at frequency f is computed by multiplying its execution cost in processor cycles with the energy per cycle at frequency f . The value of the energy per cycle is estimated as follows [START_REF] Burd | Energy efficient CMOS microprocessor design[END_REF]:

Energy per cycle at f P dynamic × cycle length (C × V 2 × f ) × (1/f )
where C is the switching capacitance. This is a typical estimation for COMS circuitry, based on the classical model of dynamic power consumption. Additionally, since the voltage can be assumed to be proportional to the frequency (i.e., V ∝ f [START_REF] Burd | Energy efficient CMOS microprocessor design[END_REF]) we may safely assume that the energy per cycle is ∝ f 2 .

Comparing DVFS schemes

Since the number of DVFS schemes grows exponentially, it is important that we determine the relative superiority of one scheme over another. We compare these schemes using their WCRT and WCEC and the notion of Pareto dominance (see Section 1.2). For the running example, the computations of WCRT and WCEC can be abstracted as max operations. Let t(n) denotes execution time of the node n and e(n) denotes its energy consumption. Then, for this running example, the WCRT and WCEC are calculated as: WCRT = max(t(B8), (t(B4) + t(B6)))

(1) WCEC = max(e(B8), (e(B4) + e(B6))) WCRT and WCEC are antagonistic in nature: reducing the WCRT requires to increase the frequencies, which in turn increases the WCEC. Similarly, reducing the WCEC requires to decrease the frequencies, which in turn increases the WCRT. Table 3 shows four of the DVFS schemes for the running example, with their corresponding WCRT and WCEC. The corresponding points are plotted in Fig 2 .  Scheme 2 is clearly worse than scheme 1 as it has the same WCRT but a larger WCEC: in other words, scheme 1 weakly dominates scheme 2. The portion of the plane highlighted in gray contains all the tradeoffs dominated by scheme 1. Scheme 3 and scheme 4 are distinct in terms of frequency values, but they have identical WCRT and WCEC: this is because we are interested in the worst case values (but they could have distinct average case values for instance). Finally, scheme 3 and scheme 1 are incomparable: none dominate the other one. Still, if we are given a strict deadline on the WCRT (as plotted in Fig. 2), then scheme 1 is the best tradeoff.

ALGORITHM FOR WCRT AND DVFS SCHEME COMPUTATION

Given a synchronous program and a timing deadline T , we have to (i) determine the frequency values for the proposed DVFS scheme as defined in Section 3, and to (ii) ensure that WCRT ≤ T . To solve this problem, we have developed a framework illustrated in Fig. cess interacts with the other one until a positive or negative answer is returned. This approach allows us to avoid under-approximation caused by abstracting the program execution.

The developed approach, though iterative, has a wider scope than the recent iterative algorithms [START_REF] Raymond | Timing analysis enhancement for synchronous program[END_REF][START_REF] Wang | ILPc: A novel approach for scalable timing analysis of synchronous programs[END_REF]. The existing techniques are timing analyses that use Implicit Path Enumeration Technique (IPET) to reduce the search space, therefore improving the scalability for WCRT computation. In our framework, the timing analysis is only one of the internal processes. We extend the classical timing analysis for WCRT with a DVFS analysis to compute DVFS schemes that have a WCRT less than a given deadline and a WCEC that is as small as possible.

Our framework starts with a TCCFG and an initial DVFS scheme where all DVFS control points are set to the slowest frequency (0.25MHz). The analysis first computes the WCRT and its corresponding execution path 1 . Then the WCRT is compared to the timing deadline 2 . If the WCRT is longer than the deadline, then the DVFS algorithm will take the corresponding execution path of the WCRT, and increase the frequencies of the DVFS control points along that path to reduce its execution time 3 . Then the analysis starts over again with the WCRT analysis 1 , to search for the other execution paths that exceed the deadline. As soon as the computed WCRT is shorter than the deadline 2 , the analysis finishes: the DVFS scheme along with the WCRT and its corresponding execution path are reported. Finally the WCEC of the reported scheme is computed 4 . If the DVFS algorithm 3 cannot make an execution path shorter than the deadline (i.e., all DVFS control points are already at the highest frequency), then the analysis terminates with an error, stating that the deadline is not achievable.

In the following sections, we will present more details on the DVFS scheme adjustments 3 , which uses a heuristic algorithm to further reduce the analysis time, as well as the analysis technique developed for computing the WCRT 1 and the WCEC 4 .

WCRT AND WCEC ANALYSIS

Computing the WCRT and WCEC follow are similar processes as they determine the program path that maximizes a system property, namely the execution time and energy consumption respectively. We propose here a novel ILP-based technique to compute the WCRT and WCEC from a given TCCFG, using execution time and energy cost for respective analysis.

In the following subsections, we will first present the existing ILP formulation for modeling the control flow of a TCCFG [START_REF] Wang | ILPc: A novel approach for scalable timing analysis of synchronous programs[END_REF], then our extension for taking into account DVFS scheme. The TCCFG in Fig. 1 will be used as the running example. For simplicity, the outflow edge of an EOT node or a start node is termed an EOT edge.

Modeling the high level control flow

As in IPET [START_REF] Raymond | A general approach for expressing infeasibility in implicit path enumeration technique[END_REF], the objective function of the formulation is to maximize the sum of all the edges multiplied by their associated costs. For example, the objective function for the example TCCFG is:

max N i=1 Ei × ci ( 3
)
where N = 12 in this example, and ci denotes either the associate cost of Ei, which is the execution time in the WCRT analysis, or the energy cost in WCEC analysis.

As for the ILP constraints, two assumptions are made based on the properties of synchronous programs and the nature of ILP. First, the value of each edge is bounded to be either '1' (active) or '0' (inactive). This assumption is valid for any synchronous program as they do not have recursion and instantaneous loops. Second, an edge is set to '1' whenever possible, since the execution/energy costs of all nodes have positive contributions toward the objective function (i.e., all costs are positive values).

Computation and condition nodes

The control flow at computation and condition nodes are modeled as inflow edges equals to the outflow edges, as the number of times the control flow enters a computation or condition node must equal to the number of time it leaves the node. For example, the following constraints are created for the running example:

B1: E1 = E2 + E9 B4: E4 = E5 B6: E7 = E8 B8: E10 = E11

EOT nodes

The ILP constraints for EOT nodes are used for emulating the execution of one tick and enforce the parent-child relationships. This is achieved by eliminating the EOT edge combinations, which are unfeasible. For example, the following ILP constraints are generated for the motivating ex-ample:

T0-T1: E1 + E4 ≤ 1 T0-T2: E1 + E7 ≤ 1 T0-T3: E1 + E10 ≤ 1 T1-T3: E4 + E10 ≤ 1 T2-T3: E7 + E10 ≤ 1
Each ILP constraint consists of a group of EOT edges which are mutually exclusive, i.e., at most one of them can be '1'. For example, the constraint "T0-T1" prevents the threads T0 and T1 to execute concurrently, since T0 is the parent thread of T1. Similarly, the other constraints eliminate other unfeasible combinations. Only T1 and T2 are allowed to execute concurrently.

Fork nodes

Fork nodes spawn child threads when the control flow reaches them. Therefore, all the outflow edges of a fork are active if any of its inflow edge is active. The ILP constraints for the fork node B2 are as follows:

B2: E2 = E3, E2 = E6

Join nodes

The execution of a join node is denoted by its outflow edge. When the last child thread terminates, the outflow edges of a join node is active; otherwise it should remain inactive. Given that the ILP solver will set the outflow edge to '1' if no ILP constraint is created for the join node (see assumptions at the beginning of Section 5.1), our formulation only needs to model the cases where the outflow edge must be '0'. For example, the following ILP constraints are generated for the join node B9:

LastChild: E5 + E8 ≥ E12 T1: (1 -E4 -E2) + (E5) ≥ E12 T2: (1 -E7 -E2) + (E8) ≥ E12
The outflow edge E12 is on the right hand side of all the ILP constraints. It is '0' if any of the left hand side of the constraints is '0'. Otherwise it is set to '1' by the ILP solver. The LastChild constraint captures the moment when the last child thread reaches the join node B12, during which at least one of E5 and E8 must be active.

The constraints for T1 and T2 capture the states of the two child threads respectively. A child thread forces the outflow edge of the join node to a value '0' if (i) it has an active edge, and (ii) it does not reach the join node. This is captured as two parts in the left hand side of the constraints (separated by parentheses). For example, if any edge in T1 is active, then either its EOT edge E4 is active or the inflow edge of the fork node E2 is active, since all the edges in T1 are initiated from them. So the first part "1-E4-E2" is '0' if this is true (E4 and E2 are mutually exclusive from each other because of the constraints for EOT nodes). The second part "E5" is '0' if the control does not reach the join node. If any of the two parts is '1', then the outflow edge E12 is set to '1' provided that T2 also meets the same condition.

Modeling the DVFS scheme

Processor cycle has been the primary unit in conventional timing analysis. The actual physical time is calculated by multiplying the number of clock cycles with the cycle length (which is 1/f ). However, when DVFS is employed, the processor frequency changes during the execution, and this effect must be taken into consideration in the timing analysis. This changes of frequency complicates the ILP formulation as it introduces multiple costs for each node.

Conventional IPET [START_REF] Raymond | A general approach for expressing infeasibility in implicit path enumeration technique[END_REF] cannot model nodes with variable costs. To cope with this limitation, we consider a TCCFG with DVFS as multiple virtual TCCFGs, where each TCCFG corresponds to one specific frequency in F. An illustration of this for our example TCCFG is shown in Fig. 4. The frequency switching is emulated as the control flow jumping between these virtual TCCFGs.

The virtual TCCFGs are identical to the original, except that the cost of each node is based on the associated frequency. The calculations of execution and energy cost are shown in Section 3.1. The switching penalty is modeled as a constant cost and is added to the start node, and all the EOT and join nodes (i.e., all the nodes in C). Now, the virtual TCCFGs introduce multiple avatars for each edge in the original TCCFG. For example, with processor frequencies 0.25MHz, 0.5MHz, 0.75MHz and 1MHz, an edge Ei in the original TCCFG becomes four avatars, one for each frequency. We differentiate these avatars using suffixes, so that the avatars of Ei are Ei_a, Ei_b, Ei_c and Ei_d, corresponding to 0.25MHz, 0.5MHz, 0.75MHz and 1MHz respectively. Each avatar, like the original edge, is a binary value.

At any time, only one avatar of an edge can be active in an execution path, since a node can only execute at one frequency. Thus the sum of all the avatars of an edge should be at most '1'. To enforce this, the following ILP constraint is created for each edge:

Ei a + Ei b + Ei c + Ei d ≤ 1
After creating all the avatars, we rebuild all the ILP constraints for a TCCFG with these avatars. We developed a systematic approach, which can preserve the modeling of the high level control flow while integrating the DVFS scheme. Each ILP formulation is treated differently depending whether the node has a DVFS control point or not.

Nodes without DVFS control points

For the computation nodes, condition nodes and fork nodes, which do not have frequency control points, we consider each virtual TCCFG independently, and generate ILP constraints for each using the corresponding avatars. For example, the following ILP constraints are generated for the running example: 0.25MHz TCCFG: to each virtual TCCFG. However, the processor frequency cannot be changed during the execution of these nodes. For example, if the condition node B1 executes at 0.5MHz, then its inflow edge E1_b is '1' (the other avatars are '0'). The ILP constraint will allow B1 to choose its outflow edge, but the processor frequency must remain to be 0.5MHz (i.e., the only choices are E2_b and E9_b).

B1: E1 a = E2 a + E9 a

Nodes with DVFS control points

For the EOT and join nodes that are in C, we consider the all the virtual TCCFGs as a whole, and create a combined set of ILP constraints to allow the control flow to jump between them. We replace each edge in the original ILP formulation with the sum of its avatars. For example, the following ILP constraints are generated for our running example:

EOT nodes:

T0-T1: (E1 a + E1 b + E1 c + E1 d)+ (E4 a + E4 b + E4 c + E4 d) ≤ 1 T0-T2: (E1 a + E1 b + E1 c + E1 d)+ (E7 a + E7 b + E7 c + E7 d) ≤ 1 . . . Join node B9: LastChild: (E5 a + E5 b + E5 c + E5 d)+ (E8 a + E8 b + E8 c + E8 d) ≥ (E12 a + E12 b + E12 c + E12 d) T1: [1 -(E4 a + E4 b + E4 c + E4 d)+ -(E2 a + E2 b + E2 c + E2 d)]+ (E5 a + E5 b + E5 c + E5 d) ≥ (E12 a + E12 b + E12 c + E12 d) T2: [1 -(E7 a + E7 b + E7 c + E7 d)+ -(E2 a + E2 b + E2 c + E2 d)]+ (E8 a + E8 b + E8 c + E8 d) ≥ (E12 a + E12 b + E12 c + E12 d)
These ILP constraints preserve the modeling of the control flow while allowing the frequency to be changed. For example, if T1 executes at 0.25MHz (i.e., E4_a and E5_a are '1') and T2 executes at 0.5MHz (i.e., E7_b and E8_b are '1'), then the control flow of the two threads are on different virtual TCCFGs, but the ILP constraints still allow their join node B9 to be active. Furthermore, the join node B9 can switch to any frequency (i.e., it is free to select any avatars of E12).

Choosing a frequency

Initially, the ILP solver is free to choose any frequency for the control points. As the analysis proceeds, the DVFS algorithm will select specific frequencies for the control point to reduce the execution time. This selection is modeled differently for the WCRT and WCEC analyses. The ILP solver always prioritizes the frequency that maximizes the objective function: Ei_a > Ei_b > Ei_c > Ei_d in the WCRT analysis, and in reverse order in the WCEC analysis. Therefore, we can choose a specific frequency by setting the unwanted ones to '0', which forces the solver to take the next best. For example, to choose 0.5MHz for the EOT node B3 (i.e., E4_b to be '1'), we create the following ILP constraints respectively for the WCRT and WCEC analyses:

In WCRT analysis:

In WCEC analysis:

E4 a = 0 E4 d = 0 E4 c = 0

THE DVFS ALGORITHM

We propose a greedy heuristic for the DVFS algorithm for adjusting the frequency of the DVFS control points along an execution path. The idea of the algorithm is to reduce the execution time of the path as much as possible during each frequency adjustment. To do so, the algorithm only increases the frequency values (i.e., frequency values are never decreased). It adjusts one DVFS control point with one frequency step at a time and it prioritizes the DVFS control points which yield the largest reduction in execution time.

For example, let us consider an execution path from the running example in Fig. 1, whose execution time is the sum of the execution time of B4 and B6. Let us assume both nodes are executing at 0.25MHz, which results in an execution time equal to 360+120. The algorithm will select B3 and increase its frequency by one step, to 0.5MHz (execution time: 180+120). Indeed, this adjustment yields the largest reduction in execution time compared with increasing the frequency of B5 by one step (execution time: 360+60). The execution time of the whole execution path is re-calculated with each frequency adjustment, and the process repeats until the execution time of the path is shorter than the deadline (or a failure is reported if it is not possible). This heuristic algorithm effectively reduces the search space.

RESULTS

Recall that for a multi-criteria optimization problem, there is no single optimal solution. In this section, we compute the estimated Pareto front for a set of benchmark programs using the proposed approach, and compare the results with two existing approaches: the fixed frequency approach, which uses a single processor frequency throughout the execution, and the linearized approach [START_REF] Xie | Compile-time Dynamic Voltage Scaling Settings: Opportunities and Limits[END_REF][START_REF] Huang | Leakage-aware intraprogram voltage scaling for embedded processors[END_REF]. We first discuss the results in terms of WCRT and WCEC, followed by an comparison in analysis time.

Benchmark method

The target platform uses a Microblaze processor with four frequencies: 1Mhz, 0.75Mhz, 0.5Mhz, and 0.25Mhz. The frequency switching penalty is assumed to be five processor cycles at 1Mhz. Our benchmark programs are industrial applications taken from [START_REF] Yoong | Implementing constrained cyber-physical systems with IEC 61499[END_REF]. For each benchmark, we build the Concurrent Control Flow Graph (CCFG), then we compile it into C code and then into binary to obtain the exact number of processor cycles of each node, therefore resulting in the TCCFG (that is, the Timed CCFG). The benchmarks are summarized in Table 4.

We applied the three approaches to each of the benchmark programs. For our bicriteria algorithm, the initial deadline is selected to be the WCRT of the program executing at 1MHz (i.e., the tightest deadline that the system can meet, where all DVFS points are set to the fastest frequency). The deadline increases by 20% of the starting value each time until the deadline is equal to the WCRT of the program executing at 0.25MHz (i.e., all DVFS points are set to the slowest frequency).

Computing the Pareto front

The fixed frequency approach serves as a baseline for comparison, which represents the case of executing a conventional synchronous program at a fixed frequency, without using DVFS. There is no frequency switching penalty, and the WCRT and WCEC are computed using the ILP technique presented in Section 5. The fixed frequency approach is benchmarked for each of the four available frequencies for each program. The results are plotted in the (WCRT,WCEC) plane in Fig. 5 where each blue triangle represents one such fixed frequency tradeoff. All these four points are strong Pareto points, meaning that there does not exist any point better in both criteria WCRT and WCEC. In Fig. 5, we also outline the portion of the (WCRT,WCEC) plane dominated by each fixed frequency point.

The linearized approach handles the bicriteria optimization problem by first linearizing the program execution path using profiling [START_REF] Xie | Compile-time Dynamic Voltage Scaling Settings: Opportunities and Limits[END_REF][START_REF] Huang | Leakage-aware intraprogram voltage scaling for embedded processors[END_REF]. To implement this principle, we use the critical path obtained through conventional WCRT analysis as the profiled path, and we optimize the energy over that path. We also restrict the allocation of DVFS control points to be at the start, EOT, and join nodes to be consistent with the proposed approach. The results are plotted in Fig. 5 as green circles. Most of the schemes produced by the linearized approach overlaps with each other, thus only three points are visible in each graph. All of these points are dominated by the points produced by the other approaches. More importantly, due to the abstraction, the linearized approach misses the given deadline most of the time.

Finally, for each benchmark program, we produce tradeoffs with our bicriteria approach by varying the deadline (Section 7.1). The results are plotted in Fig. 5 as red squares.

Compared with the fixed frequency approach, our approach is disadvantaged due to the frequency switching penalties. However, we produce a lot more tradeoffs. For instance, for the Channel Protocol (Fig. 5a), the Pareto front consists of 12 points: 4 are fixed frequency while 8 are obtained with our approach (recall that the Pareto front consists of all the points that are non-dominated). This gives a lot more choices to the end user. Consider again the Channel Protocol and suppose that the user wishes to find the best DVFS scheme such that the WCRT be less than 3500 clock cycles, where "best DVFS scheme" means "the DVFS scheme with the smallest WCEC" because the WCRT aspect is already taken care off through the 3500 clock cycles constraint. As shown in Fig. 5a, the tradeoffs with the smallest WCEC found are: linearized method: WCEC = 766 fixed frequency method: WCEC = 256 our approach: WCEC = 176

Of course, the ranking between the three approaches will depend on the deadline chosen by the user, but it is obviously better to have more points in the Pareto front, which demonstrates the usefulness of our approach. Finally, in some cases our approach fails to provide nondominated solutions in key portions of the search space. This is the case of the Flasher and the Cruise Controller, where our approach failed to provide non-dominated solutions between the two fixed frequency Pareto optima corresponding respectively to 0.5MHz and 0.75MHz (see Figs. 5b and5d). A future work direction will be to investigate how to improve the DVFS schemes found by our algorithm by a local search in the search space.

Analysis time

Analysis time is an important quantitative consideration. Off-line DVFS algorithms are hard to scale due to the large search space [START_REF] Huang | Efficient and Scalable Compiler-directed Energy Optimization for Realtime Applications[END_REF]. The analysis time for the largest benchmark program, the Cruise Controller, is recorded in Table 5. All three approaches are practical to use. The fixed frequency approach does not require any analysis, and both the proposed approach and linearized approach finish the analysis in matter of seconds.

The linearized approach has a constant analysis time as its search space does not change. On the other hand, the analysis time of the proposed approach increases as the deadline get tighter. This is because the starting point of our algorithm is the DVFS scheme where all control points are set to the smallest frequency (0.25MHz), so a lot of iterations must be performed before reaching a DVFS scheme that meets a tight deadline. A future work direction will be to investigate whether it is doable to initiate our algorithm with another fixed frequency scheme, e.g., 0.75MHz for deadlines less than 2500 clock cycles in the case of the Robot Sonar. However, given the fact that the search space of the proposed approach is much larger than the linearized approach, it is nice to observe that they have comparable analysis time. We conclude that the greedy approach that performs no backtracking for computing a DVFS scheme (λ) significantly reduces the analysis time.

CONCLUSIONS AND FUTURE WORK

Synchronous languages provide efficient programming of embedded safety critical systems. Recent results have focused on computing their Worst-Cases Reaction Time (WCRT) with ILP based methods improved with IPET [START_REF] Kuo | Efficient WCRT analysis of synchronous programs using reachability[END_REF][START_REF] Wang | ILPc: A novel approach for scalable timing analysis of synchronous programs[END_REF][START_REF] Raymond | Timing analysis enhancement for synchronous program[END_REF]. This is perfectly relevant to validate whether a given system meets its timing constraint. However, modern embedded systems must also satisfy strict energy consumption constraints, especially for tiny devices such as sensor networks, biomedical implants, and ICD devices. Indeed, premature battery wear-off is a life threatening issue in such devices [START_REF] Alemzadeh | Analysis of safety-critical computer failures in medical devices[END_REF]. Regarding synchronous programs, the energy equivalent of the WCRT is the Worst-Cases Energy Consumption (WCEC), which is a guaranteed bound of the energy consumption of the system during one reaction.

We have proposed an original bicriteria (WCRT,WCEC) static analysis optimization method to tackle this. It combines a classical WCRT analysis with an iterative DVFS algorithm for synchronous programs. The proposed approach provides a lot more (WCRT,WCEC) tradeoffs than the previous techniques. Because the two criteria are antagonistic, we advocate that our algorithm must provide as many as possible non-dominated points (in the Pareto dominance sense) in the (WCRT,WCEC) space. To achieve this, we transform the bicriteria optimization problem into a single criterion optimization problem by turning the WCRT criterion into a constraint and by minimizing the WCEC under this constraint. By varying the WCRT constraint as in the ε-constraint method [START_REF] Haimes | On a bicriterion formulation of the problems of integrated system identification and system optimization[END_REF], we are able to produce a whole set of non-dominated Pareto optima, that is, the Pareto front. To the best of our knowledge, this is the first bicriteria (WCRT,WCEC) optimization method for synchronous programs.

The originality of our approach is to use a new iterative algorithm that explores the search space incrementally. Thanks to benchmarking with real-life programs, we have proved that our algorithm can produce a lot more points than the fixed frequency and the linearized approaches [START_REF] Xie | Compile-time Dynamic Voltage Scaling Settings: Opportunities and Limits[END_REF][START_REF] Huang | Leakage-aware intraprogram voltage scaling for embedded processors[END_REF]. The benchmarks also show that, even though the search space is exponential, the analysis time required by our algorithm is linear, thanks to a fast heuristic search. Regarding future work, we plan to investigate several avenues. First, we wish to improve the DVFS schemes found by our algorithm by a local search in the state space, for instance by lowering the frequencies in the control paths that are not the critical path. Second, we wish to accelerate our heuristic by segmenting the search space into intervals bounded by fixed frequencies, and by starting our algorithm with the fixed frequency scheme corresponding to the targeted WCRT deadline.

Figure 1 :

 1 Figure 1: An example of a TCCFG.

  Scheme2

Figure 2 :

 2 Figure 2: DVFS schemes for the running example plotted in the (WCRT,WCEC) plane.

Figure 3 :

 3 Figure 3: Overview of the proposed framework.

B2:Figure 4 :

 4 Figure 4: Converting the TCCFG in Fig 1 into four virtual TCCFGs for emulating DVFS.

Figure 5 :

 5 Figure 5: Benchmark results. Comparing three approaches in the (WCRT,WCEC) plane).

Table 1 :

 1 The execution time of the computation nodes in Fig.1at four different frequencies (Unit: cycle length at 1Mhz).

			frequency
	node 0.25Mhz 0.5Mhz 0.75Mhz 1Mhz
	B4	360	180	120	90
	B6	120	60	40	30
	B8	440	220	146.7	110

Table 2 :

 2 The estimated energy consumption of the computation nodes in Fig.1(Unit: Joule/C × 10 12 ).

			frequency
	node 0.25Mhz 0.5Mhz 0.75Mhz 1Mhz
	B4	5.62	22.5	50.62	90
	B6	1.87	7.5	16.87	30
	B8	6.87	27.5	61.87	110

Table 3 :

 3 DVFS scheme examples.

	Scheme	DVFS values (MHz)	WCRT WCEC (×10 12 )
	1	fB4=0.5, fB6=0.75, fB8=0.5	220	39.37
	2	fB4=0.75, fB6=0.5, fB8=0.5	220	58.12
	3	fB4=0.75, fB6=1, fB8=1	150	110
	4	fB4=1, fB6=0.5, fB8=1	150	110

Table 4 :

 4 Information of the benchmark programs.

	Name	LOC Threads DVFS ctrl points
	Channel Protocol	591	7	20
	Flasher	816	7	24
	Robot Sonar	962	7	18
	Cruise Controller 2302	25	60

Table 5 :

 5 Analysis time comparison between the proposed approach and the linearized approach.

		Analysis Time (s)
	Multiplier	Proposed	Linearized (Excluding profiling
			and WCRT analysis)
	1.0	3.03	1.04
	1.2	2.19	1.15
	1.4	2.25	1.12
	1.6	2.59	1.10
	1.8	2.39	1.14
	2.0	3.11	1.08
	2.2	2.54	1.14
	2.4	2.83	1.10
	2.6	3.12	1.15
	2.8	2.27	1.15
	3.0	2.00	1.11
	3.2	2.50	1.18
	3.4	1.58	1.16
	3.6	1.36	1.23
	3.8	0.45	1.24
	4.0	0.19	1.11

* This work has been partially supported by the RIPPES Inria International Lab. † Partha Roop acknowledges the German Research Foundation DFG (ME 1427/6-2, HA 4407/6-2) and the research and study leave from Auckland University.