Luca Rizzi
email: luca.rizzi@univ-grenoble-alpes.fr

Ulysse Serres
email: ulysse.serres@univ-lyon1.fr

ON THE CUT LOCUS OF FREE, STEP TWO CARNOT GROUPS

On the cut locus of free, step two Carnot groups

Introduction

We recall some basic facts about sub-Riemannian manifolds and their geodesics (see [START_REF]Introduction to Riemannian and sub-Riemannian geometry[END_REF][START_REF] Rifford | Sub-Riemannian geometry and optimal transport[END_REF] for details). Let M be a smooth manifold and D ⊂ T M be a smooth distribution with constant rank k = rank D, satisfying the Hörmander condition:

(1) Lie(Γ(D)) q = T q M, ∀q ∈ M, where the l.h.s. denotes the smallest Lie algebra generated by smooth sections of D. If g is a smooth scalar product defined on D, then (D, g) is a sub-Riemannian structure on M . A Lipschitz path γ : [0, T ] → is horizontal if γ(t) ∈ D γ(t) for a.e. t ∈ [0, T ]. For any horizontal path γ, we define its length as

(2) L(γ) := T 0 γ(t) g dt.
The length is invariant by Lipschitz reparametrization, so we can always reparametrize a horizontal curve in such a way that it has constant speed. Furthermore, we define the sub-Riemannian distance d : M × M → R as d(q, q ) := inf L(γ), where the infimum is taken over all horizontal curves that join q with q . Thanks to (1), d is finite and continuous, and so (M, d) is a locally compact length metric space.

A geodesic is a non-trivial horizontal curve γ : [0, T ] → M , with constant speed, that locally minimizes the length between its endpoints. It is maximal if it is not the restriction of a geodesic defined on a larger interval [0, T ]. The cut time of a maximal geodesic is (3) t cut (γ) := sup{t > 0 | γ| [0,t] is a minimizing geodesic} > 0.

Assuming that (M, d) is complete, the cut locus of q ∈ M is the set of cut points γ(t cut (γ)):

(4) Cut q := {γ(t cut (γ)) | γ is a maximal geodesic starting at q}.

2010 Mathematics Subject Classification. 53C17, 49J15.

The cut locus, together with the so-called abnormal geodesics, play an important role in the regularity properties of the sub-Riemannian distance [CR08, RT09, Agr09, Agr16, MM16b] and of the heat kernel of sub-Laplacians [START_REF] Barilari | Small-time heat kernel asymptotics at the sub-Riemannian cut locus[END_REF][START_REF] Barilari | On the heat diffusion for generic riemannian and sub-riemannian structures[END_REF][START_REF]Heat kernel asymptotics on sub-Riemannian manifolds with symmetries and applications to the bi-Heisenberg group[END_REF]. Its properties are quite different with respect to the Riemannian setting, for example, q is an accumulation point for Cut q , and Cut q ∪ {q} might be not closed. Moreover, while in Riemannian geometry q ∈ Cut q if and only if q is a critical value of the exponential map or there are two distinct minimizing geodesics joining q with q , this characterization is no longer true in the sub-Riemannian case due to the occurrence of abnormal geodesics.

The cut time is known explicitly for a handful of left-invariant structures on 3D Lie groups [START_REF] Boscain | Invariant Carnot-Caratheodory metrics on S 3 , SO(3), SL(2), and lens spaces[END_REF][START_REF] Sachkov | Cut locus and optimal synthesis in the sub-Riemannian problem on the group of motions of a plane[END_REF][START_REF] Butt | Cut locus and optimal synthesis in sub-riemannian problem on the lie group sh(2)[END_REF], on some Stiefel manifolds [START_REF] Autenried | Sub-Riemannian geometry of Stiefel manifolds[END_REF], and for the following nondisjoint classes of Carnot groups: contact [START_REF] Agrachev | On the Hausdorff volume in sub-Riemannian geometry[END_REF], corank 1 [START_REF] Rizzi | Measure contraction properties of Carnot groups[END_REF], corank 2 [START_REF] Barilari | On 2-step, corank 2, nilpotent sub-Riemannian metrics[END_REF], the Engel group [START_REF] Ardentov | Cut time in sub-Riemannian problem on Engel group[END_REF], the bi-Heisenberg groups [START_REF]Heat kernel asymptotics on sub-Riemannian manifolds with symmetries and applications to the bi-Heisenberg group[END_REF] (that is, corank 1 Carnot groups of dimension 5), and H-type groups [START_REF] Autenried | The sub-Riemannian cut locus of H-type groups[END_REF]. From the knowledge of the cut time, the cut locus can be computed via (4). In a few symmetric cases, this yields an elegant and compact description for the cut locus.

1.1. Summary of the results. In this paper, we focus on the hierarchy G k of free, step 2 Carnot groups of rank k 2. Their first appearance traces back to the seminal works of Gaveau [START_REF] Gaveau | Principe de moindre action, propagation de la chaleur et estimées sous elliptiques sur certains groupes nilpotents[END_REF] and Brockett [START_REF] Brockett | Control theory and singular Riemannian geometry[END_REF], and for this reason the corresponding minimization problem is called the Gaveau-Brockett problem by Liu and Sussmann [START_REF] Liu | Shortest paths for sub-Riemannian metrics on rank-two distributions[END_REF].

The groups G k have the largest sets of symmetries among all Carnot groups of the same rank and step, and their geodesics can be computed quite explicitly. These features lead to expect the existence of an explicit formula describing their cut locus. This is indeed the case for k = 2 (the Heisenberg group), where the closure of the cut locus is described succinctly as the zero set of an algebraic function. Surprisingly, such an algebraic description also exists for the free step 2 Carnot group of rank 3 -the so-called (3, 6) group -as proven independently and with different strategies in [START_REF] Myasnichenko | Nilpotent (3, 6) sub-Riemannian problem[END_REF] and [START_REF] Montanari | On the subRiemannian cut locus in a model of free two-step Carnot group[END_REF].

Extrapolating from the known results in these cases, in [START_REF] Myasnichenko | Nilpotent (3, 6) sub-Riemannian problem[END_REF][START_REF]Nilpotent (n, n(n + 1)/2) sub-Riemannian problem[END_REF][START_REF] Montanari | On the subRiemannian cut locus in a model of free two-step Carnot group[END_REF], the authors conjectured precise algebraic formulas for the cut loci of G k , for all k 2. In this note, we disprove these conjectures. In particular, we exhibit an explicit hierarchy of sets C k ⊂ G k of cut points (see Definition 5), which coincide with the corresponding cut loci for k = 2, 3, but are strictly larger than the conjectured ones for k 4. While the previously conjectured cut loci were, respectively, smooth algebraic sets of codimension Θ(k 2 ) [START_REF] Myasnichenko | Nilpotent (3, 6) sub-Riemannian problem[END_REF][START_REF]Nilpotent (n, n(n + 1)/2) sub-Riemannian problem[END_REF] and algebraic sets of codimension Θ(k) [START_REF] Montanari | On the subRiemannian cut locus in a model of free two-step Carnot group[END_REF], the set C k is semi-algebraic of codimension 2, for all k (see Theorem 11).

Furthermore, we study the relation between the cut locus and the so-called abnormal set. In the low dimensional cases, it is known that

(5) Abn 0 (G k ) = Cut 0 (G k ) \ Cut 0 (G k ), k = 2, 3.
Instead, for each k 4, we show that the cut locus always intersects the abnormal set, and that there are plenty of abnormal geodesics with finite cut time (see Proposition 12).

Finally, and as a straightforward consequence of our results, we obtain an explicit lower bound for the small time heat kernel asymptotics at the points of C k (see Theorem 13).

Free step 2 Carnot groups of rank k

We refer to [START_REF] Montgomery | A tour of subriemannian geometries, their geodesics and applications[END_REF][START_REF] Jean | Control of nonholonomic systems: from sub-Riemannian geometry to motion planning[END_REF], for the definition of Carnot group. Here, we only deal with the specific free, step 2 case. Let G k := R k ⊕ ∧ 2 R k . We identify ∧ 2 R k with the vector space of skew-symmetric real matrices, that is v ∧ w = vw * -wv * for v, w ∈ R k . We denote points (x, Y ) ∈ G k , where x ∈ R k and Y is a skew-symmetric matrix. The free, step 2 Carnot group of rank k is the sub-Riemannian structure on G k generated by the set of global orthonormal vector fields: (6)

X i := ∂ x i - 1 2 1 <m k (e i ∧ x) m ∂ Y m , i = 1, . . . , k,
where {e 1 , . . . , e n } is the standard basis of R k . More precisely, the horizontal distribution is defined by D := span{X 1 , . . . , X k } and the sub-Riemannian metric by g(X i , X j ) = δ ij . For all i < j, we have [X i , X j ] = ∂ Y ij . In particular, the vector fields (6) generate the free, nilpotent Lie algebra of step 2 with k generators:

(7) g = g 1 ⊕ g 2 , where g 1 = span{X 1 , . . . , X k }, g 2 = span{∂ Y ij } i<j .
There exists a unique Lie group structure on G k such that the vector fields X i are leftinvariant, given by the polynomial product law

(8) (x, Y ) (x , Y ) = x + x , Y + Y + 1 2 x ∧ x .
Thus, (G k , ) is a connected, simply connected Lie group of dimension k(k + 1)/2, such that its Lie algebra g of left-invariant vector fields is isomorphic to the free, nilpotent, stratified Lie algebra of step 2 with k generators, and such that the first stratum g 1 is equipped with a left-invariant scalar product. Any Lie group (G k , ) with the same properties is isomorphic to (G k , ) and they carry isometric sub-Riemannian structures.

Carnot groups are equipped with a one-parameter family of dilations. For G k , it is given by δ ε (x, Y ) := (εx, ε 2 Y ), for ε > 0. As a consequence of this fact, the metric spaces (G k , d) are complete, and there exists a minimizing geodesic joining any given pair of points.

Example 1. The case k = 2 is the well-known Heisenberg group. Indeed, we can identify (x, Y ) ∈ R 2 ⊕ ∧ 2 R 2 with (x, z) ∈ R 2 ⊕ R, so that the generating vector fields (6) read (9)

X 1 = ∂ x 1 - x 2 2 ∂ z , X 2 = ∂ x 2 + x 1 2 ∂ z .
Example 2. The case k = 3 can be dealt with by identifying (x, Y ) ∈ R 3 ⊕ ∧ 2 R 3 with (x, t) ∈ R 3 ⊕ R 3 . More precisely, any 3 × 3 skew-symmetric matrix can be written as Y = v ∧ w (in a non-unique way), and is identified with the cross product t = v × w. Under this identification, the tautological action of Y on R 3 reads (10)

Y x = (v ∧ w)x = x × (v × w) = x × t, ∀x ∈ R 3 ,
and the generating vector fields (6) are

(11) X 1 = ∂ x 1 + x 3 2 ∂ t 2 - x 2 2 ∂ t 3 , X 2 = ∂ x 2 + x 1 2 ∂ t 3 - x 3 2 ∂ t 1 , X 3 = ∂ x 3 + x 2 2 ∂ t 1 - x 1 2 ∂ t 2 .
2.1. Geodesics. In this case, a Lipschitz path γ :

[0, T ] → G k is horizontal if there is a control u ∈ L ∞ ([0, T ], R k ) such that, for almost every t ∈ [0, T ], we have (12) γ(t) = k i=1 u i (t)X i (γ(t)).
The standard method to solve the length minimization problem is based on the Pontryagin maximum principle. In the step 2 case, thanks to Goh condition, one can rule out all the so-called strictly abnormal curves. This analysis (see e.g. [LS95, Appendix A]) yields that geodesics are all the horizontal curves γ(t) = (x(t), Y (t)) with control (13)

u(t) = e -tΩ p, (p, Ω) ∈ R k ⊕ ∧ 2 R k .
In particular, all geodesics are real-analytic curves and can be extended on the maximal interval [0, +∞). The pair (p, Ω) is also referred to as the initial covector in the Hamiltonian formalism. Taking into account the explicit fields (6), for any fixed (p, Ω), the actual geodesic can be recovered integrating

(14) ẋ(t) = u(t), Ẏ (t) = 1 2 x(t) ∧ u(t).
Remark 3. We point out a useful observation from [GJZ09, Prop. 3], which follows by putting Ω in real normal form. For any control u(t) as in (13), there is a unique Ω ∈ ∧ 2 R k with all simple non-zero eigenvalues and largest possible kernel such that e -tΩ p = e -t Ωp.

Most of the progress in the study of free Carnot groups of low dimension (k = 2, 3) is due to the fact that ( 14) can be explicitly integrated in terms of trigonometric functions. Geodesics for the case k = 2, corresponding to the Heisenberg group, are well known. For the case k = 3, explicit formulas appeared first in [START_REF] Myasnichenko | Nilpotent (3, 6) sub-Riemannian problem[END_REF]. An explicit integration of the geodesic flow for the general case G k appeared first in [Mya06, Thm. 4.1]. However, if the latter were correct, it would imply that the geodesic flow can only reach points (x, Y ) with rank(Y ) = 2. To our best knowledge, the only complete and correct integration of the geodesic flow for general G k appears in [START_REF] Monroy-Pérez | The step-2 nilpotent (n, n(n+1)/2) sub-Riemannian geometry[END_REF], in terms of the spectral projectors of Ω, together with explicit worked out examples for all k 5.

The conjectured cut loci

By the left-invariance of the sub-Riemannian structure on G k , one can recover Cut q (G k ) for any q ∈ G k by left translation of Cut 0 (G k ), where 0 = (0, 0) is the identity of G k . In the following, with the term "cut locus", we will refer to Cut 0 (G k ).

The cut locus of G 2 (the Heisenberg group) is well known, and consists in the set of points (0, Y ) ∈ R 2 ⊕ ∧ 2 R 2 , with Y = 0. The computation of the cut time and cut locus for the case k = 3 is much harder, and requires a careful manipulation of geodesic equations and symmetries. This has been done independently in [START_REF] Myasnichenko | Nilpotent (3, 6) sub-Riemannian problem[END_REF] and [START_REF] Montanari | On the subRiemannian cut locus in a model of free two-step Carnot group[END_REF]. There, it was proved that, for the step 2 free Carnot groups of rank k = 2, 3, one has

(15) Cut 0 (G k ) = {(x, Y ) | Y = v ∧ w = 0, Y x = 0}, k = 2, 3.
The characterization in (15) being dimension-free, the authors were naturally led to two closely related conjectures.

Conjecture 1 [START_REF] Myasnichenko | Nilpotent (3, 6) sub-Riemannian problem[END_REF][START_REF]Nilpotent (n, n(n + 1)/2) sub-Riemannian problem[END_REF]). The cut locus of step 2, free Carnot group of rank k is

(16) P k := {(x, Y ) | Y = v ∧ w = 0, Y x = 0}.

Conjecture 2 ([MM16a]

). The cut locus of step 2, free Carnot group of rank k is

(17) Σ k := {(x, Y ) | Y = 0, Y x = 0}.
Clearly, Σ k = P k for k = 2, 3, but P k Σ k for k 4. These conjectures imply analogous claims on the size of the cut loci, in particular codim(P k ) = Θ(k 2 ) and codim(Σ k ) = Θ(k), for large k (see the forthcoming Proposition 4).

Semi-algebraic sets.

A set A ⊂ R n is semi-algebraic if it is the result of a finite number of unions and intersections of sets of the form {f = 0}, {g > 0}, where f , g are polynomials on R n . We recall some of their basic properties, referring to [START_REF] Bochnak | Real algebraic geometry[END_REF] for details. If A is smooth, its dimension as a semi-algebraic set [BCR98, Def. 2.8.1] is equal to its dimension as a smooth manifold. Moreover, if A is the finite union of semi-algebraic sets A 1 , . . . , A p , then dim(A) = max i dim(A i ). Finally, if A is a semi-algebraic set, then its closure, its interior and any Cartesian projection of A are semi-algebraic sets.

Proposition 4. P k ⊆ G k is a smooth semi-algebraic set of codimension (k 2 -5k + 10)/2. Σ k ⊆ G k is a semi-algebraic set of codimension 2 k/2 .
Proof. It is clear that both sets are semi-algebraic. Notice also that

(18) Σ k = k/2 r=1 Σ 2r k , Σ 2r k = {(x, Y ) | Y x = 0, rank(Y ) = 2r}.
Indeed P k = Σ 2 k . We now prove that each Σ 2r k is a smooth manifold and we compute its codimension. The set S 2r ⊂ ∧ 2 R k of rank 2r skew-symmetric matrices is a smooth submanifold of codimension dim(so(k

-2r)) = (k -2r)(k -2r -1)/2. Let φ 2r : R k × S 2r → R k be the smooth map φ 2r (x, Y ) = Y x. Then, Σ 2r k = φ -1 2r (
0). We claim that φ 2r has rank 2r. This claim concludes the proof since φ -1 2r (0) is then a codimension 2r submanifold of a codimension (k-2r

)(k-2r-1)/2 submanifold of G k , that is a submanifold of codimension (k -2r)(k -2r -1)/2+2r = (k 2 -(4r +1)k +2r(2r +3))/2.
To prove the claim, we drop the subscript from the notation φ 2r , and we assume, without loss of generality, that (x,

Y ) ∈ Σ 2r k is of the form (19) Y = Ȳ 0 k-2r , x = 0 x 0 , rank( Ȳ ) = 2r, x 0 ∈ R k-2r .
The smooth variations

γ i (ε) = (x + εe i , Y ) ∈ R k × S 2r , for i = 1, . . . , 2r, are such that (20) ∂ t φ(γ i )(0) = Ȳ e i .
Thus, the image of the differential of φ at (x, Y ) contains span{e 1 , . . . , e 2r }, and rank(φ) 2r. Suppose that rank(φ) > 2r. Therefore, there exists

γ(t) = (x+tx , Y +tY ) ∈ R k ×S 2r , such that ∂ t φ(γ(0)) and span{e 1 , . . . , e 2r } are independent. Splitting Y as in (19), (21)
Y = Y 1 Y 2 Y 3 Y 4 , we obtain ∂ t φ(γ(0)) = * Y 4 x 0 ,
where the latter is independent of span{e 1 , . . . , e 2r }. But this means that, for sufficiently small t, we have rank(Y + tY ) > 2r, which contradicts the fact that γ ∈ R k × S 2r .

A larger set of cut points

The orthogonal group O(k) acts smoothly on G k . In particular, any M ∈ O(k) induces an isometry of Carnot groups ρ M :

G k → G k , given by ρ M (x, Y ) = (M x, M Y M * ). Definition 5. We define C k as the set of points (x, Y ) ∈ R k ⊕ ∧ 2 R k , with Y = 0, such that there exists a non-trivial M ∈ O(k) stabilizing (x, Y ), i.e. (22) M x = x, M Y M * = Y,
and such that M | ker Y = 1.

For the Heisenberg group,

C 2 = Σ 2 = P 2 = {(0, Y ) | Y = 0}. For rank k = 3, we can write Y = v ∧ w for some (non-unique) pair v, w ∈ R k . In this case, any M ∈ O(3) such that M Y M * = Y is
an orthogonal transformation in the plane span{v, w} = (ker Y ) ⊥ . By the last condition, such a rotation is non-trivial. Thus, the condition M x = x implies that x ⊥ span{v, w}, and C 3 = Σ 3 = P 3 . As it will be evident from the forthcoming discussion about normal forms, this is no longer true for k 4. We now clarify the shape of C k .

Y =      α 1 J 2m 1 . . . α J 2m 0 m 0      , J 2m = 0 1 m -1 m 0 .
Accordingly, we write x = (x 1 , . . . , x , x 0 ), with x i ∈ R 2m i for i = 1, . . . , , and

x 0 ∈ R m 0 . If we insist that (x, Y ) ∈ C k , the conditions M | ker Y = 1 and M Y M * = Y imply that (24) M =      M 1 . . . M 1 m 0      , M i ∈ O(2m i ) ∩ Sp(2m i ) U(m i ).
Thus, the conditions M x = x, and M = 1 yield that (x, Y ) ∈ C k if and only if, for some i ∈ {1, . . . , }, either x i = 0 or m i > 1. In the first case, there is no further restriction on M i , while in the second we must have

M i x i = x i . Thus, if (x, Y ) ∈ C k , the set (25) M (x, Y ) := {M ∈ O(k) satisfying the property of Definition 5} ∪ {1 k }, is a non-trivial Lie subgroup of O(k) of dimension (26) dim M (x, Y ) = i∈I 0 dim U(m i ) + i / ∈I 0 dim U(m i -1),
where I 0 ⊆ {1, . . . , } is the set of indices i such that x i = 0. With this notation, we have the following.

Proposition 6. For any geodesic γ joining the origin with (x, Y ) ∈ C k , there exists a nonempty family of mutually distinct geodesics with the same endpoints and length, obtained by the action of ρ M on γ, with M ∈ M (x, Y ). The dimension of such a family is

(27) N (x, Y ) = i∈I 0 m 2 i + i / ∈I 0 (m i -1) 2 > 0.
Proof. For all M ∈ M (x, Y ), the curves ρ M (γ) are geodesics joining the same endpoints and with the same length of γ. By contradiction, assume that ρ M 1 (γ(t)) = ρ M 2 (γ(t)) for all t ∈ [0, T ] and for

M 1 = M 2 ∈ M (x, Y ). In particular, 1 = M * 2 M 1 ∈ M (x, Y ). Let E R n be the eigenspace of M *
2 M 1 corresponding to the eigenvalue 1. Since M 1 x(t) = M 2 x(t), and ẋ(t) = u(t), both u(t), x(t) ∈ E for all t ∈ [0, 1]. The geodesic equations ( 14) imply ( 28)

Y (t) ∈ ∧ 2 E ∧ 2 R k . Let v ⊥ E, with v = 0. By (28), v ∈ ker Y . By Definition 5, however, M * 2 M 1 v = v, hence v ∈ E and v = 0.
This is a contradiction, so the ρ M (γ) are mutually distinct. The dimension of such a family is N (x, Y ) = dim M (x, Y ) given by (26).

Remark 7. By construction, the family in Proposition 6 is certainly composed of distinct geodesics obtained from a given one by isometric transformations. This might not exhaust all geodesics joining the origin with a given point (x, Y ) ∈ C k . For example, in corank 1 Carnot groups there are continuous families of non-isometric geodesics joining the origin with a given (non-generic) point, all having the same length [START_REF] Lerario | How many geodesics join two points on a contact sub-Riemannian manifold?[END_REF].

Proposition 8. For all k 2, the set C k ⊂ G k is semi-algebraic of codimension 2.
M | ker Y = 1, that is M W = W .
Thus, C k is the projection on the first two factors of the set V k ⊂ R k ×∧ 2 R k ×O(k)×R k×k , defined by the following formula:

(29) V k := {(x, Y, M, W ) | M x = x, M Y M * = Y, rank(W ) = k -rank(Y ), Y W = 0, M W = W } .
The rank conditions being semi-algebraic, V k is a semi-algebraic set, and so is C k .

To compute the dimension of C k , we follow a strategy inspired by [AGL15, App. A]. We decompose C k as the disjoint union of a finite number of smooth semi-algebraic sets, each one containing points (x, Y ) with the same "type" of normal form as above:

(30) C k = k/2 r=1   r-1 =1 m 1 +•••+m =r C 2r|m 1 ,...,m k I 0 =∅ C 2r|1,...,1|I 0 k   .
Each component in (30) is labeled by the rank 2r of Y , and the multiplicities m 1 , . . . , m of its non-zero distinct eigenvalues, ordered according to their absolute values 0 < α 1 < • • • < α . In the special case of all simple eigenvalues = r, that is m 1 = • • • = m r = 1, we further decompose according to the set I 0 ⊆ {1, . . . , r} such that x i = 0 for i ∈ I 0 , and x i = 0 for i / ∈ I 0 . This decomposition exhausts all points of C k . We now prove that the component C 2r|1,...,1|I 0 k , corresponding to the case in which rank(Y ) = 2r = 2 k/2 is maximal, with simple eigenvalues, and where I 0 = {i} is a singleton, is smooth of codimension 2. A similar argument, which we omit, proves that all the other components in the decomposition are smooth, and with larger codimension.

Assume, for simplicity, that k is even and that

I 0 = {1}. Let R n * = R n \ {0} and R n ord denote the open subset of ordered, distinct n-tuples 0 < α 1 < • • • < α n . Consider the map Φ : R k/2 ord × R k/2-1 * × O(k)/U(1) → C 2k|1,...,1|I 0 k given by (31) Φ(α 1 , . . . , α k/2 , b 2 , . . . , b k/2 , M ) = ρ M           0 b 2 . . . b k/2      ,      α 1 J 2 α 2 J 2 . . . α k/2 J 2          
, where O(k)/U(1) is the homogeneous space (a smooth manifold) of left cosets M U , where U ∈ U(1) ⊂ O(k) is the subgroup that stabilizes the given normal form of (x, Y ), that is

(32) U =      M 1 1 2 . . . 1 2      , M 1 ∈ O(2) ∪ Sp(2) U(1).
The map Φ is a bijection. This gives C 2k|1,...,1|I 0 k a smooth structure, with dimension

(33) k 2 + k 2 -1 + dim O(k)/U(1) = k(k + 1) 2 -2.
The same argument yields the same formula for odd k. Thus, since dim(G k ) = k(k + 1)/2, the codimension of the larger component of C k is 2. Notice that, for k = 2, 3, there is only one component in (30), hence in these cases C k is a smooth semi-algebraic set.

The next statement holds true, in particular, for all step 2 Carnot groups.

Lemma 9. Let (D, g) be a complete sub-Riemannian structure on M . Assume that all minimizing curves are real-analytic. If there exist two distinct minimizing geodesics joining q with q , then q ∈ Cut q .

Proof. Let γ i : [0, T ] → M , i = 1, 2, be distinct minimizing geodesics joining q with q . In local coordinates, ∃n ∈ N such that their n-th derivatives at T are different. Assume that we can extend γ 1 to a minimizing geodesic on the interval [0, T + ε], with ε > 0. Then, the curve γ : [0, T + ε] → M such that γ| [0,T ] = γ 2 , and γ| [T,T +ε] = γ 1 is minimizing on the interval [0, T + ε], but not real-analytic. This is a contradiction, then q ∈ Cut q .

Remark 10. The conclusion of Lemma 9 holds true, although with a different proof, if the real-analytic hypothesis is replaced by the assumption that there are no abnormal minimizing curves. Otherwise, one cannot a priori exclude the case of "branching" minimizing curves. More precisely, one can prove that the segment γ| [T,T +ε] defined as in the proof of Lemma 9 must be an abnormal minimizing curve.

It is now clear that the set C k of Definition 5 is a set of cut points, which for k 4 is strictly larger than the sets P k and Σ k of Conjectures 1 and 2.

Theorem 11. The following chain of inclusions holds true, for all k 4: (34)

P k Σ k C k ⊆ Cut 0 (G k ). In particular, if Cut 0 (G k ) is semi-algebraic, codim(Cut 0 (G k )) 2, for all k.
Proof. The first strict inclusion is trivial. The second is a consequence of the normal forms in Section 4.1. Notice that all minimizing curves of G k are real-analytic. Hence, the third inclusion follows from Proposition 6, applied to minimizing geodesics, and Lemma 9. The claim on the codimension follows from Proposition 8.

4.2. Relation with the abnormal set. Let End : L ∞ ([0, 1], R k ) → G k be the end-point map, which associates with a given control u(t) the end-point of the corresponding horizontal curve (12). Trajectories corresponding to critical points of End are called singular curves, and the set Abn 0 (G k ) of critical values of End is the so-called abnormal set. Determining the "size" of Abn 0 (G k ) on a general sub-Riemannian manifold is one of the hard open problems in the field [START_REF]Some open problems[END_REF]. A computation using the formula for the differential of the end-point map and the chronological calculus (see [START_REF]Introduction to Riemannian and sub-Riemannian geometry[END_REF]) yields

(35) Abn 0 (G k ) = W ∈Gr(k,k-2) W ⊕ ∧ 2 W, ∀k 2,
where Gr(k, k -2) is the Grassmannian of codimension 2 subspaces of R k . We refer to [DMO + 15] for a proof of (35) not resorting to chronological calculus, and several alternative characterizations. We point out that, if dim(ker Y ) 3, then (x, Y ) ∈ Abn 0 (G k ).

In the low dimensional cases, it turns out that the boundary of the cut locus coincides with the abnormal set:

(36) Abn 0 (G k ) = Cut 0 (G k ) \ Cut 0 (G k ), k = 2, 3.
For k 4, formula (36) is no longer true. More precisely, we have the following result, yielding also the first example of abnormal geodesic with finite cut time.

Proposition 12. For all k 2, we have

(37) C k \ C k ⊆ Abn 0 (G k ) C k .
For each k 4, the first inclusion is strict, and thus Cut 0 (G k ) ∩ Abn 0 (G k ) = ∅. Moreover, there exist abnormal geodesics with finite cut time.

Proof. Let (x, Y ) ∈ C k \ C k .
By the discussion in Section 4.1, we assume that

(38) x =      x 1 . . . x x 0      , Y =      α 1 J 2 . . . α J 2 0 m 0      , 0 = x i ∈ R 2 , i = 1, . . . , .
Let (x n , Y n ) ∈ C k be a sequence converging to (x, Y ) as n → +∞. Since the nonzero eigenvalues of Y are simple, the same holds for the non-zero eigenvalues of Y n . Consequently, the associated real projectors π n i , for i = 1, . . . , , are continuous [Kat95, Chap. 2, Thm. 5.1], and so are the projections π n ⊥ := 1 -i=1 π n i . Observe that lim n π n ⊥ = π 0 is the projection on ker Y , but π n ⊥ is not necessarily the projection on ker Y n (a pair of eigenvalues ±iα n = 0 of Y n could coalesce to zero for n → +∞). In particular, we stress that m 0 = dim(ker Y ) 2, otherwise no such a sequence could exist.

If x 0 = 0, then (x, Y ) ∈ Abn 0 (G k ) by the characterization (35) of the abnormal set. Thus, let x 0 = 0. By our assumption on the sequence (x n , Y n ), we have π n i (x n ) = 0, for all i ∈ {1, . . . , , ⊥}, and large n. Since (x n , Y n ) ∈ C k , and by the characterization of Section 4.1, we have m 0 = dim(ker Y ) 3. In this case, (x, Y ) ∈ Abn 0 (G k ), by (35). This proves that

C k \ C k ⊆ Abn 0 (G k ).
We now prove that Abn 0 (G k ) C k . By the characterization of (35), let (x, Y ) ∈ W ⊕ ∧ 2 W , where W ⊂ R k is a codimension 2 subspace. Up to an orthogonal transformation, we can assume that

W = {(x 1 , . . . , x k-2 , 0, 0) * ∈ R k } R k-2 , that is (39) x = x 1 0 , Y = Y 1 0 2 , x 1 ∈ W, Y 1 ∈ ∧ 2 W.
Let (x n , Y n ) be the sequence with

x n = x,
x = 0, Y = αJ 2 0 k-2 , α > 0.
It is not hard to prove that the geodesics γ :

[0, 1] → G k with control u(t) = e -tΩ p, with (41) p
= p 0 0 , Ω = 2πJ 2 0 k-2 , p 0 ∈ R 2 , |p 0 | 2 = 4πα,
joins the origin with the given point and is abnormal. The formula for the distance in Appendix A yields d((0, 0), (x, Y )) 2 = 4πα = L(γ) 2 , yielding the minimality property of γ. By Theorem 11, all elements of C k are cut points, therefore t cut (γ) = 1.

Small time heat kernel asymptotics

Any sub-Riemannian manifold equipped with a smooth measure supports an intrinsic hypoelliptic operator, the sub-Laplacian, playing the role of the Laplace-Beltrami operator of Riemannian geometry. In the case of Carnot groups equipped with a left-invariant measure -and in particular for G k -the sub-Laplacian is the "sum of squares" given by (42

) ∆ := k i=1 X 2 i .
Let K t (q , q) denote the heat kernel of ∆, that is the fundamental solution of the heat equation ∂ t ψ = ∆ψ. Its existence and smoothness for t > 0 on complete sub-Riemannian manifolds is classical, see [START_REF] Strichartz | Sub-Riemannian geometry[END_REF]. A well known result due to Léandre [START_REF] Léandre | Majoration en temps petit de la densité d'une diffusion dégénérée[END_REF][START_REF]Minoration en temps petit de la densité d'une diffusion dégénérée[END_REF] yields the relation between the heat kernel and the sub-Riemannian distance:

(43) lim

t→0 + 4t log K t (q, q ) = -d(q, q ) 2 .
Furthermore, the singularity of the heat kernel for t → 0 depends on the "amount" of minimizing geodesics joining q with q (see [START_REF] Barilari | Small-time heat kernel asymptotics at the sub-Riemannian cut locus[END_REF] and references therein).

For what concerns Carnot groups of step 2, an integral representation for K t (q) := K t (q, 0) has been obtained in [START_REF] Beals | The Green function of model step two hypoelliptic operators and the analysis of certain tangential Cauchy Riemann complexes[END_REF]Thm. 4] (the general case can be easily recovered by left-invariance). For G k , and in our notation, such a formula reads (44)

K t (x, Y ) = 2 (4πt) k 2 /2 ∧ 2 R k det B sin B cos Tr(BY * ) 2t e -x * B cot(B)x 4t dB,
where dB is the Lebesgue measure of ∧ 2 R k R k(k-1)/2 , and analytic functions of a matrix are defined by their power series. As an application of our results, combined with the ones in [START_REF] Barilari | Small-time heat kernel asymptotics at the sub-Riemannian cut locus[END_REF], we mention the following small time asymptotic result.

Theorem 13. Let (x, Y ) ∈ C k \ Abn 0 (G k ).
Then there exists a constant C > 0 such that,

(45) K t (x, Y ) C + O(t) t (n+N (x,Y ))/2 e -d 2 0 (x,Y )/4t , as t → 0,
where d 0 is the sub-Riemannian distance from the origin, n = dim(G k ) = k(k + 1)/2, and N (x, Y ) is given by Proposition 6 in terms of the normal form of the point (x, Y ).

Proof. Let (p, Ω) be the initial covector corresponding to a minimizing curve with control u(t) = e -Ωt p, joining the origin with (x, Y ) ∈ C k . By Proposition 6, there exists a continuous family of distinct minimizing geodesics with the same endpoints. The set of initial covectors of this family is Since (x, Y ) / ∈ Abn 0 (G k ), then all geodesic γ : [0, 1] → M of this family are normal. The sub-Riemannian structure on G k being real-analytic, then these geodesics are also strongly normal (that is, any restriction γ| [0,t] , for t ∈ (0, 1], is not abnormal, see [ABR13, Prop. 3.12]). A straightforward modification of [BBN16, Thm. 16] yields the result.

(46) O = {(M p, M ΩM * ) | M ∈ M (x, Y )} ⊂ T * (x,Y ) G k G k , where M (x, Y )

Closing remarks and open problems

6.1. Historical remarks. The seminal papers by Gaveau [START_REF] Gaveau | Principe de moindre action, propagation de la chaleur et estimées sous elliptiques sur certains groupes nilpotents[END_REF] and Brockett [START_REF] Brockett | Control theory and singular Riemannian geometry[END_REF] mentioned in the introduction contain general statements on the minimization problem on G k some of which are left unproven. These results might be useful to determine the true shape of Cut 0 (G k ), and we list them here, commenting their implications.

• The third (unproven) statement of [Bro82, Thm. 2] asserts that any two optimal controls u i (t), i = 1, 2, such that the associated geodesics γ i join (0, 0) with (x, Y ) are related by u 1 (t) = θu 2 (t) for some orthogonal matrix such that θY θ * = Y .

We were not able neither to prove nor to disprove this claim.

= Cut 0 (G k ).
Proof. For a general sub-Riemannian manifold, let II q be the set of points reached by at least two distinct minimizing geodesics, issued from q. By [ABB16, Thm. 8.54], we have, (47) Cut q ⊆ II q ∪ (Cut q ∩ Abn q ).

Assume that the third statement of [Bro82, Thm. 2] is true. In particular, if (x, Y ) ∈ II 0 (G k ), there exists a non-trivial θ ∈ O(k) that stabilizes (x, Y ). Then, the discussion about normal forms of points of C k in Section 4.1 and the characterization (35) imply that II 0 (G k ) ⊆ C k ∪ Abn 0 (C k ). Taking the closure, and using Proposition 12, II 0 (G k ) ⊆ C k . This, together with (47), imply Cut 0 (G k ) ⊆ C k ⊆ Cut 0 (G k ), concluding the proof.

6.2. The definition of cut locus. The sub-Riemannian cut locus is sometimes defined in the literature as the set of points joined by two distinct minimizing geodesics, see e.g. [START_REF] Autenried | Sub-Riemannian geometry of Stiefel manifolds[END_REF][START_REF] Autenried | The sub-Riemannian cut locus of H-type groups[END_REF][START_REF]Tangent Hyperplanes to Subriemannian Balls[END_REF]. While this definition might be useful for some statements, it does not correspond to the classical one in the Riemannian case.

6.3. Codimension. Generically, one would expect that the cut locus of a complete sub-Riemannian structure has Hausdorff codimension 1 (with respect to some, and then any, auxiliary Riemannian metric). This is indeed the case if there are no abnormal minimizing curves [START_REF] Rifford | On the stabilization problem for nonholonomic distributions[END_REF]. However, Carnot groups are certainly not "generic", and in fact codim(Cut 0 (G k )) = 2 for k = 2, 3. By Theorem 11, codim(Cut 0 (G k )) 2 for k 4. It is reasonable to expect that codim(Cut 0 (G k )) = 2 for all k, but this problem remains open.

6.4. Symmetries and geodesic equations. Note that one can prove that Cut 0 (G 2 ) = C 2 with a purely symmetry argument. Indeed, the invariance of the cut locus by dilations and isometries implies that if (x, z) ∈ Cut 0 (G 2 ) and x = 0, then (εM x, ε 2 z) ∈ Cut 0 (G 2 ), for all ε > 0 and M ∈ O(2). This parabola separates G 2 into two disconnected regions. Since no minimizing path can cross the parabola at an intermediate time, this yields a contradiction, and implies that Cut 0 (G 2 ) ⊆ {(x, 0) | x ∈ R 2 * } = C 2 . By Theorem 11 -the proof of which, we stress, never uses the integrated geodesic equations -we have C 2 ⊆ Cut 0 (G 2 ), which proves the claim.

Unfortunately, this argument does not carry on for rank k 3. In order to prove whether C k = Cut 0 (G k ), for all k, it seems that a fine analysis of the integrated geodesic equations is required (as for the cases k = 2, 3). Any attempt in this direction had proved inconclusive, due to the very complicated structure of the geodesics for k 4.

Acknowledgments. This research has been supported by the European Research Council, ERC StG 2009 "GeCoMethods", contract n. 239748. The first author was partially supported by the Grant ANR-15-CE40-0018 of the ANR, by the iCODE institute (research project of the Idex Paris-Saclay), and by the SMAI project "BOUM". This research benefited from the support of the "FMJH Program Gaspard Monge in optimization and operation research" and from the support to this program from EDF. The second author was partially supported by the Grant ANR-12-BS03-0005 LIMICOS of the ANR.

The first author wishes to thank A. Gentile for helpful discussions on semi-algebraic sets, U. Boscain for hinting at the implications of Proposition 6 to heat kernel asymptotics, and Y. Sachkov for carefully reading a preliminary version of the paper.

4. 1 .

 1 Normal forms. Let (x, Y ) ∈ G k . Let 0 < α 1 < • • • < α be the absolute values of the non-zero distinct eigenvalues of Y , each having multiplicity m i , and m 0 = dim(ker Y ). Up to an isometry ρ O , we can assume that Y has the real normal form (23)

Proof.

  By Definition 5, (x, Y ) ∈ C k if and only if there exist M ∈ O(k) stabilizing (x, Y ), a matrix W ∈ R k×k whose columns generate ker Y , and

  is the Lie subgroup of isometries discussed in Section 4.1. In particular, O is an N (x, Y )-dimensional manifold, and C k are critical values of the sub-Riemannian exponential map exp 0 : T * 0 G k → G k , which sends (p, Ω) to the point at time 1 of the curve with control u(t) = e -tΩ p. Moreover, dim(ker d λ exp 0 ) N (x, Y ), for all points λ ∈ O.

  and Y n obtained from Y by replacing the lower right block 0 2 with 1 n J 2 . Then (x n , Y n ) ∈ C k for all n, and the sequence converges to the given abnormal point (x, Y ). This proves Abn 0 (G k ) ⊆ C k . The inclusion is strict, for all k 2, since (x, Y ) ∈ C k does not belong to Abn 0 (G k ) if Y has maximal rank. To conclude the proof, for each k 4, we give an example of cut point (x, Y ) ∈ C k ∩ Abn 0 (G k ), reached by an abnormal minimizing geodesic. Let

	(40)

  However, if true, it would imply that C k = Cut 0 (G k ) (see Proposition 14); • By Proposition 6, points of C k are critical values of the sub-Riemannian exponential map (see [ABB16,Rif14]). Consequently, any point (x, Y ) ∈ C k with x ⊥ ker Y is a counter-example to the "only if" part of the fourth claim of [Bro82, Thm. 2]; • For completeness, we recall that, according to [LS95, Appendix A], the result in [Gav77, Thm. 1, p. 133] is false. If the third statement of [Bro82, Thm. 2] is true, then C k

	Proposition 14.

Thus, the non-zero eigenvalues of Ω are of the form ±i2πφ, with φ ∈ N. By Remark 3, and up to an orthogonal transformation, we assume that Ω = diag(2πφ 1 J 2 , . . . , 2πφ J 2 , 0 k-2 ), with all simple eigenvalues, 2 = rank(Ω), and with distinct φ i ∈ N. We split accordingly p = (p 1 , . . . , p , p 0 ), with p i ∈ R 2 for i = 1, . . . , and p 0 ∈ R k-2 . Under these assumptions, it is not hard to integrate the vertical part of the geodesic equations (14). We obtain

Then, 2 = 2r, and |p j | 2 = 4πφ j α j for all j = 1, . . . , r. The squared length of γ is

The minimum of the length is obtained when φ j = r -j + 1, for all j = 1, . . . , r.


Bernard Prade

Aurélien Houard

Jean Larour

Michel Pellet

A Mysyrowicz
email: andre.mysyrowicz@ensta.fr.

Transfer of microwave energy along a filament plasma column in air

ouverte pluridisciplinaire HAL, est destinée au dépôt

A sufficiently intense femtosecond laser pulse can propagate in air in the form of a filamentary pulse that maintains a high intensity over several hundreds of meters [1]. This non linear propagation regime results from a complex interplay between diffraction, the optical Kerr effect, high field ionization, plasma defocusing and other nonlinear effects. A column of weakly ionized plasma, of free electron density n e ~ 10 17 cm -3 , is left in the wake of the filamentary pulse. This plasma column can lead to many applications. It is the source of a forward oriented conical emission of pulsed THz radiation [2]. It can form a low-jitter fast spark gap [3]. It can also trigger and guide meter long electric discharges [4], allowing contactless transfer of large currents [5]. In turn, long-lived guided discharges can act as virtual antennas [6] or could improve the aerodynamics of fast moving object through atmosphere [7].

In this article we use an oscillating RF electric field to excite an electromagnetic TM wave with maximum electric field amplitude located at the interface between a small section of a multi-filament plasma column and the non ionized external medium [8][9]. This surface wave propagates out of the excitation region, allowing the transfer, over a significant distance, of microwave energy to the plasma in the form of longitudinal oscillations. The propagation distance should reach several meters in the case of lower frequency electromagnetic waves. Guided propagation of microwave radiation by filaments using a different approach has been discussed previously both theoretically [10][11][12] and experimentally [13][14][15]. We suggest that a self-sustained plasma antenna could be realized if sufficient microwave or VHF power is injected.

Fig. 1 shows the experimental arrangement: the plasma channel is created by filamentation of intense laser pulses at 800 nm. Such laser pulses with 200-mJ energy and adjustable duration between 50 fs and 1 ps are delivered by a multiterawatt CPA laser source (Enstamobile) at a rate of 10 Hz. When focused in air with a lens of 5-m focal length, they create a bundle of ionized filaments merging into a cylindrical plasma of average density n e ~ 10 17 cm -3 , roughly one millimeter in diameter and 1-2 meter in length (see Fig. 2) [16]. A microwave field is delivered to a small section of the plasma using a well-known source called surfatron [17]. The surfatron consists of a commercial microwave generator (magnetron) emitting continuously at 2.45 GHz with a maximum power of 500 Watt (GMP 12 KED by SAIREM S.A.). The microwave is injected as a TM mode via a capacitive coupler in a 40-mm long, annular brass cavity surrounding the filament bundle, itself contained in a fused silica tube. Transfer of microwave to the plasma region occurs via a leak in the annular cavity, an oscillating longitudinal E-field being imposed across the launching gap. The diameter of the surfatron hole is 6 mm. In our experiment, the microwave radiation was produced continuously while the laser operated in single pulse mode. No provision was made to synchronize the time of arrival of the laser pulse with a cycle of the microwave field. To measure the transfer of microwave energy to the filament plasma, we detected the induced longitudinal electronic currents in the plasma by using a fast I-dot probe placed outside the surfatron cavity, at a distance D. The I-dot probe allows the contactless measurement of longitudinal currents in the plasma [18]. It was connected with a shielded 50-ohm cable to an oscilloscope of 10-GHz bandwidth. We checked that no signal was detected by the I-dot probe in the absence of laser filament, even at maximum microwave power. Fig. 3 Longitudinal current signal measured by the I-dot probe as a function of microwave power in the presence of filament.

The probe is placed away from the surfatron exit (D = 6.3 cm).

In Fig. 3, we show the signals registered in the presence of the filaments by the I-dot probe for different microwave powers at a fixed position D = 6.3 cm. Since there is no synchronization between the laser pulse and the surfatron, different laser shots arrived randomly in the cavity at different phases of the microwave field. We have selected from a large number of shots the highest signals recorded for a given microwave value. They correspond to the arrival of the pulse at the positive or negative peak of a microwave cycle. The weak signal recorded in the absence of microwave (black curve) corresponds to the longitudinal oscillation that is imparted to the plasma during its formation, due to the ponderomotive force of the laser pulse [2].

We have measured the maximum signal detected by the I-dot probe as a function of distance D. The maximum signal divided by the maximum signal in the absence of microwave field is plotted in Fig. 4 for a microwave injected power of P = 100 W. It shows an amplification of the longitudinal current over a length longer than 15 cm beyond the exit of the surfatron. To interpret the flow of excitation along the filament, we consider the propagation of a surface guided TM wave. The filament is modeled as a homogeneous column of plasma of radius r = a and complex susceptibility p. The index of refraction of air and plasma are noted n 0 (r > a) and n p (r < a) respectively. For cold plasmas the index of refraction n p can be expressed by

22 2 / 1, 1 p pp e n j        (1)
where  p is the plasma circular frequency and  e is the frequency collision between electrons and neutral molecules.

In the formulas e is the charge of the electron, m e its mass, and n e is the electronic density of the plasma. The field of a forward propagating mode can be written in the form:

( , ) exp( ) exp( ) ( , ) exp( ) exp( )

eff eff E r j t jk n z H r j t jk n z          (3)
Here  is the circular frequency of the microwave emitted by the surfatron, k is the vacuum wave number, n eff is the complex effective index and (r, θ, z) are the coordinates in cylindrical coordinates, where the Z axis is defined by the filament axis. Due to axial symmetry we restrict the solution to a TM 0 mode. By adopting non dimensional notations we let:

2 2 2 2 0 . eff eff p W k a n n and U k a n n     (4)
With these notations the microwave field can be expressed for r > a (in air)

0 1 0 0 0 ( , ) 0 ( , ) r eff z ka Ur E r j A n I U a E Ur E r A I a                      , (5) 2 1 0
( , )

pp rz p p HH ka Ur H r j A n I Ua             , (6)
and for r < a (inside the plasma)

0 1 0 0 0 0 () ( , ) () 0 () ( , ) () p r eff p p z IU ka Wr E r j A n K W K W a E IU Wr E r A K K W a                      , ( 7
) 2 0 1 0 0 () ( , ) () p p r z p g H H IU ka Wr H r j A n K W K W a              , (8)
where I and K are the classical notations for Bessel functions and A is a constant in V/m. Note that for simplicity in the expressions above µ 0 does not appear so that H and E are both expressed in electric unities. The field continuity at the plasma boundary leads to the dispersion equation that reads [18]:

0 2 11 2 00 ( ) ( ) 0. ( ) ( ) p n I U K W U I U n W K W  (9)
The effective index n eff is obtained by solving numerically the dispersion equation in the complex plane. The imaginary part of n eff gives the attenuation of the field during propagation. The pointing vector S z takes the form

** 0 1 , 4 z r r S c E H E H   
  (10) and is used to calculate the radiated power and longitudinal electric field in the plasma with . The first case: a = 100 m models a single filament. The second case a = 450 m models a bundle of focused filaments and corresponds to the realistic case of our experimental conditions. The amplitude of the propagating microwave field calculated for a mm-diameter plasma is similar to the experimental curve of Fig. 4 and supports the hypothesis of a surface wave coupling. Assuming a perfect coupling for a = 450 µm, we find n eff = 1.015-0.3203i and the magnitude of the longitudinal component of the electric field of the TM 0 mode at the interface is found to be E z = 1.2 kV/cm for 100 W of microwave power. What is the fate of the amplified longitudinal current? As discussed now, it leads to the amplified emission of a characteristic short conical THz pulse by the plasma. This pulsed emission of a few ps duration is due to the heavily damped longitudinal free plasma oscillations [18]. It presents a peak at 1 THz with a 3-THz bandwidth, ranging from low frequencies up to several THz. As shown in ref. [20], a DC electric field leads to an amplification of this broadband THz emission. Since in the present case the duration of the laser pulse (~100 fs) and of the damped longitudinal oscillations of the plasma (1-2 ps) are very short as compared to a microwave cycle (~1 ns), the plasma can be viewed as under the influence of a constant electric field. We have measured the amplification of the THz emission through three spectral windows, around 3 GHz, 12 GHz and 91 GHz.

The I-dot detects the low frequency components of this amplified short THz pulse between 1 and 5 GHz. Fig. 6 shows the normalized spectra calculated by taking the Fourier transform of the signals of Fig. 3. Also shown is the injected microwave spectrum. To measure the microwave line, a small metallic wire was introduced in the surfatron along z axis in the absence of filament. In this case, the microwave signal recorded is quasi monochromatic and centered at 2.45 GHz (see violet curve). In the presence of filaments, the microwave signal covers the spectral bandwidth of the I-dot probe, as expected.

We have also detected an amplification of the transverse THz emission from the filament at 12 GHz with a homodyne detector and at 91 GHz with a selective heterodyne detector. By comparing the amplification of the signal obtained with the surfatron to the amplification of THz radiation in the presence of a longitudinal DC field [20], one can estimate the magnitude of the longitudinal component of the electric microwave field acting on the plasma to be ~1 kV/cm for an incident microwave power of 100 W, in good agreement with the calculated value of the longitudinal component of the electric field at the interface between plasma and air. It is instructive to evaluate the flow of excitation with other injected RF frequencies. In table I we show the calculated propagation distance for different frequencies and plasma diameters. As can be seen, the propagation distance is between one half and one third of the excitation wavelength, thus the excitation transfer distance covers several meters with MHz frequencies. With a sufficiently strong field amplitude E ~ 15 kV/cm, corresponding to a RF wave of a few tens of kW power, an avalanche regime would be reached, leading to a self-sustained virtual antenna [21]. An avalanche induced breakdown would lead to a high density hot plasma that would behave like a metallic antenna, as demonstrated in ref. [6]. In conclusion, we have shown evidence of the coupling of microwave radiation to a multi-filament plasma column and its propagation out of the excited region of the plasma in the form of a TM wave. This leads to the amplification of longitudinal oscillations of the plasma over a distance exceeding 10 cm. We predict that this transfer could reach several meters if lower frequency is injected in a larger plasma radius.

Fig. 1

 1 Fig. 1 (a) General set-up. (b) Photograph of the experiment.

Fig. 2

 2 Fig. 2 Side view (yz, x = 0 plane) of the intensity map of a multifilament bundle. The laser is impinging from the left. Dark colors (color scale > 100) correspond to ionized regions.

Fig. 4

 4 Fig. 4 Maximum signal measured by the I-dot probe divided by the maximum signal in the absence of microwave field as a function of the position along the filament measured for P = 100 Watt.

Fig. 5

 5 Fig. 5 shows the calculated attenuation of the electric field component for two typical radii of the plasma. Here

Fig. 5 1 .

 51 Fig. 5 Normalized amplitude of the electric field of the TM surface wave as a function of propagation length calculated with: n e = 10 17 cm -3 ; n 0 = 1;  e = 15x10 12 s -1 . In black: a = 100 m (attenuation length of 7 mm). In red: a = 450 m (attenuation length of 6 cm). The dashed horizontal line indicates the attenuation length defined at 1/e.

Fig. 6

 6 Fig. 6 Normalized spectral amplitude of longitudinal current as a function of frequency for different injected microwave powers.

TABLE I .

 I Calculated attenuation length of the RF longitudinal electric field for various frequencies and plasma diameters.

	Frequency ν	Wavelength 	Plasma radius a	Attenuation length L
	(GHz)	(m)	(mm)	(m)
	2.45	0.122	0.45	0.06
	0.245	1.22	1.3	0.66
	0.0245	12.2	3	4

Acknowledgments

This project has been supported by the French Direction Generale de l'Armement (Grant N° 2013.95.0901 and N°2009.34.0031). The authors acknowledge technical support by Yves-Bernard André, Jérôme Carbonnel, Guillaume Point, Bruno Honnorat and fruitful discussions with Sylvain Flajolet and Besma Larbi.


Jean-Guillaume Dumas
email: jean-guillaume.dumas@imag.fr

Philippe Elbaz-Vincent

Pascal Giorgi
email: pascal.giorgi@univ-perp.fr

Anna Urba
email: anna.urbanska@imag.fr

Parallel computation of the rank of large sparse matrices from algebraic K-theory

This paper deals with the computation of the rank and some integer Smith forms of a series of sparse matrices arising in algebraic K-theory. The number of non zero entries in the considered matrices ranges from 8 to 37 millions. The largest rank computation took more than 35 days on 50 processors. We report on the actual algorithms we used to build the matrices, their link to the motivic cohomology and the linear algebra and parallelizations required to perform such huge computations. In particular, these results are part of the first computation of the cohomology of the linear group GL7(Z).

INTRODUCTION 1.1 Motivation from K-theory

Numerous problems in modern number theory could be solved or at least better understood, if we have a good knowledge on the algebraic K-theory (or motivic cohomology) of integers of number fields or the cohomology of arithmetic groups (i.e. subgroups of finite index of GLN (Z)). As a short list, we could mention:

• modular forms and special values of L functions,

• Iwasawa theory and understanding of the "cyclotomy",

• Galois representations (or automorphic representations).

Let us explain first what is algebraic K-theory: to a commutative ring R we can associate (functorially) an infinite family of abelian groups Kn(R) which encodes a huge amount of information on its arithmetic, geometric and algebraic structures. These groups extend some classical notions and give higher dimensional analogues of some well known results. For instance if R is a ring of integers or a polynomial ring over a field, we have • K0(R) is the classical Grothendieck group (classifying finitely generated R-modules),

• K1(R) is the group of invertibles of R,

• K2(R) classifies the universal extensions of SL(R) and is related to the Brauer group in the case of a field.

We can give a general abstract definition of Kn for n > 0,

Kn(R) = πn(BGL(R) + ) ,
where BGL(R) + is the Quillen +-construction applied to the classifying space BGL(R) and πn denotes the nth homotopy group [START_REF] Rosenberg | Algebraic K-Theory and its applications[END_REF]. We also can see the K-groups as a way to understand GL(R). For instance, K1(R) is isomorphic to GL(R) modulo elementary relations. For a more detailed background on K-theory and its applications see [START_REF] Rosenberg | Algebraic K-Theory and its applications[END_REF].

Fact: these groups are hard to compute.

The Vandiver conjecture as an illustration

The Vandiver conjecture plays an important role in the understanding of the "cyclotomy" in number theory. Its statement is as follows; Conjecture of Vandiver: Let p be an odd prime, ζp = e 2πi/p , C the p-Sylow subgroup of the class group of Q(ζp)

and C + the subgroup fixed by the complex conjugation. The Vandiver conjecture is the statement that C + = 0.

To illustrate its interplay with K-theory, we have Fact [START_REF] Kurihara | Some remarks on conjectures about cyclotomic fields and K-groups of Z[END_REF] [START_REF] Kurihara | Some remarks on conjectures about cyclotomic fields and K-groups of Z[END_REF]: If K4n(Z) = 0 for all n > 0, then the conjecture of Vandiver is true. Some partial results on this conjecture and its connection with the cohomology of SLN (Z) are given in [START_REF] Soulé | Perfects forms and the vandiver conjecture[END_REF].

General problem: Find explicit methods for computing (co)homologies of arithmetic groups and the K-theory of number fields (or their ring of integers).

Our first task will be the computation of the (co)homologies of linear groups (mainly GLN (Z) and SLN (Z)). We can show [START_REF] Elbaz-Vincent | Quelques calculs de la cohomologie de GLN (Z) et de la k-théorie de Z[END_REF][START_REF] Elbaz | Perfects lattices, homology of modular groups and algebraic k-theory[END_REF][START_REF] Elbaz-Vincent | Perfect forms, cohomology of modular groups and k-theory of integers[END_REF] that the computation of those groups is the key point for computing K-groups with our method.

We will begin to recall some facts from topology.

Topological Excursion: Cellular complexes, or how to simply describe the "combinatorial" structure of a topological space

The notion of cellular complex is a generalization of a graph to several dimensions. We call n-cell a topological space homeomorphic to the open unit ball of R n and such that its closure is homeomorphic to the closed unit ball. A cellular complex (or cell complex or also cellular decomposition or cellular space) is a family of sets X n (with n ∈ N), such that each X n is a collection (eventually infinite) of n-cells.

Usually we work with cell complexes with a finite number of cells. A classical result [START_REF] Spanier | Algebraic Topology[END_REF] shows that any (reasonable) topological space can be approximated by such cell complexes.

A cellular decomposition of the cube

Computation of the homology of a cell complex

To a cell complex, we can associate a family Cn (n ∈ N) of free Z-modules and a family dn : Cn → Cn-1 of linear maps.

The module Cn is the free module with basis the n-cells (modulo a choice of orientation).

If the complex is finite, then all the modules are of finite rank and we will denote by {b n λ } λ∈Λn a basis of Cn, Λn being an index set for the n-cells. Then the map dn is defined by

dn(b n λ ) = µ [b n λ : b n-1 µ ]b n-1 µ ,
and the integer number

[b n λ : b n-1 µ
] is called the incidence number of the cell e n-1 µ inside the cell e n λ . The relation dn • dn+1 = 0 (i.e., formally d 2 = 0) should hold for any n. If the complex is regular (i.e., always at most one cell of dimension n + 1 between two cells of dimension n), then we can build the incidences inductively starting from the 0-cells up to the maximal cells using the d 2 = 0 condition and moreover the incidence numbers will be 0, or ±1. The nth homology group of the complex is defined as the quotient of Ker(dn) by Im(dn+1). This construction is functorial (in the category of cell complexes). As a consequence, we can determine the homology groups effectively by computing the Smith form of the integral matrices of the dn (relatively to the fixed basis). Notice that the Smith form gives both the rank of the free part and the explicit description of the torsion part. In case the computation of the torsion is unnecessary (or to difficult to achieve), we can tensorize by -⊗ Z Q and the homology groups become Q-vector spaces with their dimensions given by the ranks of the matrices of the differentials.

In general, the matrices of the differentials can be very large even for a relatively simple cell complex. However, they are also very sparse and we may look for some other favorable properties which would enable the computation despite the size of the problem.

How to use such settings for the computation of linear groups ?

If G is a group acting on a cell space X (i.e., G sends n-cells to n-cells), then, under some technical assumptions on X and on the action, we can show [START_REF] Brown | Cohomology of groups[END_REF] that roughly computing the homology of G (as group homology) is the same as computing the homology of the cell space X/G. Hence, if X/G can be calculated effectively, we can compute explicitly its homology, and from this the homology of G (similarly for the cohomology). Notice that in general the space X/G will not be regular anymore. The main difficulty is to find a cell space X such that X/G will be effective. We will discuss in section 2.1 how we can construct such cellular space for linear groups.

Parallelism motivations

This first idea to deal with very large sparse matrices is to use them as blackboxes, i.e. only using the matrix-vector product. This will let the matrix remain sparse all along the algorithm where Gaussian elimination for instance would fill it up. To compute the rank, the fastest black box algorithm is Wiedemann's as shown in [START_REF] Dumas | Computing the rank of sparse matrices over finite fields[END_REF]. This algorithm computes a sequence of scalars of the form u t A i v (u and v are vectors) with i matrix vector products and dot products. It has been shown to successfully deal with large sparse matrix problems from homology, see e.g. [START_REF] Dumas | On efficient sparse integer matrix Smith normal form computations[END_REF]. Nevertheless, when matrices are very large (e.g millions on non-zero entries) computations would require months or years. This is due to the low practical efficiency of the computation of a sparse matrix-vector product. For instance, in our case of homology computation, one would need 300 days of CPU to compute the sequence involving matrix of GL7(Z) with n = 19 (GL7d19 matrix). To achieve computations of many homologies in a realistic time we then need to parallelize the computation of the sequence. Then the algorithm candidate is the block Wiedemann method, which computes a sequence X T A i Y where X and Y are blocks of vectors. This step can be easily parallelized by distributing vectors of block Y to several processors. We thus have several objectives with regard to the parallelism in this paper:

• We want to solve large problems coming from homology computation.

• We want to experimentally validate our parallel implementation of the block Wiedemann rank algorithm.

• We want to show the parallel scaling of block Wiedemann approaches.

Summary of the paper

In section 2, the algorithms and optimizations we used to generate the matrices and compute with them are discussed. Then, section 3 shows our experimental results on these large sparse matrices coming from homology.

ALGORITHMS AND OPTIMIZATIONS

As seen in section 1.1.3, we can effectively compute the homology of a cellular space and of a group which acts "nicely" on a cellular space. The main difficulty remains to find such explicit cellular space. In section 2.1 we present the process of matrix generation and optimization. Then in the following sections we give a description of the algorithms used for the computation of the rank and the Smith form of those matrices.

Matrices generation

In the case of subgroups of GLN (Z), we have an "obvious" action on Z N . We can then capture the topology by regarding Z N not as a free Z-module but as a lattice (or equivalently as quadratic forms), and see if this leads to some interesting topological construction. We will describe this approach below and the results that we can get.

Voronoï's reduction theory

Let N 2 be an integer. We let CN be the set of positive definite real quadratic forms in N variables. Given h ∈ CN , let m(h) be the finite set of minimal vectors of h, i.e. vectors

v ∈ Z N , v = 0, such that h(v) is minimal. A form h is called perfect when m(h) determines h up to scalar: if h ′ ∈ CN is such that m(h ′ ) = m(h), then h ′ is proportional to h.
Example 2.1. The form q(x, y) = x 2 +y 2 has minimum 1 and minimal vectors ±(1, 0) and ±(0, 1). Nevertheless this form is not perfect, because there is an infinite number of definite positive quadratic forms having these minimal vectors. On the other hand, the form q(x, y) = x 2 + xy + y 2 has also minimum 1 and has exactly 3 minimal vectors (up to sign), the one above and ±(1, -1). This form is perfect, the associated lattice is the "honeycomb lattice" (with optimal spheres packing in the plane), it is the only one.

h • γ = γ t h γ , γ ∈ Γ , h ∈ C * N ,
where h is viewed as a symmetric matrix and γ t is the transposed of the matrix γ. Voronoï proved that there are only finitely many perfect forms modulo the action of Γ and multiplication by positive real numbers ( [START_REF] Voronoi | Nouvelles applications des paramètres continus à la théorie des formes quadratiques i[END_REF], Th. p. 110).

Given v ∈ Z N -{0} we let v ∈ C * N be the form defined by

v(x) = (v | x) 2 , x ∈ R N ,
where (v | x) is the scalar product of v and x. The convex hull of a finite subset B ⊂ Z N -{0} is the subset of X * N image by π of the elements j λj vj , vj ∈ B, λj 0. For any perfect form h, we let σ(h) ⊂ X * N be the convex hull of the set m(h) of its minimal vectors. Voronoï proved in [32, § 8-15], that the cells σ(h) and their intersections, as h runs over all perfect forms, define a cell decomposition of X * N , which is invariant by the action of Γ. We endow X * N with the corresponding CW -topology. If τ is a closed cell in X * N and h a perfect form with τ ⊂ σ(h), we let m(τ ) be the set of vectors v in m(h) such that v lies in τ . Any closed cell τ is the convex hull of m(τ ) and m(τ

) ∩ m(τ ′ ) = m(τ ∩ τ ′ ).

Voronoï's complex

Let d(N ) = N (N + 1)/2 -1 be the dimension of X * N and n d(N ) a natural integer. We denote by Σ * n a set of representatives, modulo the action of Γ, of those cells of dimension n in X * N which meet XN , and by Σn ⊂ Σ * n the cells σ such that the stabilizer Γσ of σ in Γ preserves its orientation. Let Vn be the free abelian group generated by Σn. We define as follows a map dn : Vn → Vn-1 .

For each closed cell σ in X * N we fix an orientation of σ, i.e. an orientation of the real vector space R(σ) of symmetric matrices spanned by the forms v, v ∈ m(σ). Let σ ∈ Σn and let τ ′ be a face of σ. Given a positive basis B ′ of R(τ ′ ) we get a basis B of R(σ) by adding after B ′ a vector v, v ∈ m(σ)m(τ ′ ). We let ε(τ ′ , σ) = ±1 be the sign of the orientation of B in the oriented vector space R(σ) (this sign does not depend on the choice of v).

Next, let τ ∈ Σn-1 be the cell equivalent to τ ′ = τ • γ. We define η(τ, τ ′ ) = 1 (resp. η(τ, τ ′ ) = -1) when γ is com- patible (resp. incompatible) with the chosen orientations of R(τ ) and R(τ ′ ).
Finally we define

dn(σ) = τ ∈Σ n-1 τ ′ η(τ, τ ′ ) ε(τ ′ , σ) τ , (1)
where τ ′ runs through the set of faces of σ which are equivalent to τ . It is shown in [START_REF] Elbaz-Vincent | Perfect forms, cohomology of modular groups and k-theory of integers[END_REF], that up to p-torsions with p N + 1, the homology of this complex computes the cohomology of G.

For N = 5, 6, 7 we get the following results for Σn. The previous result gives the precise size of the matrices involved in the computation of the homology.

The main challenge was then the computation of the ranks of the matrices of the differential (this gives the free part of the homology) and the computation of the Smith forms (which gives the relevant arithmetical information of the homology), in particular for N = 7, knowing that such matrices are particularly sparse. We can emphasize the fact that what we want to detect is the "high torsion" in the homology (i.e. the prime divisors > 7 of the Smith invariants).

In the following paragraphs we will discuss the different methods chosen for the computations and to take up the challenge 1 .

Coppersmith Block Wiedemann

One successful approach to deal with linear algebra computations on large sparse matrices is to rely on Lanczos/Krylov black-box methods. In particular, block versions of Wiedemann method [START_REF] Wiedemann | Solving sparse linear equations over finite fields[END_REF] are well suited for parallel computation. This technique has been first proposed by Coppersmith in [START_REF] Coppersmith | Solving homogeneous linear equations over GF[2] via block Wiedemann algorithm[END_REF] for computation over GF (2) and then analyzed and proved by Kaltofen [START_REF] Kaltofen | Analysis of Coppersmith's block Wiedemann algorithm for the parallel solution of sparse linear systems[END_REF] and Villard [START_REF] Villard | Further analysis of Coppersmith's block Wiedemann algorithm for the solution of sparse linear systems[END_REF][START_REF] Villard | A study of Coppersmith's block Wiedemann algorithm using matrix polynomials[END_REF]. The idea is for a matrix A ∈ F n×n to compute the minimal generating matrix polynomial of the matrix sequence {XA i Y } ∞ i=0 ∈ F s×s , 1 All the matrices are available on line in the "Sparse Integer Matrix Collection" (ljk.imag.fr/membres/ Jean-Guillaume.Dumas/simc.html) where X, Y are blocks of s vectors (instead of vectors in the original Wiedemann's algorithm). Therefore, an intuitive parallelization of this method is to distribute the vectors of the block X, Y to several processors. Thus, the computation of the sequence, which is the major performance bottleneck of this method, can be done in parallel and then allow for better performance.

Lots of implementations and practical experimentations has been developed on parallel block Wiedemann. For instance, in 1996, Kaltofen and Lobo [START_REF] Kaltofen | Distributed matrix-free solution of large sparse linear systems over finite fields[END_REF] have proposed a coarse grain implementation to solve homogeneous linear equations over GF (2). They have thus been capable to solve a system of 252 252 linear equations with about 11.04 million non-zero entries, in about 26.5 hours using 4 processors of an SP-2 multiprocessor.

Lately, in 2001, Thomé in [START_REF] Thomé | Fast computation of linear generators for matrix sequences and application to the block Wiedemann algorithm[END_REF] improved Coppersmith's algorithm by introducing matrix half-gcd's computation, and its implementation [START_REF] Thomé | Subquadratic computation of vector generating polynomials and improvement of the block Wiedemann algorithm[END_REF] was able to outperform Kaltofen-Lobo's software. One may remark that introduction of matrix gcd was first suggested by Villard in [START_REF] Villard | Further analysis of Coppersmith's block Wiedemann algorithm for the solution of sparse linear systems[END_REF] who relied on the work of Beckermann and Labahn [START_REF] Beckermann | A uniform approach for the fast computation of matrix-type Padé approximants[END_REF] on power Hermite Padé approximation. Finally, Giorgi, Jeannerod and Villard have generalized in [START_REF] Giorgi | On the complexity of polynomial matrix computations[END_REF] block Wiedemann algorithms by introducing σ-basis computation and then reducing the complexity to polynomial matrix multiplication. A sequential implementation of this algorithm is now available in the LinBox library (www.linalg.org).

Block symmetry

In order to reduce the number of dot products, we used a symmetric projection. In other words, we set X = Y T in the XA i Y sequence. Indeed, in this case the probability of success is reduced but the obtained block is symmetric as soon as A is symmetric. This is always the case when the preconditioners of [START_REF] Dumas | Computing the rank of sparse matrices over finite fields[END_REF] are used (they are of the form A T A). This reduces the dot product part of the computation of the sequence by a factor of two. For instance column i can be deduced from its top i elements and row i. This induces some load balancing issues when one process owns the computation of one column. Note also that we use BLAS level-2 for the computation of this dot products. In other words we perform them by blocks.

σ-basis computations

In order to efficiently compute σ-basis we rely on algorithm PM-Basis of [START_REF] Giorgi | On the complexity of polynomial matrix computations[END_REF] which reduces this computation to polynomial matrix multiplication. One can multiply two polynomial matrices A, B ∈ F n×n [x] of degree d in O(n 3 d + n 2 d log d) finite field operations if d-th primitive roots of unity are available in F. Consequently, we decided for our computations to define F as a prime field with primes of the form c × 2 k + 1 such that c × 2 k ≡ 0 mod d. These primes are commonly called FFT primes since they allow the use of FFT evaluation/interpolation for polynomials of degree d. We refer the reader to [START_REF] Cantor | On fast multiplication of polynomials over arbitrary algebras[END_REF][START_REF] Bostan | Polynomial evaluation and interpolation on special sets of points[END_REF] and references therein for further informations on fast polynomial matrix arithmetic.

When finite fields not having d-th primitive roots of unity are used, polynomial matrix multiplication is still be done efficiently by using Chinese Remainder Theorem with few FFT primes. Let be F a prime field of cardinality p, then the multiplication of A, B ∈ F n×n [x] of degree d can be efficiently done by using CRT with FFT primes pi satisfying pi > d × n × p 2 . This is equivalent to perform the multiplication over the integers and then reduce the result in F. The overall performance of the multiplication, and then of the σ-basis , is dependent on the numbers of FFT primes needed.

Rank

Our main interest in the block Wiedemann approach is to compute the rank of large sparse matrices given by the Homology group determination problem explained in section 1.1. Hence, we rely on Kaltofen-Saunder's rank algorithm [START_REF] Kaltofen | On Wiedemann's method of solving sparse linear systems[END_REF] and its block version [START_REF] Turner | A block wiedemann rank algorithm[END_REF] to achieve efficient parallel computation.

The Kaltofen-Saunders approach is based on the fact that if à is a good preconditioned matrix of A then its rank is equal to the degree of its minimal polynomial minus its valuation (or co-degree) [START_REF] Kaltofen | On Wiedemann's method of solving sparse linear systems[END_REF]. Thus, by using well chosen preconditioners and Wiedemann algorithm one can easily compute the rank of a sparse matrix over a finite field. The block version of this method is presented e.g. in [29, §4]. We recall now the basic outline of this algorithm.

Block Wiedemann Rank Algorithm

: let A ∈ F n×n ,
1 form à from A with good preconditioners (e.g. those of [START_REF] Dumas | Computing the rank of sparse matrices over finite fields[END_REF]). Note that if the minimal matrix of step 3 is in Popov form (e.g. computed using the σ-basis of [START_REF] Giorgi | On the complexity of polynomial matrix computations[END_REF]), then the degree of det(F à Y ) is simply the sum of the row degrees of the matrix F à Y . Then the co-degree is zero if the determinant of the constant term of F à Y , seen as a matrix polynomial, is nonzero. In the latter case the computation of the determinant of the whole polynomial matrix can be avoided.

When this fails, this determinant is computed by a massively parallel evaluation/interpolation. It could be interesting, though, to interpolate only the lower coefficients of this polynomial incrementally. This was not required for the matrices we considered and we therefore did not investigate more on these speed improvements.

Note that to probabilistically compute the rank over the integers, it is sufficient to choose several primes at random and take the largest obtained value, see e.g. [START_REF] Dumas | On efficient sparse integer matrix Smith normal form computations[END_REF] for more details. Moreover, one can choose the primes at random among primes not dividing the determinant (and thus preserving the rank). In order to ensure this property it it sufficient to select primes not dividing the valence or last invariant factor computed by one of the methods of next section.

Smith form

The computation of the Smith form for the matrices of GL7(Z) turned out to be a very challenging problem.

Smith for via the Valence

Prior experience with sparse homology matrices led us to try the SmithViaValence algorithm of [START_REF] Dumas | On efficient sparse integer matrix Smith normal form computations[END_REF]. The idea is to compute the minimal valence (the coefficient of the smallest non zero monomial of the minimal polynomial) of the product A T A to determine the primes p which divide the invariant factors of the Smith form of A. When the primes have been found, one can compute the local Smith forms of A at each p separately and return the resulting Smith form S as the product of the local Smith forms Sp over all p. Local Smith form computation can be done by a repeated Gauss elimination modulo p e where the exponent e is adjusted automatically during the course of the algorithm. This algorithm works very efficiently for sparse matrices provided that the minimal polynomial of the product AA T has a small degree. Unfortunately, the latter condition does not hold in the case of GL7(Z) matrices. Moreover, some early experiments with small matrices showed that much more primes occur in the computed valence than in the Smith form of the original matrix.

Saunders and Wan's adaptive algorithm

Thus, we decided to apply the adaptive algorithm of Saunders and Wan [START_REF] Saunders | Smith normal form of dense integer matrices, fast algorithms into practice[END_REF] which is a modified version of Eberly-Giesbrecht-Villard algorithm [START_REF] Eberly | On computing the determinant and Smith form of an integer matrix[END_REF]. In [START_REF] Eberly | On computing the determinant and Smith form of an integer matrix[END_REF] the authors proposed a procedure OneInvariantFactor(i, A) (OIF) which computes the ith invariant factor of a n × n matrix A. Then the binary search for distinct invariant factors allows them to find the Smith form of A. OIF reduces the ith factor computation to the computation of the last (nth) invariant factor (LIF) of a preconditioned matrix A + UiVi, where Ui, V T i are random n × (ni) matrices. In [START_REF] Saunders | Smith normal form of dense integer matrices, fast algorithms into practice[END_REF] the method was extended to handle the rectangular case of m × n matrix. It is done by computing the last (ith) invariant factor of a preconditioned matrix LiARi where Li is a m × i and Ri is a i × n matrix.

The procedure OIF is of Monte Carlo probabilistic type where the probability of correctness is controlled by repeating the choice of preconditioners. Assuming the correctness of LIF computation, it gives a multiple of the ith invariant factor. In practice, LIF is also of randomized Monte Carlo type. The idea is to get a divisor of the last invariant factor by solving a linear equation M x = b with random righthand side. After several solvings we get the last invariant factor with large probability. Thus, the overall situation is more complex and we cannot exclude the possibility that some primes are omitted or unnecessary in the output of OIF. However, the probability that a prime is omitted or is unnecessary in this output can be controlled for each prime separately and is smaller for bigger primes.

Therefore in [START_REF] Saunders | Smith normal form of dense integer matrices, fast algorithms into practice[END_REF] the authors introduce a notion of smooth and rough parts of the Smith form. The idea is to compute the local Smith form for smaller primes by for example the SmithViaValence or OIF algorithm and to recover only large primes with the invariant factor search of [START_REF] Eberly | On computing the determinant and Smith form of an integer matrix[END_REF]. When we consider large primes, a sufficient probability of correctness can be obtained by a smaller number of repetitions.

More adaptiveness

As we did not want to compute the valence, we introduced some minor changes to the algorithm, which at the end is as follows:

1. r = rank(A)

For primes 1 < p < 100 compute the local Smith form

Sp of A;

3. Compute sr(A) by OneInvariantFactor algorithm;

4. P = all primes p > 100 which divide sr(A);

5. If P = ∅ return S = Π p Sp;
One advantage of this method is that we get the information on the smooth form of the matrix very quickly. Moreover, the OIF computation acts as a certification phase which allows us to prove that no other primes are present with a sufficiently large probability. This probability is explicit in the following theorem:

Theorem 2.3. The probability that there exists a prime p > P that divides the ith invariant factor but does not divide the output of OIF which uses M random preconditioners Li, Ri and N random vectors b, with b ∈ {0, 1, . . . β -1} and β > si(A), in the LIF procedure is bounded by

M ∞ p>P 2 p N .
Proof. As we take the gcd of the result with different preconditioners Li and Ri, is suffices that the LIF computation fails in one case to spoil the computation. We are free to choose a large bound for b such that b > si(A) without increasing the complexity of LIF computation. Then the probability that a prime p < b is omitted in LIF is less than or equal to 1 β ⌈ β p ⌉ < 2 p , see [START_REF] Abbott | Fast deterministic computation of determinants of dense matrices[END_REF]. Finally, we bound the probability that any prime p > P, p | si(A) is missing by taking a sum over all primes.

The choice M = N = 2 suffice to obtain a small probability 0, 015 of omitting an important prime, and at the same time to exclude all primes that are not in the ith invariant factor. In our experiments, there was no need to perform the computation for any additional prime p > 100 as all the primes were excluded by the OIF computation. This is one of the most important advantages over the valence computation.

The algorithm [START_REF] Saunders | Smith normal form of dense integer matrices, fast algorithms into practice[END_REF] was stated in the case of dense matrices. We slightly modified it in order to exploit the sparse structure of the matrix. In particular, we used the sparse local Smith form computation of [9, Algorithm LRE] but stick to the dense Dixon solver [START_REF] Dixon | Exact solution of linear equations using p-adic expansions[END_REF] as long as the memory was sufficient. Any other solver, including the new sparse solver of [START_REF] Eberly | Solving sparse rational linear systems[END_REF] could potentially be used for larger matrices.

The limits of this method are imposed by the available memory. For example, it was possible to use the dense Dixon solver only for the six smallest matrices. Furthermore, sparse elimination reached its limits for matrices of size greater than 171375 × 47271 and 21074 × 105054 when the filling of the matrices started to be impossible to handle. For the 460261×171375 matrix GL7d15 and 105054×349443 matrix GL7d23, specialized space-efficient elimination procedures mod 2, 3 and 5 allowed us to compute the rank mod 2, 3 and 5 respectively.

Chain Reductions

The encountered problems have shown a need for a more elaborated reduction algorithm. We focused our attention on the algebraic reduction algorithm for chain complexes of [START_REF] Kaczyński | Homology computation by reduction of chain complexes[END_REF]. We implemented a simplified version of the algorithm in the language of matrices using the LinBox library. The heuristic behind this algorithm is that Gaussian elimination can propagate from one matrix of a chain complex to the next thanks to the exactness of the differential map (i.e. d 2 = 0 condition).

The motivations come from the geometric properties of homologies. By a free face we refer to a (k -1) cell a which is in the differential of only one k-dimensional cell b. By removing the pair (a, b) we obtain a retract of the initial cell complex (viewed as a geometrical object), see Figure 3. The process can be repeated. From the homology theory we know that the groups of homology are the same for the set and its retract. The removal of pairs leads to a reduction of the basis of the cell complex. We refer to [START_REF] Kaczyński | Computational Homology[END_REF]Ch.4] for a full description of the procedure.

If the differential map is represented by matrices whose rows represent the cells of dimension k -1 and columns -dimension k, the removal of a pair (a, b) can be interpreted as the

dk v = d k (v) -λ -1 [v : a]d k (b). (2)
Thus, in the basic case of free face removal no modification is needed. In the case of matrices, the formula describes a step of Gauss elimination where the reduced row is removed and not permuted. This proves that the Smith form (or the rank) of the initial and reduced matrix will be the same, provided we add a number of trivial invariant factors equal to the number of rows reduced to the reduced Smith form.

The important characteristic of this methods is that, thanks to the exactness of the matrix sequence, we can also remove row b and column a from the neighboring matrices. In this way, elimination in one matrix can propagate on the others.

Due to the format of data (large files with matrix entries) we decided to implement only the simplest case of 1-rows removal which led to entries removal but no modifications. We removed the empty rows/columns at the same time and performed the whole reduction phase at the moment of reading the files. This led to vast matrix reduction in the case of GL7(Z) matrices from the beginning of the sequence. The propagation of reductions unfortunately burned out near GL7d14 matrix and stopped completely on GL7d19. Applying the process for the transposed sequence did not improve the solution. Next step would be to implement the propagation of Gauss elimination steps as in Eq. ( 2). It would be interesting to examine whether the burn-out can be connected to the loss of regularity for X/G and/or a huge rectangularity of the input matrix GL7d14.

EXPERIMENTAL RESULTS

Parallel implementation

For the sake of simplicity in our experimental validation, we decided to develop our parallel implementation using shell tasks distribution on SMP architecture. Thus, a simple script code is used to distribute all different tasks over all the processors and files are used to gather up the computed results. Our parallel implementation has been done as follow:

1. The block of vectors Y and the sparse matrix A are broadcasted on every processors.

2. Each process takes one column of the blackbox A T .A and compute the corresponding column's sequence using the first ith columns of Y . Each process writes the result in a file labeled with the corresponding index of the column sequence.

3. When all previous processes have terminated, the σbasis computation is sequentially performed after loading the sequence from the generated files.

Despite the na ïve approach used for the parallelization, our implementation authorized us to perform very large computations as show in next section. However, our experiments show a need for at least a more robust parallel computation scheduler.

Rank and Smith form

All our computation have been done on a SGI Altix 3700 gathering 64 Itanium2 processors with 192Gb memory and running SuSE Linux Enterprise System 10. Further informations on this platform are available at http://www.math. uwaterloo.ca/mfcf/computing-environments/HPC/pilatus.

In Table 1 we include the information about the dimensions of the GL7(Z) matrices and their sparsity. The matrices are very sparse which is illustrated by the fact that less than 1% of the entries are non-zero except for matrices GL7d10 and GL7d11. This value drops to less than 0,2% in the case of the largest matrices. Also in Table 1 we give the results for the rank and the Smith form computations. We have obtained a full information on the rank of GL7(Z) matrices. For the computation of the Smith form, full result has been obtained in the case of matrices 10,11,12,13,25,26. For matrices 14 and 24 only the smooth part of the Smith form has been computed. For matrices 15 and 23 we have proved the existence of a non-trivial local Smith form at 2 and 3 and a triviality of the local Smith form at 5. As the result for these matrices we give the number of invariant factors divisible by 2 and 3.

In Table 2 we give the times for the Smith form algorithms used. For cases with * no data are available or relevant.

The rank computation for GL7(Z) matrices was performed modulo 65537. This FFT prime allowed us to use both BLAS routines and sigma-basis reconstruction using fast polynomial multiplication. In Table 3 we give the timings for different operations used in the computation i.e. the sparse matrix-vector product, BLAS-based matrix-vector product and, for the sake of comparison, the time of scalar dot product equivalent to the BLAS computation.

In Table 4 we give the estimation of sequential and parallel cpu time of rank computation and compare it with the real time of parallel computation. The times are estimated based on the number of iterations and the times of one step which can be computed from Table 3 (notice, that in the scalar case we use 1 dot product instead of BLAS). The real time of computation includes the time of writing and reading the data which was considerable. The difference of the real and estimated running times may also be due to the overload of the computation cluster. Moreover, long computations suffered from system crashes and/or shutdowns.

A Some restoration scripts were used to recover the data which unfortunately required re-running some part of the computation. Thus, the real time given in Table 4[Col.5] should be treated as a rough approximation.

CONCLUSION

Using the previous methods and computations, we get the following new result for the rational cohomology of GL7(Z).

Theorem 4.1. (Elbaz-Vincent/Gangl/Soulé) [START_REF] Elbaz-Vincent | Perfect forms, cohomology of modular groups and k-theory of integers[END_REF] We have H m (GL7(Z), Q) = Q if m = 0, 5, 11, 14, 15, 0 otherwise.

Clearly the simple parallelization we used together with the highly optimized routines were the key to enable these computations.

To go further and solve even larger problems, it is mandatory to improve the parallelism. On SMP we can split the matrices into blocks and perform the matrix-vector products with different threads. This, and the unbalanced load we had when we choose to assign one vector to one process, advocates for the use of more advanced scheduling. We are experimenting KAAPI2 but were not ready for the computation of GL7.

Other improvements are of algorithmic type. They include the use of the sparse projections of [START_REF] Eberly | Solving sparse rational linear systems[END_REF] for the matrix sequence. But then we loose the symmetry of the projections and therefore must pay a factor of two for the number of iterations. We could also use an early termination strategy to stop the iteration earlier, but up to now this require to loose the fast algorithm for the sigma bases. Then if a good structure for the sparsity of the matrices could be found, e.g. an adapted reordering technique, this would enable an efficient clustering and therefore faster and more scalable matrix-vector products.

In order to have the relevant part of the torsion of the integral cohomology of GL7(Z), we would need the complete description of the Smith forms of all the matrices described above. Our experiments have shown that this can be an enormously difficult task. The computation remains to be done but would have applications in number theory.

  Denote by C * N the set of non negative real quadratic forms on R N the kernel of which is spanned by a proper linear subspace of Q N , by X * N the quotient of C * N by positive real homotheties, and by π : C * N → X * N the projection. Let XN = π(CN ) and ∂X * N = X * N -XN . Let Γ be either GLN (Z) or SLN (Z). The group Γ acts on C * N and X * N on the right by the formula

Figure 2 :

 2 Figure 2: Cardinalities of Σn and Σ * n for N = 7

2

  choose random block Y ∈ F n×s and compute the matrix sequence S = {Y T Ãi Y } for i = 0 . . . 2n/s + O(1).3 compute the minimal matrix generator F Ã Y ∈ F s×s [x] of the matrix sequence S.

4

  return the rank r as r = deg(det(F Ã Y ))-codeg det(F Ã Y ).

Figure 3 :

 3 Figure 3: Retraction for a square

Table 1 :

 1 Results of the rank and Smith form computation for GL7(Z) matrix A of dimension n × m with Ω non-zero entries. For (*) the information is incomplete -only divisors of the invariant factors were determined based on the rank mod 2 and 3 computation.

		Ω	n		m	rank	ker		Smith form
	GL7d10	8		60	1	1	59			1
	GL7d11	1513	1019	60	59	960			1 (59)
	GL7d12	37519	8899	1019	960	7939		1 (958), 2 (2)
	GL7d13	356232	47271	8899	7938	39333		1 (7937), 2 (1)
	GL7d14	1831183	171375	47271	39332	132043		1 (39300),2 (29),4 (3)
	GL7d15	6080381	460261	171375	132043	28218 1 (131993), 2•? (46), 6•? (4) (*)
	GL7d16 14488881	955128	460261	328218	626910
	GL7d17 25978098 1548650	955128	626910	921740
	GL7d18 35590540 1955309 1548650	921740 1033569
	GL7d19 37322725 1911130 1955309 1033568	877562
	GL7d20 29893084 1437547 1911130	877562	559985
	GL7d21 18174775	822922 1437547	559985	262937
	GL7d22	8251000	349443	822922	262937	86506
	GL7d23	2695430	105054	349443	86505	18549	1 (86488), 2•? (12), 6•? (5) (*)
	GL7d24	593892	21074	105054	18549	2525		1 (18544),2 (4),4 (1)
	GL7d25	81671	2798	21074	2525	273		1 (2507), 2 (18)
	GL7d26	7412	305	2798	273	32		1 (258), 2 (7), 6 (7), 36 (1)
	A	ñ	m	r	Red	RAdaptive SmoothSF AdaptiveSF SFValence
	GL7d11	39	8	52	0.01s	< 10 -2 s		0.09s	0.26s	4.84s
	GL7d12	289	58	909	0.30s	0.16s		9.75s	218.68s	4.04h
	GL7d13	7938	740 7250	3.12s	159.16s		0.76h	*	2526.65h
	GL7d14 165450 35741 4279 21.62s	*		796h	*	*
	GL7d25	2797 20990	0	1.74s	*		17.67s	4.40h	52.13h
	GL7d26	302	2748	0	0.14s	*		0.29s	26.81s	274.35s

Table 2 :

 2 Times for Smith Form computation for GL7(Z) matrices. From left to right: the dimensions of the matrix after reductions, rank approximation by reductions, time of reading and reducing the matrix, time for the adaptive algorithm for a reduced matrix; times for the original matrix: smooth form computation, adaptive algorithm; valence computation in parallel -sequential time equivalent.

	A	iter [1]	time app	iter [p]	time	time app σ-basis
	GL7d11	120	0.02s	6 [30]	0.01s	< 10 -2 s	0.58s
	GL7d12	1922	7.53s	66 [30]	0.32s	0.26s	12.16s
	GL7d13	15878	880.28s	532 [30]	51.65s	29.49s 249.17s
	GL7d14	78666	5.90h	2625 [30]	0.56h	0.20h	0.45h
	GL7d15	264088	66.98h	8805 [30]	2.25h	2.23h	2.45h
	GL7d16	656438	509.80h 21884 [30]	27.29h	17.00h	6.03h
	GL7d17 1253822	113.90d 41796 [30]	14d	3.80d	0.57d
	GL7d18 1843482	256.50d 46089 [40]	28d	6.41d	1.00d
	GL7d19 2067138	321.89d 41345 [50]	35d	6.44d	1.56d
	GL7d20 1755126	212.82d 36568 [48]	10d	4.43d	1.41d
	GL7d21 1119972	75.01d 37335 [30]	5d	2.50d	0.55d
	GL7d22	525876	293.60h 17532 [30]	16.47h	9.79h	5.85h
	GL7d23	173012	21.18h	5769 [30]	1.17h	0.71h	1.09h
	GL7d24	37100	3172.79s	1239 [30] 188.78s	105.96s 666.83s
	GL7d25	5052	40.21s	171 [30]	1.56s	1.36s	41.47s
	GL7d26	548	0.40s	21 [30]	0.03s	0.02s	2.03s

Table 4 :

 4 A summary of large-scale parallel rank computation. From left to right: number of iteration in the scalar case, time estimation in this case, number of iterations on p processors computed as 2 + 2 • r/p, average (real) time of sequence computation, estimated time on p processors, the time of the σ basis computation for a sequence of length iter [p] of p × p matrices.

	A	1A T Au [s] 1U T v [s] 30u T v [s]
	GL7d11	0.0002	< 10 -4	< 10 -4
	GL7d12	0.0038	0.0001	0.0002
	GL7d13 GL7d14	0.0550 0.2677	0.0005 0.0025	0.0036 0.0190
	GL7d15	0.9048	0.0082	0.0708
	GL7d16	2.7724	0.0234	0.2641
	GL7d17 GL7d18	7.8003 11.9457	0.0485 0.0759	0.5052 0.8710
	GL7d19	13.3591	0.0948	1.0710
	GL7d20 GL7d21	10.4056 5.7461	0.0711 0.0408	0.8587 0.4604
	GL7d22	1.9919	0.0180	0.2082
	GL7d23	0.4354	0.0052	0.0459
	GL7d24 GL7d25	0.0843 0.0078	0.0012 0.0002	0.0085 0.0008
	GL7d26	0.0007	< 10 -4	< 10 -4

Table 3 :

 3 CPU timings (in sec.) for different operations used in large-scale parallel rank computation. All times in seconds. From left to right: time of a matrix-vector product, a BLAS multiplication of a vector and a 30 × min(n, m) matrix U and 30 dot products.

Kernel for Adaptive, Asynchronous Parallel and Interactive programming, kaapi.gforge.inria.fr

ACKNOWLEDGMENTS

We are grateful to Arne Storjohann and to the Computer Science Computing Facilities of the University of Waterloo for letting us fill up their SMP machine to perform our parallel computations.