Tetsuya Ito

Bert Wiest
email: bertold.wiest@univ-rennes1.fr

LAWRENCE-KRAMMER-BIGELOW REPRESENTATIONS AND DUAL GARSIDE LENGTH OF BRAIDS

Keywords: 2010 Mathematics Subject Classification. Primary 20F36, Secondary 20F10, 57M07 Lawrence-Krammer-Bigelow representation, Braid group, curve diagram, dual Garside length

We show that the span of the variable q in the Lawrence-Krammer-Bigelow representation matrix of a braid is equal to the twice of the dual Garside length of the braid, as was conjectured by Krammer. Our proof is close in spirit to Bigelow's geometric approach. The key observation is that the dual Garside length of a braid can be read off a certain labeling of its curve diagram.

Introduction

The question whether the braid group B n is linear or not was a long standing problem. At the end of the 20th century, the problem was solved affirmatively by Krammer [Kra2] and Bigelow [Big] independently. They showed that a certain linear representation of the braid group first constructed by Lawrence [Law] and now called the Lawrence-Krammer-Bigelow representation (LKB representation, for short) is faithful. Interestingly, the two proofs of the faithfulness of the LKB representations are completely different: Bigelow's proof is geometric whereas Krammer's proof is algebraic.

Our main result is (Theorem 1.1) that the span of the variable q in the image of the LKB representation of a braid β is equal to twice the dual Garside length of β, as conjectured by Krammer in [Kra1]. This is an analogue of Krammer's theorem [START_REF] Krammer | Braid groups are linear[END_REF] that the span of the variable t in the image of the LKB representation is equal to the classical Garside length. It is remarkable that, despite the analogy with Krammer's result, our proof is rather based on Bigelow's techniques.

One of our main tools is a labeling of curve diagrams of braids called the wall crossing labeling. In Theorem 3.3, we will prove that the wall crossing labeling tells us the dual Garside length of braids. On the other hand, in Lemma 4.1 we will observe that the wall crossing labeling is also related to the variable q in the noodle-fork pairing, a homological intersection pairing appearing in the LKB representation. Thus wall crossing labelings of curve diagrams serve as a bridge connecting two seemingly unrelated objects, namely the LKB representation and the dual Garside structure.

In order to state our main Theorem, we set up some notation. In this paper we denote the LKB representation by L : B n → GL n(n-1)

2

The matrix representative L(β) depends on the choice of the basis, and there are variations of the matrices of the LKB representation: all matrices given in the three papers [Big], [START_REF] Krammer | The braid group B 4 is linear[END_REF], [START_REF] Krammer | Braid groups are linear[END_REF] are different!

In this paper we use Bigelow's expression of the matrices given in [Big, Theorem 4.1], with the correction of a sign error from [START_REF] Bigelow | The Lawrence-Krammer representation, Topology and geometry of manifolds[END_REF]. This expression has geometric origin, and is essentially the same as the matrices in [START_REF] Krammer | The braid group B 4 is linear[END_REF], except that the variable t of [START_REF] Krammer | The braid group B 4 is linear[END_REF] corresponds to -t of [Big].

Let us take the basis {F i,j } 1≤i<j≤n of R n(n-1)/2 , defined by the standard forks -see Section 2 for details. Using this basis, the matrix representative of the LKB representative is given as follows.

[L(σ

i )](F j,k ) =                F j,k i ∈ {j -1, j, k -1, k} qF i,k + (q 2 -q)F i,j + (1 -q)F j,k i = j -1 F j+1,k i = j = k -1 qF j,i + (1 -q)F j,k + (q -q 2 )tF i,k i = k -1 = j F j,k+1 i = k -q 2 tF j,k i = j = k -1
The image L(β) is a n(n-1)

2

× n(n-1)

2

-matrix whose entries are in Z[q ±1 , t ±1 ]. For a Laurent polynomial a ∈ Z[q ±1 , t ±1 ], let M q (a) and m q (a) be the maximal and the minimal degree of the variable q, respectively. For a matrix A = (a ι,κ ) ∈ GL( n(n-1)

2

; Z[q ±1 , t ±1 ]) we define

M q (A) = max 1 ι,κ n(n-1) 2 {M q (a ι,κ )}, m q (A) = min 1 ι,κ n(n-1) 2 {m q (a ι,κ )}.
We denote the supremum, the infimum and the length function of a braid β for the dual Garside structure by inf Σ * (β), sup Σ * (β), and l Σ * (β), respectively -see Section 2.1 for precise definitions.

The main result in this paper is the following, which was already conjectured by Krammer in [START_REF] Krammer | The braid group B 4 is linear[END_REF].

Theorem 1.1 (Variable q in the LKB representation and the dual Garside length).

For β ∈ B n , (1) 2 sup Σ * (β) = M q (L(β)). (2) 2 inf Σ * (β) = m q (L(β)). (3) 2l Σ * (β) = max{0, M q (L(β))} -min{0, m q (L(β))}.

Acknowledgement

The first author was supported by the JSPS Institutional Program for Young Researcher Overseas Visits.

Preliminaries

First of all we set up our notation and conventions. Let D 2 = {z ∈ C||z| 1} be the unit disc in the complex plane and D n = D 2 -{p 1 , . . . , p n } be the n-punctured disc, where each puncture p i is put on the real line so that -1 < p 1 < p 2 < • • • < p n < 1 holds. The braid group B n is identified with the mapping class group of D n .

Thoroughout this paper, we adopt the following conventions. Braids acts on the left. A positive standard generator σ i is identified with the left-handed, that is, the clockwise half Dehn-twist which interchanges the punctures p i and p i+1 . Also, when calculating winding numbers, the positive direction of winding is taken as the clockwise direction. This convention is the opposite to Bigelow's.

2.1. The dual Garside length. First we review the definition of the supremum, the infimum and the length functions of the dual Garside structure. We do not need detailed Garside theory such as normal forms or algorithms. We only need a length formula, Proposition 2.1 below.

Recall that the braid group B n is presented as

B n = σ 1 , . . . , σ n-1 σ i σ j σ i = σ i σ j σ i , |i -j| = 1 σ i σ j = σ j σ i , |i -j| > 1
The elements of the generating set Σ = {σ 1 , . . . , σ n-1 } are called Artin generators or Garside generators.

In the dual Garside structure, we use a slightly different generating set which contains Σ as a subset. For 1 i < j n, let a i,j be the braid

a i,j = (σ i σ i+1 . . . σ j-2 )σ j-1 (σ i σ i+1 . . . σ j-2 ) -1
The generating set Σ * = {a i,j | 1 i < j n} was introduced in [BKL], and its elements are called the dual Garside generators, or band generators, or Birman-Ko-Lee generators. In our conventions, the braid a i,j is represented by the left-handed half Dehn-twist along an arc connecting p i and p j in the lower-half of the disc {z ∈ D 2 | Im z < 0}.

A dual-positive braid is a braid which is written by a product of positive dual Garside generator Σ * . The set of dual-positive braids is denoted by B + * n . The dual Garside element is a braid δ given by δ = a 1,2 a 2,3 . . . a n-1,n Let Σ * be the subword partial ordering with respect to the dual Garside generating set Σ * : β 1 Σ * β 2 if and only if β -1 1 β 2 ∈ B + * n . For a given braid β, the supremum sup Σ * (β) and the infimum inf Σ * (β) is defined by

sup Σ * (β) = min{m ∈ Z | β Σ * δ m } and inf Σ * (β) = max{M ∈ Z | δ M Σ * β} respectively.
A dual-simple element is a dual-positive braid s which satisfies 1 Σ * s Σ * ∆. The set of dual-simple element is denoted by [1, δ]. The Garside length l Σ * is the length function with respect to the generating set [1, δ].

The next formula relates the supremum, infimum and the length.

Proposition 2.1 (Dual Garside length). For a braid β ∈ B n we have the following equality

l Σ * (β) = max{0, sup Σ * (β)} -min{inf Σ * (β), 0}.
See [BKL] for proof.

2.2. The Lawrence-Krammer-Bigelow representation of the braid groups.

We review a definition of the Lawrence-Krammer-Bigelow representation. Our description of the LKB representation is homological and mainly follows Bigelow, but it is slightly modified so as to agree with our conventions. For details, see [Big].

Let C = {(z 1 , z 2 ) ∈ D n × D n | z 1 = z 2 }/S 2
be the configuration space of two unordered points of D n . We take base points d 1 = exp( 3 2 πi-ε) and d 2 = exp( 3 2 πi+ ε) in ∂D n so that d 2 lies on the right side of d 1 -see Figure 1. We take {d 1 , d 2 } as a base point of C. Let q, t be the free abelian group of rank two generated t and q, and define a homomorphism φ : π 1 (C) → q, t as follows: let {γ, γ ′ } : [0, 1] → C a loop representing an element x ∈ π 1 (C). We define the number a as the sum of the winding numbers along each puncture p j :

a = - 1 2πi n j=1 γ dz z -p j + γ ′ dz z -p j .
We define the number b as twice the winding number of the path γ -γ ′ (i.e. as twice the relative winding number of the points):

b = - 1 πi γ-γ ′ dz z
Now we define φ(x) = q a t b . (Here we add the minus signs since we adapted the convention that the positive winding direction is the clockwise direction.) Let π : C → C be the covering of C associated to ker φ. Fix the lift of the base point {d 1 , d 2 }. By abuse of notation, we use the same symbol {d 1 , d 2 } to represent the base point both in C and C. Then q and t are regarded as deck translations and the second homology group

H 2 ( C; Z) is a Z[q ±1 , t ±1 ]-module.
Let Y be the Y -shaped graph shown in Figure 1, having one distinguished external vertex r, two other external vertices v 1 and v 2 , and one internal vertex c. We orient the edges of Y as shown in Figure 1.

v 1 v 2 c r d 1 d 2 d 1 d 2 F F ′ N Figure 1. Fork and Noodle A fork F is an embedded image of Y into D 2 such that: • All points of Y \ {r, v 1 , v 2 } are mapped to the interior of D n .
• The distinguished vertex r is mapped to the base point d 1 .

• The other two external vertices v 1 and v 2 are mapped to two different puncture points.

• The edge [r, c] and the arc [v 1 , v 2 ] = [v 1 , c] ∪ [c, v 2 ] are both mapped smoothly. The image of the edge [r, c] is called the handle of F . The image of [v 1 , v 2 ] = [v 1 , c] ∪ [c, v 2 ],
regarded as a single oriented arc, is called the tine, denoted T (F ). The image of c is called the branch point of F .

A parallel fork F ′ is a fork as depicted by a dotted line in Figure 1:

F ′ is parallel to F and the distinguished vertex is mapped to d 2 .
A noodle is an oriented smooth embedded arc which begins at d 1 and ends at d 2 . Let γ, γ ′ : [0, 1] → D n be the handles of the forks F and of its parallel

F ′ , respectively. Let { γ, γ ′ } : [0, 1] → C be the lift of the path {γ, γ ′ } : [0, 1] → C taken so that { γ, γ ′ }(0) = {d 1 , d 2 }. Consider the surface Σ(F ) = {{x, y} ∈ C | x ∈ T (F ), y ∈ T (F ′ )} in C. Let Σ(F ) ⊂ C be the component of π -1 Σ(F ) which contains the point { γ, γ ′ }(1).
The surface Σ(F ) in C defines an element of H 2 ( C; Z). By abuse of notation, we will use the same symbol F to represent the 2nd homology class defined by the surface Σ(F ). In a similar way, a noodle N defines a surface Σ(N ) which defines an element of H 2 ( C, ∂ C; Z). Again, by abuse of notation we denote the homology class [ Σ(N )] by N .

For 1 i < j n, let us take a fork F i,j as shown in Figure 2. These forks are called standard forks. We call a standard fork of the form F = F i,i+1 a straight fork. Similarly, let N i be the noodle which encloses the i-th puncture point p i as shown in Figure 2. We call such a noodle a standard noodle. Bigelow showed that H 2 ( C; R) is a free R[q ±1 , t ±1 ]-module and {F i,j } is a basis of H 2 ( C; R). From now on, by using this basis we always identify

H 2 ( C; R) with R[q ±1 , t ±1 ] n(n-1)/2 . p i F i,j p j p i N i

Figure 2. Standard forks and standard noodles

The braid group B n acts on the covering C so that it commutes with the deck translations t and q. Hence we get a representation B n → GL(H 2 ( C; R)) = GL( n(n-1)

2

; R[q ±1 , t ±1 ]). The exact matrices are given in the Introduction. This is the Lawrence-Krammer-Bigelow (LKB) representation.

Remark 2.2. As shown in [PP], the standard forks do not form a basis of the Z[q ±1 , t ±1 ]-module H 2 ( C; Z). Thus in the above description it was important to use real coefficients, even if the matrix coefficients actually all lie in Z[q ±1 , t ±1 ].

Noodle-fork pairings. The Noodle-Fork pairing is a homology intersection pairing H

2 ( C; ∂ C)×H 2 ( C) → Z[q ±1 , t ±1 ].
As Bigelow showed in his so-called Basic Lemma ([Big, Lemma 2.3]), the pairing is calculated as follows.

Let z 1 , . . . , z m be the intersection points of T (F ) with N , and let z ′ i be the intersection of T (F ′ ) with N which corresponds to z i .

Observe that a pair of intersection points {z i , z ′ j } ∈ C corresponds to an intersection point of the surfaces Σ(N ) and Σ(F ). Hence it contributes to the total intersection pairing of N and F as a monomial ε i,j m i,j = ε i,j q ai,j t bi,j where ε i,j denotes the sign of the intersection at {z i , z ′ j }.

The monomial ε i,j m i,j is computed as follows. First we define c i,j by

c i,j = +1 if d 2 and z ′ j belong to the same component of N -z i . -1 if d 1 and z ′ j belong to the same component of N -z i .
Take three paths A, B and C in D n as follows:

• A is a path from d 1 to the branch point of F along the handle of F .

• B is a path from the branch point to z i along the tine T (F ).

• C is a path from z i to d k along the noodle N . Here k = 1 if c i,j = +1 and k = 2 if c i,j = -1. In other words, C goes along N starting from z i , choosing the direction so as to avoid z ′ j . Similarly, we take three paths A ′ , B ′ and C ′ by

• A ′ is a path from d ′ 1 to the branch point of F ′ along the handle of F ′ . • B ′ is a path from the branch point to z ′ j along the tine T (F ′ ). • C ′ is a path from z ′ j to d k ′ along the noodle N . Here k ′ = 2 if c i,j = +1 and k ′ = 1 if c i,j = -1.
Now the concatenation of the three paths {C, C ′ }{B, B ′ }{A, A ′ } defines a loop l in C. The monomial m i,j is given by φ(l). The sign of the intersection ε i,j is given by

ε i,j = -ss ′ c i,j
where s (respectively s ′ ) is the sign of the intersection of N and T (F ) (respectively T (F ′ )) at z i (respectively z ′ j ). In summary, the noodle-fork pairing is given by the following sum, where we recall that m denotes the number of intersections of T (F ) with N :

N, F = 1 i,j m ε i,j m i,j ∈ Z[q ±1 , t ±1 ].
By direct computations we observe the following, which will play an important role in the proof of main theorems.

Lemma 2.3. Let N be a standard noodle and let F be a standard fork. Then

M q ( N, F ) 1, m q ( N, F ) -1.
Proof. Here we give the calculation of the paring N, F , for the most complicated case, the other cases are treated similarly. Let F = F i,j and N = N k be a standard fork and noodle and assume that i < k < j. Then T (F ) and N intersect at two points, z 1 and z 2 hence the two surfaces Σ(N ) and Σ(F ) intersects at four points. Now a i,j , b i,j and ε i,j are calculated as shown in Figure 3. The paths CBA and C ′ B ′ A ′ are depicted by a black and gray line, respectively. Thus, we conclude N, F = qt + 1 -t -q -1 .

The wall crossing labeling of curve diagrams

In this section we introduce the wall crossing labeling on curve diagrams and show that this labeling reflects the dual Garside length of braids. This result is interesting in its own right. (1) The number of intersections of φ(E) with the walls is as small as possible within the diffeotopy class of φ.

z 1 z ′ 1 z 2 z ′ 2 N, F = qt + 1 -t -q -1 {z 1 , z ′ 1 }    a 1,1 = 1 b 1,1 = 1 ε 1,1 = 1 {z 2 , z ′ 1 }    a 2,1 = 0 b 2,1 = 1 ε 2,1 = -1 {z 1 , z ′ 2 }    a 1,2 = 0 b 1,2 = 0 ε 1,2 = 1 {z 2 , z ′ 2 }    a 2,2 = -1 b 2,2 = 0 ε 2,2 = -1
(2) Near the puncture points, the image φ(E) coincides with the real line.

We denote the curve diagram of β by D β and the total curve diagram by D β . See the right side of Figure 4 for an example.

To introduce the wall-crossing labeling and make the correspondence between fork and curve diagram explicit, we use a modified version of the curve diagrams. For each puncture p i other than β(p n ), we take a small disc neighborhood of p i , say B i , and let B = ∪ i B i . Around each puncture p i , we modify the curve diagram D β as shown in Figure 5. We denote the resulting (total) curve diagram by M D β (M D β ), and call it the (total) modified curve diagram. The right side of Figure 5 shows the total modified curve diagram of σ 1 . Notice that in Definition 3.2 we used the largest and smallest labels only of curve diagram D β , not the total curve digram D β . However, in order to determine the wall crossing labelings we need to consider the total curve diagram.

p i B i p i B i M D σ1
We now show that the dual Garside (Birman-Ko-Lee) length of a braid can be read off the wall crossing labeling of its curve diagram.

Theorem 3.3. For a braid β ∈ B n we have the following equalities:

(1) sup Σ * (β) = LWcr(β).

(2) inf Σ * (β) = SWcr(β). 2 σ 1 ) 2 , and its wall crossing labeling. Among the labels of the solid (not dashed) arcs, the smallest one is -2 and the largest one is 2.

(3) l Σ * (β) = max(LWcr(β), 0) -min(SWcr(β), 0).

Example 3.4. Let us consider the braid β = (σ -1 2 σ 1 ) 2 = (a -1 2,3 a 1,2 ) 2 -see Figure 6. Any word representing this braid with letters belonging to {a ±1 1,2 , a ±1 2,3 , a ±1 1,3 , δ ±1 } has at least two negative letters, so inf Σ * (β) = -2. Similarly, any word representing β using these letters has at least two positive letters, so sup Σ * (β) = 2. (Indeed, the dual normal form of β is a -1 1,2 .a -1 1,3 .a 1,2 .a 1,2 .) Now Theorem 3.3 asserts that the smallest and largest wall crossing labelings occurring in the curve diagram of β should be -2 and 2, respectively, and Figure 6 shows that this indeed the case.

Proof of Theorem 3.3. First of all, we show that for a dual-positive braid β the equality LWcr(β) = l Σ * (β) holds.

To treat the dual Garside structure, we temporarily isotope our curve diagram and walls so that the all punctures sit on the circle |z| = 1 2 and walls are disjoint from the subdisc |z| 1 2 as shown on the left side of Figure 7. This isotopy does not affect the wall crossing labelings, since the wall crossing labeling is defined by using algebraic intersection of arcs and walls.

As shown in [BKL], the set [1, δ] of dual-simple braids is in bijection with the set of disjoint collections of convex polygons in D n whose vertices are punctures. This bijection is given as follows: to any such collection of polygons we can associate a dance of the puncture points, moving each puncture which belongs to some polygon P in the clockwise direction along the boundary of P , to the position of the adjacent vertex. In this way, each dual-simple element can be represented by some disjoint convex polygons in D n -see Figure 7.

We remark that polygons may be degenerate, having only two vertices; in the associated braids, the two corresponding punctures are interchanged by a clockwise half Dehn-twist. Given a curve diagram D β , equipped with the wall crossing labeling, and given a collection P = {P 1 , . . . , P k } of disjoint polygons with vertices in the punctures representing some dual-simple braid β + , let us describe in detail how to obtain the curve diagram D β + β and its wall crossing labeling. Consider the disjoint collection of annuli A 1 , . . . , A k in D n as follows (see Figure 8):

(1) The outer and inner boundary component of A i are both homotopic in D n to the boundary of a regular neighborhood of P i .

(2) The boundary components of A i are in reduced relative position (no bigons) with respect to the walls W and also with respect to the curve diagram D β .

(3) Among all collections of annuli satisfying ( 1) and ( 2), we choose the one where the annuli contain as many intersection points D β ∩ W as possible.

Roughly speaking, we push as much twisting of the curve diagram around the polygons as possible into the annuli.

Let us denote the components of D n \ A i containing the polygon P i by N (P i ).

Note that the intersection of D β with A i is simply a spiral, and that the labels on this spiral interpolate linearly between the labels on the outer and inner boundary component. Now the β + -action on D β is simple to describe: on

D n \ (∪ k i=1 (A i ∪ N (P i ))
), the diagram and its labeling is unchanged. On N (P i ) the diagram is turned one notch in the clockwise sense, and all labels are increased by one. On the annuli, we have some twisting, but the labels still just interpolate -see Figure 8.

The action of a negative dual-simple braid β -is similar, the only difference being that the twisting is in the counterclockwise direction and labels are decreased by one.

From the above description for a dual-positive braid β and a dual-simple element β + , we have the inequality LWcr(β + β)

LWcr(β) + 1. This implies the inequality LWcr(β) l Σ * (β)
The converse inequality is now implied by the following lemma:

Lemma 3.5. The dual-positive braid β can be written as the product of LWcr(β) positive dual-simple braids.

P i A i 4 3 4 3 2 1 0 1 β + Figure 8. Action of dual-simple elements
Proof of Lemma 3.5. From the above description of the action of a dual-positive braid, we see that SWcr(β) 0, as no negative labels can ever be created from non-negative ones.

In order to prove the lemma, we distinguish two cases. Firstly, if SWcr(β) = LWcr(β), then no arc of the modified curve diagram M D β crosses a wall. This implies that β is a power of δ, more precisely, β = δ SWcr(β) , and the lemma is true.

Secondly, if 0 SWcr(β) < LWcr(β), then we proceed inductively. We shall construct a negative dual-simple braid β -such that LWcr(β -β) < LWcr(β) and SWcr(β -β) SWcr(β).

Consider the arc segments having the maximal wall crossing labeling. Each such segment connects two walls. Let I = {(i, j)} be the pairs of walls which are connected by some maximal labeled arcs. Let P be the minimal (with respect to inclusion) collection of convex polygons which contains all straight lines connecting p i and p j for (i, j) ∈ I. Let β -be the inverse of the dual-simple braid that corresponds to P . Garside theoretically speaking, β -= ( (i,j)∈I a i,j ) -1 .

According to our description above of the β --action on the curve diagram D β , the action of β -decreases all the maximal wall crossing labelings by one, so LWcr(β -β) = LWcr(β) -1.

On the other hand, we claim that SWcr(β -β) SWcr(β), i.e., contrary to the largest label, the smallest label does not decrease during the β --action.

Let us prove this claim. First we observe that, roughly speaking, minimally labeled arcs are S-shaped, whereas maximally labeled arcs are S-shaped. More precisely, when both endpoints of a minimally labeled arc lie in the interior of walls, then the initial and terminal segment of the arc lie on the counterclockwise sides of the wall, whereas maximally labeled arcs begin and terminate on clockwise sides of the respective walls (see Figure 9(a)). Now, in order to prove the claim, we have to rule out the existence of minimally labeled arcs (i.e. arcs labeled SWcr(β)) which, under the β --action, give rise to arcs with an even smaller label.

First we observe that a minimally labeled arc α cannot intersect the interior of any of the polygons of P . Indeed, assume that an arc segment α enters into the interior of one of these polygons, say P 1 . Let us assume in addition that both endpoints of α lie in the interior of walls, not in punctures. Then α cuts P 1 into two components, and thus separates the vertex punctures of P 1 into two families (drawn white and black in Figure 9(a)). Now we observe that no maximally labeled arc can connect a wall belonging to a black puncture to a wall belonging to a white puncture. But, by construction of P , that means that the white and black punctures do not belong to the same polygon P 1 , which is a contradiction. Similarly, when a minimally labeled arc α intersects the interior of P but has one or both endpoints in punctures, then the same argument applies, we only have to decide in which color to paint a puncture at the extremity of α. The choice which works is to group such a puncture with the following punctures in the clockwise direction -see Figure 9(a).

Finally, we have to consider a minimally labeled arc segment α which ends in a vertex of P 1 without intersecting the interior of P 1 . Here we have to distinguish two cases. If an extremal segment of α lies on the counterclockwise side of the wall corresponding to its terminal puncture, as in Figure 9(b), then there is nothing to worry about since the β --action does not decrease its Wcr-label. If, on the other hand, a terminal segment of α lies on the clockwise side of the wall corresponding to its terminal puncture, as in Figure 9(c), then this wall cannot be connected to a wall of any of the other punctures of P 1 by a maximally labeled arc. Again, this contradicts the construction of P 1 .

In summary, no minimally labeled arc can generate an even smaller label under the β --action. This proves the claim, and thus Lemma 3.5.

To summarize, we have now proved that for a dual-positive braid β, the following equality holds:

LWcr(β) = l Σ * (β) For a dual-positive braid β, inf Σ * (β) 0 so by Proposition 2.1 we conclude that

LWcr(β) = sup Σ * (β).
In a similar way, we prove the equality SWcr(β) = inf Σ * (β) for any dual-negative braid β.

In order to prove the same results for arbitrary braids, we recall that the dual Garside element δ acts as the clockwise 2π n -rotation of the n-gon with vertices in all punctures. Thus the curve diagram D δβ is obtained from D β simply by a 2π nrotation, and the wall crossing labeling on each arc segment of D δβ is obtained from the label of the corresponding arc of D β by adding one. On the other hand, by definition of inf and sup, left multiplication by δ increases both inf and sup by one. Therefore for a general braid β, we get an equality

LWcr(β) = LWcr(δ -inf Σ * (β) β) + inf Σ * (β) = sup Σ * (δ -inf Σ * (β) β) + inf Σ * (β) = sup Σ * (β).
A similar calculation also yields

SWcr(β) = inf Σ * (β).
for a general braid β. Finally, Proposition 2.1 implies the third equality claimed in Theorem 3.3. This completes the proof of Theorem 3.3.

Remark 3.6. In [W] the second author defined another labeling on the curve diagram called the winding number labeling, and proved the similar formula for the winding number labeling and the usual Garside length ([W, Theorem 2.1]). The proof of Theorem 3.3 given here is a direct generalization of the proof of [START_REF] Wiest | How to read the length from its curve diagram[END_REF]Theorem 2.1].

Noodle-fork pairing and wall crossing labeling

In this section we make the observation (in Lemma 4.1) that the wall crossing labeling reflects the exponents of the variable q in the noodle-fork pairing.

We need some technical preparation for this. We first show that for any straight fork F i,i+1 , and for any braid β, the tine T (β(F i,i+1 )) can in a natural way be equipped with the wall crossing labeling. For this, we will have to relate curve diagrams and forks.

Let us consider the part of the curve diagram D β that is the image of the line segment between the i-th and (i + 1)-st punctures. We identify this part of the curve diagram with β(T (F i,i+1 )), the image of the tine of the straight fork F i,i+1 . Moreover, up to moving the footpoint from -1 to d 1 , part of the modified curve diagram can naturally be regarded as the handle of β(F i,i+1 ), as shown in Figure 10. This identification also induces the desired wall crossing labelling on each arc segment of β(F i,i+1 ). Let N be a standard noodle, and let F be a straight fork with tine T (F ). Consider the image β(F ) of a straight fork F . From now on, we always assume that N is isotoped so that it intersects β(T (F )) minimally. For an intersection point z i ∈ β(T (F )) ∩ N , we denote the wall crossing labeling of the arc segment of β(T (F )) containing z i by Wcr(z i ). We recall that the union of the walls is denoted W . We also recall that the intersection pairing of N and β(F ) is N, β(F ) = 1 i,j m ε i,j q ai,j t bi,j where the sum ranges over all pairs of intersection points N ∩ β(T (F )).

Lemma 4.1. Let F be a straight fork and let N be a standard noodle. Let {z i } be the set of intersection points of β(T (F )) and N . Then a i,i , the exponent of q at the intersection point {z i , z ′ i } is given as follows.

a i,i = 2Wcr(z i ) + 1 if z i and d 1 belong to the same component of N -W 2Wcr(z i ) -1 if z i and d 2 belong to the same component of N -W.
Proof. Recall that a i,i is defined as the sum of the winding numbers of the paths CBA and C ′ B ′ A ′ around each puncture:

a i,i = - 1 2πi n j=1 CBA dz z -p j + C ′ B ′ A ′ dz z -p j
where A, A ′ , B, B ′ , C, C ′ are the paths defined in Section 2.3. See Figure 11. is equal to the algebraic intersection number of γ and W j . Thus, the fork part of the path {CBA, C ′ B ′ A ′ }, namely the subarc {BA, B ′ A ′ }, contributes to the total winding number a i,i by 2Wcr(z i ). Finally, we observe that the rest of the loop {C, C ′ } (the noodle N part of the path {CBA, C ′ B ′ A ′ }) contributes to the total winding number a i,i by +1 (respectively -1) if z i and d 1 belongs to the same (respectively to different) components of N -(Z ∪ W ) -see Figure 11. This completes the proof.

A A ′ B B ′ C C ′ z i z ′ i d 1 d 2 A A ′ B B ′ C C ′ z i z ′ i d 1 d 2 W W
Remark 4.2. For the winding number labeling introduced in [W], and for b i,i , the exponent of t at the intersection {z i , z ′ i }, one can get a similar formula by similar arguments. Thus, schematically speaking, we have the following correspondence among three objects in braid groups:

LKB representation Curve diagram Garside structure Variable t

Winding number labeling Usual Garside structure Variable q

Wall crossing labeling Dual Garside structure Unfortunately, we did not manage to reprove Krammers result from [START_REF] Krammer | Braid groups are linear[END_REF], i.e. the analogue of Theorem 1.1 for the variable t and the usual Garside length, using our geometric techniques.

The dual Garside length formula

In this section we prove Theorem 1.1. For monomials q a t b and q a ′ t b ′ we define the lexicographical ordering q,t by q a t b q,t q a ′ t b ′ if a < a ′ , or if a = a ′ and b b ′ . The next lemma is the crucial result in Bigelow's proof of faithfulness.

Lemma 5.1. [START_REF] Bigelow | Braid groups are linear[END_REF]Bigelow's Key Lemma 3.2 and Claim 3.4] Assume that a noodle N and a fork F have the minimal geometric intersection. Then all intersection points {z i , z ′ j } of N with F which attain the < q,t -maximal monomial m i,j in N, F have the same sign ε i,j .

Roughly speaking, Lemma 5.1 states that the < q,t -maximal contributions to N, F do not cancel.

Lemma 5.2. Let β ∈ B n be a braid. Let F be a straight fork F such that β(T (F )) contains an arc segment having the largest wall crossing labeling LWcr(β). Then there exists a standard noodle N such that M q ( N, β(F ) ) 2LWcr(β) + 1.

Proof. Throughout proof, we assume that β(T (F )) intersects each standard noodle N i minimally. It is sufficient to show that there is a standard noodle N i and an intersection point z ∈ β(T (F )) ∩ N i such that Wcr(z) = LWcr(β) and such that the two points z and d 1 lie in the same component of N i -(W i ∩N i ): by Lemma 4.1, the intersection point {z, z ′ } contributes to the pairing N i , β(F ) by ±q 2LWcr(β)+1 t b for some b ∈ Z, and then Lemma 5.1 completes the proof. Let γ be an arc segment of β(T (F )) whose wall crossing labeling is LWcr(β). Take a standard noodle N i which intersects γ in a point z. Assume that z and d 1

Figure 3 .Figure 4 .

 34 Figure 3. Calculation of the noodle-fork pairing: an example

Figure 5 .

 5 Figure 5. Modified curve diagrams

Figure 6 .

 6 Figure 6. The total curve diagram of the braid β = (σ -12 σ 1 ) 2 , and its wall crossing labeling. Among the labels of the solid (not dashed) arcs, the smallest one is -2 and the largest one is 2.

Figure 7 .

 7 Figure7. The left picture shows the basic curve diagram (bold lines) and the walls (thin dashed lines). The pictures in the middle and on the left illustrate the action of a dual-simple braid P given by two disjoint polygons P 1 and P 2 .

Figure 9 .

 9 Figure 9. Fat dashed lines indicate minimally Wcr-labeled arcs, and fat solid lines indicate maximally Wcr-labeled arcs

Figure 10 .

 10 Figure 10. Viewing a curve diagram as a union of tines of forks, and viewing initial segments of modified curve diagrams as tines.

Figure 11 .

 11 Figure 11. Contribution of the noodle part

Figure 12 .

 12 Figure 12. Intersections achieving q 2LWcr(β)+1 exist

lie on the two different components of N i -(W i ∩ N i ). If γ does not fall into the puncture p i , then we can find another intersection point z * of γ with N i that lies on the same component of N i -(W i ∩ N i ) as d 1 (Figure 12 (a)) and the proof is complete. Assume that γ falls into the puncture p i . Then γ must pass under the adjacent puncture p i+1 because otherwise we either find another arc segment having strictly larger wall crossing labeling (Figure 12 (b)) or contradict the hypothesis that N i and β(T (F )) have the minimal intersection (Figure 12 (c)). Then the standard noodle N i+1 and γ have an intersection z * with the desired property (Figure 12(d)).

Proof of Theorem 1.1. For a braid β, let us take a straight fork F and a standard noodle N in Lemma 5.2. Thus we have M q ( N, β(F ) ) 2LWcr(β) + 1. Let us write the fork β(F ) as the linear combination of the standard forks

Since F is a straight fork, x i,j is an entry of the matrix L(β). Now we have an equality N, β(F ) =

Since in Lemma 2.3 we observed that

where the second inequality follows from Lemma 5.2 and the last equality from Theorem 3.3. Thus, we get an inequality M q (L(β)) 2 sup Σ * (β).

On the other hand, for each dual-simple element s ∈ [1, δ] M q (L(s)) 2

holds. Hence we get the converse inequality M q (L(β)) 2 sup Σ * (β).

We conclude that M q (L(β)) = 2 sup Σ * (β) The proof of (2) is similar, and (3) follows from Proposition 2.1.


Naudet V Gourry

3D Electrical Resistivity Tomography to locate DNAPL contamination in an urban environment

This study presents results from electrical resistivity campaigns performed around a housing estate located downstream to an old coke tar site in France. This coke tar has been previously studied with geochemical analyses and geophysical prospections. The previous results have shown a possible migration of the contaminant plume toward the housing estate. As geophysical measurements are difficult to perform in such an urban environment, the electrical array has been deployed all around the housing estate with an innovative arrangement of surface electrodes in C-shape and acquisition geometries in order to achieve a real 3D imaging of the subsoil and locate the contaminant plume above the housing estate. The electrical resistivity data were inverted with the ERTLabTM 3D inversion software developed by Multi-Phase Technologies and Geostudi Astier. Results show, at the depth of the aquifer, a very conductive plume emanating from old tar ponds and a slag heap and spreading through the housing estate.

Introduction

Dense nonaqueous phase liquids (DNAPLs) are separate-phase hydrocarbon liquids that are denser than water, such as chlorinated solvents, wood preservative wastes, coal tar wastes, and pesticides. Groundwater can be contaminated by DNAPLs as they can be introduced into the subsurface by spills, uncontrolled releases and leaking storage and disposal facilities. Upon release, DNAPLs exist within the subsurface in either a free-phase or a residual form. The free-phase form spreads downward under the influence of gravity until they are distributed to such an extent that residual levels are attained and capillary trapping prevents further movement, or an impermeable layer is reached, at which point pooling and lateral spreading may occur. The residual form whereby becomes trapped in soil pores or rock fractures. Some of the more commonly-used techniques to locate DNAPL in the subsurface include use of monitoring wells, multi-level samplers, organic chemical analyses of soil samples at different depths, and soils borings or cone penetrometers to determine site stratigraphy. However, these techniques provide limited spatial information of the subsurface at discrete intervals leading to incomplete site assessments and inadequate remedial designs. In addition, some of them require drilling into the subsurface, which can result in the creation of new pathways for continued vertical migration of free-phase DNAPLs. Therefore, complementary tools given spatial information can help to reduce the ambiguity of interpretation and raise the probability of success. Non-invasive geophysical methods can be very useful as they survey large areas in a short time and at a lower cost than drilling techniques, they avoid infiltration of the contaminant to non-contaminated layers, and they enable more precise information between existing boreholes or piezometers, often separated by large distances. Over the last decades, a variety of geophysical methods have been proposed for the detection of DNAPLs in the shallow subsurface including ground-penetrating radar [START_REF] Daniels | Ground penetrating radar for the detection of liquid contaminants[END_REF][START_REF] Bermejo | Geophysical discovery of a new LNAPL plume at the former Wurtsmith AFB[END_REF][START_REF] Cassidy | GPR attenuation and scattering in a mature hydrocarbon spill: a modeling study[END_REF], resistivity methods [START_REF] Newmark | Monitoring DNAPL pumping using integrated geophysical techniques[END_REF][START_REF] Benson | Mapping groundwater contamination using dc resistivity and VLF geophysical methods-a case study[END_REF][START_REF] Atekwana | Biodegradation and mineral weathering controls on bulk electrical conductivity in a shallow hydrocarbon contaminated aquifer[END_REF][START_REF] Chambers | Hydrogeophysical imaging of deposit heterogeneity and groundwater Near Surface 2011 -17 th European Meeting of Environmental and Engineering Geophysics[END_REF], induced polarization [START_REF] Chambers | Characterisation of a NAPLcontaminated former quarry site using electrical impedance tomography[END_REF] and self-potential measurements (e.g., [START_REF] Minsley | Three-dimensional self-potential inversion for subsurface DNAPL contaminant detection at the Savannah River Site, South Carolina[END_REF]. In this paper, we present two electrical resistivity campaigns performed around a housing estate located downstream an old coke tar site. The first campaign consisted in twelve 2D profiles with an electrode spacing of 5 m acquired in 2006 and 2008. The second campaign consisted in 3D crossprofiles acquisitions with a C-shape configuration and an electrode spacing of 10 m all around the housing estate. The originality of this study relies on the innovative arrangement of surface electrodes for the 2010 campaign and the processing procedure used to increase the resolution of the resistivity survey and therefore to better determine the contaminant plume extension under the housing estate.

Site presentation

The studied site is located downstream a former coke factory (North of France) that operated from 1921 to 1973. During this period, the factory produced 100 000 tons of coke, 4 000 tons of crude tar, 1 100 tons of ammonium sulfate and 400 tons of benzoyl. Since closure, the surface structures have been demolished and it was estimated that 2 000 to 20 000 tons of acid tar, a Dense Non-Aqueous Phase Liquid (DNAPL), as well as 3.10 5 m 3 of process water containing phenols and ammonium sulfate, were stored in two tar ponds. These tar ponds and the slag heap (fig. 1) constitute the main contaminant sources. One of the tar pond was emptied in 1992 and the other in 2004. The geology of the site comprises a 4 to 6 m thick layer of quaternary silt and of tertiary clayey-sand formation due to chalk alteration above a 50 m thick chalk secondary formation. The chalk aquifer is present at a depth varying between 25 to 30 m depending on the season. The impermeable substratum consists of green shale. The regional groundwater flow is predominantly northward. Since 1991, the tar ponds impact on the groundwater has been evidence and a monitoring borehole array has been installed inside and around the old coke tar (fig. 1). Geochemical analyses show elevated concentration of PAH15 and BTEX around the tar ponds (fig. 1). The hydraulic conductivity of the aquifer is very weak (K=10 -6 m/s), but possible fractures inside the chalk have been mentioned and the pollution sources, made of acid tar, may also have weathered the lower part of the aquifer close to the source zone, improving the downward migration of coal tar.

Figure 1 Location of the studied area. The produced wastes of the coal factor were stored in a slag heap and the coal tar inside tar ponds mixed with acids. Red lines indicate the location of the twelve 2D profiles (5m electrode spacing) and blue lines the location of L-shape electrodes used to perform 3D acquisition (10 m electrode spacing). Table presents some geochemical analyses of June 2009.

A housing estate is located just at the north of the contaminated site in the direction of the groundwater flow and no borehole has yet been realized in this zone. In 2009, the piezometer PzNW has been impacted by coal tar, suggesting an extension of the contaminant plume (fig. 1). As many boreholes for agriculture are present in this zone with drinking water borehole tree kilometers away, it was crucial to determine the extension of the contaminant plume under the housing estate in order to help in the location of future boreholes.

Data acquisition and processing

Base on a previous 2D electrical resistivity measurements performed in 2006 and 2008 all around the housing estate, we performed in October 2010 a 3D electrical resistivity campaign. To achieve a real 3D imaging of the subsoil, the specific working conditions lead to assume non-conventional geometries of the electrodes layout to assess if the contaminant plume spreads under the housing estate. A number of 226 electrodes were installed all around the housing estate each 10 meters along 4 straight profiles (fig. 1). Acquisitions were performed using a C-shape geometry of the electrodes (the electrodes can be approximately divided into two branches). A pole-dipole configuration in common cable (transmitting and receiving electrodes belonging to the same branch) and cross-cables (all possible combinations where transmitting and receiving dipoles are not on the same line) has been used with the -infinite‖ electrode placed 2 km away in the north-west direction. After data filtering on the geometric factor (<15000) and on the rms error during an injection cycle (<3%) at least 3% of the data have been removed leading to 33 271 measurements used for the inversion. Data were processed using the ERTLab TM 3D inversion software developed by Multi-Phase Technologies and Geostudi Astier. This code uses a Finite Elements (FEM) approach to model the subsoil by adopting mesh of hexaedrons to correctly incorporate complex terrain topography. Therefore, data-sets collected 2D profiles, 3D surface grids, boreholes, surface-to-hole arrays, or any 2D or 3D four-electrode array can be inverted [START_REF] Zhou | Finite element three-dimensional direct current resistivity modelling: accuracy and efficiency considerations[END_REF]. The reliable and robust algorithm for 3D resistivity inversion assures accuracy and convergence, even when high levels of noise are present in the dataset. The inversion procedure is based on a smoothness constrained least-squared algorithm [START_REF] Labrecque | Occam's inversion of 3-D ERT data[END_REF] with a Tikhonov model regularization. Noise is appropriately managed using a data weighting algorithm [START_REF] Morelli | Advances in ERT inverse modelling[END_REF]. A 5% standard deviation estimate for noise was assumed to invert the data set with a robust inversion. The data-set was inverted by using a 10x10x20 meters mesh size and starting from an initial homogeneous model of 43 Ohm.m. The model comprises 150 000 cells that require around 5 hours on a standard laptop to be completed.

Results

Figure 2.a shows electrical resistivity map at the depth of the aquifer obtained with the twelve 2D electrical profiles acquired in 2006 and 2008 all around the housing estate. These data have been inverted with Res2dInv [START_REF] Loke | [END_REF]) and then interpolated in 3D with surfer. A conductive zone with electrical resistivity values below 20 Ohm.m was clearly identified below the contaminant sources (slag heap and tar ponds) with an extension towards the north. This interpolated map also suggested a possible spreading of the contaminant plume towards the housing estate. Nevertheless, these 2D data were not able to image below the housing estate. To achieve a real 3D imaging of the subsoil below the housing estate, a non-conventional geometry of the electrodes have been used in L-shape and Cshape in October 2010. As this is a non-common acquisition configuration, the effectiveness of the 3D coverage has been studied with a sensitivity analysis. Such sensitivity analysis allows also investigating the resolution resulting from the geometry of electrodes. Resolution is a complex function of numerous factors (e.g. electrode layout, measurement schedule, data quality, imaging algorithm, electrical conductivity distribution) and, in general, varies significantly across the image plane. The sensitivity analysis performed on some selected quadrupoles shows that, at the aquifer depth, a satisfactory degree of sensitivity can be achieved not only below the electrodes but also below the housing estate.

Conclusions

A contaminant plume due to a former coke factory that operated from 1921 to 1973 is suspected to spread towards a housing estate according to 2D electrical resistivity profiles performed in 2006 and 2008 all around the housing estate and geochemical analyses performed in 2009. The new electrical resistivity campaign performed in 2010 consisted in innovative 3D electrical resistivity acquisitions with a C-shape configuration of the electrodes installed all around the housing estate. This original electrodes arrangement and the processing procedure used to increase the resolution of the resistivity survey helped us to better determine the contaminant plume extension under the housing estate. These geophysical results are of high interest for the placement of a new monitoring borehole planned to be drilled in 2011 in the housing estate.

Figure 2

 2 Figure 2 a. Electrical resistivity distribution between 20-30 m depth obtained with a 3D interpolation of 2D inverted data acquired in 2006 and 2008. b. Electrical resistivity distributing at 20-30 m depth obtained with 3D inverted data acquired in October 2010.

Acknowledgements

The authors acknowledge the financial support of ADEME (Agence De l'Environnement et de la Maitrise de l'Energie) for this study (agreement n°0572C0064).


Kamil Adiloglu
email: kamil.adiloglu@inria.fr

Emmanuel Vincent
email: emmanuel.vincent@inria.fr

A GENERAL VARIATIONAL BAYESIAN FRAMEWORK FOR ROBUST FEATURE EXTRACTION IN MULTISOURCE RECORDINGS

Keywords:

We consider the problem of extracting features from individual sources in a multisource audio recording using a general source separation algorithm. The main issue is to estimate and propagate the uncertainty over the separated source signals, so as to robustly estimate the features despite source separation errors. While stateof-the-art techniques estimate the uncertainty in a heuristic manner, we propose to integrate over the parameter space of the source separation algorithm. We apply variational Bayes to estimate the posterior probability of the sources and subsequently derive the expectation of the features by moment matching. Experiments over stereo mixtures of three or four sources show that the proposed method provides the best results in terms of the root mean square (RMS) error on the estimated features.

INTRODUCTION

Many applications in the field of audio information retrieval are typically solved by extracting features describing the audio content and exploiting them for e.g. speech recognition, speaker recognition, cover version detection etc. However, most audio signals consist of a mixture of several sound sources, which have their own characteristics. Applying source separation prior to feature extraction can increase retrieval accuracy. For increased robustness, the uncertainty over the separated sources must be estimated in the complex-valued time-frequency domain and propagated to the features [START_REF] Kolossa | Independent component analysis and time-frequency masking for speech recognition in multitalker conditions[END_REF].

A heuristic approach is to assume that the uncertainty is proportional to the squared difference between the separated sources and the mixture [START_REF] Kolossa | Independent component analysis and time-frequency masking for speech recognition in multitalker conditions[END_REF][START_REF] Delcroix | Static and dynamic variance compensation for recognition of reverberant speech with dereverberation preprocessing[END_REF]. In [START_REF] Astudillo | A MMSE estimator in mel-cepstral domain for robust large vocabulary automatic speech recognition using uncertainty propagation[END_REF], a more principled approach is taken whereby the separated sources are assumed to follow a Gaussian posterior distribution, whose mean and variance are those of the Wiener filter used for separation. Propagation to the features is then achieved either by moment matching [START_REF] Astudillo | A MMSE estimator in mel-cepstral domain for robust large vocabulary automatic speech recognition using uncertainty propagation[END_REF] or unscented transform [START_REF] Kolossa | Independent component analysis and time-frequency masking for speech recognition in multitalker conditions[END_REF]. This approach remains mathematically inaccurate however, since the parameters of the Wiener filter are fixed to a certain value instead of being integrated over in a fully Bayesian approach.

In a preliminary study using a simple local Gaussian source model, we proposed a Gibbs sampling algorithm and a variational Bayes (VB) algorithm to address this integration and showed that the latter decreased the RMS error over the resulting Mel frequency cepstral coefficients (MFCC) [START_REF] Adiloglu | An uncertainty estimation approach for the extraction of source features in multisource recordings[END_REF]. In this paper, we extend this approach to the general modeling framework for source separation recently introduced in [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF]. This framework generalizes a wide class of existing source separation algorithms, including nonstationarity-based frequency-domain independent component analysis (FDICA) and single-or multi-channel nonnegative matrix factorization (NMF). We propose a VB algorithm to estimate the posterior distribution of the source time-frequency coefficients and subsequently derive the expectation of the features by moment matching. Figure 1 illustrates the workflow of the proposed approach.

This paper is organized as follows. Section 2 introduces the source separation framework. Section 3 presents the proposed VB inference algorithm for the estimation of the posterior distribution of the sources. Section 4 presents the uncertainty propagation method. In Section 5, we evaluate this framework over convolutive mixtures. We conclude in Section 6.

GENERAL SOURCE SEPARATION FRAMEWORK

Classically, we address the source separation problem in the timefrequency domain by means of the Short-Time Fourier Transform (STFT). For J sources and I channels, the mixing equation writes

x f n = A f s f n + ǫ f n , (1)
where

x f n = [x 1,f n , . . . , x I,f n ] T denotes the mixture STFT coef- ficients, s f n = [s 1,f n , . . . , s J,f n ] T the source STFT coefficients. A f = [A 1,f , . . . , A J,f
] is the complex valued mixing matrix, where A j,f is the vector of mixing coefficients for source j and ǫ f n is the noise. In this formulation, f is the frequency index, n the time frame index, i the channel index and j the source index. Note that this framework also works for diffuse or reverberated sources by modeling each source as a subspace spanned by several point sources [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF].

We adopt a local Gaussian model [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF] for the source coefficients. We set a zero-mean complex Gaussian prior over the source coefficients s j,f n with variance v j,f n :

s j,f n ∼ N (0, v j,f n ). (2)
The source variances v j,f n , which encode the spectral power are decomposed via an excitation-filter model [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF]:

v j,f n = v ex j,f n v ft j,f n . (3)
The excitation spectral power v ex j,f n is decomposed into characteristic spectral patterns modulated by time activation coefficients. Finally, the characteristic spectral patterns are defined as the sum of narrowband spectral patterns w ex j,f l with associated weights u ex j,lk . Similarly, the time activation coefficients are represented as a sum of time-localized patterns h ex j,mn with their weights g ex j,km . The same decomposition applies to the filter spectral power v ft j,f n . This framework makes it possible to incorporate a wide range of constraints about the sources. For instance, harmonicity can be enforced by choosing w ex j,f l as narrowband harmonic spectra and letting the spectral envelope and the active pitches be inferred from the data via the other parameters. For more details about how to constrain spectral and temporal structures, see [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF]. As a result, the complete factorization scheme is as follows:

v ex j,f n = K ex j X k=1 M ex j X m=1 L ex j X l=1 h ex j,mn g ex j,km u ex j,lk w ex j,f l , (4)
v ft j,f n = K ft j X k ′ =1 M ft j X m ′ =1 L ft j X l ′ =1 h ft j,m ′ n g ft j,k ′ m ′ u ft j,l ′ k ′ w ft j,f l ′ . (5)

PROPOSED VARIATIONAL INFERENCE ALGORITHM

For the separation of the sources, we propose a VB algorithm [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF].

The Bayesian Approach

Let us denote the set of all model parameters as

θ = {A f , h ex j,mn , g ex j,km , u ex j,lk , w ex j,f l , h ft j,m ′ n , g ft j,k ′ m ′ , u ft j,l ′ k ′ , w ft j,f l ′ }. (6)
For the sake of simplicity, we define all prior probabilities to be non-informative, even-though informative priors could be defined as well. Hence, A f follows a flat prior and the non-negative factors {h ex j,mn , g ex j,km , . . . , u ft j,l ′ k ′ , w ft j,f l ′ } follow a Jeffreys prior. Finally, we assume Gaussian noise ǫ f n with fixed variance σ 2 b I, where I is the identity matrix.

Under these assumptions, we aim to estimate the posterior probability of the source coefficients, which is given by

p(s|x) ∝ Z p(x|s, θ)p(s|θ)p(θ)dθ. ( 7
)
This integral is intractable. By contrast with previous approaches [START_REF] Kolossa | Independent component analysis and time-frequency masking for speech recognition in multitalker conditions[END_REF], which consist of estimating the maximum likelihood (ML) value θ for the parameters and considering p(x|s, θ), we resort to more accurate VB inference [START_REF] Bishop | Pattern Recognition and Machine Learning[END_REF].

Algorithm

VB minimizes the Kullback-Leibler (KL) divergence between the true posterior distribution p(s, θ|X) and some approximation q(s, θ), which is typically specified by assuming some factorization. Unfortunately, direct application of this rule does not yield a closed form solution. Therefore we pursue a decomposition of the sources as in [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF] for the above model.

A source s j,f n is sub-divided into Λj = Λ ex j Λ ft j sub-components such that Λ ex j = K ex j M ex j L ex j and Λ ft j = K ft j M ft j L ft j .
Let us define λ to be a joint index of {k, m, l} and λ ′ is a joint index of {k ′ , m ′ , l ′ }. We then decompose the sources as

s j,f n = X λ X λ ′ c j,f n,λ,λ ′ , (8)
where each sub-component c j,f n,λ,λ ′ follows a complex Gaussian distribution

c j,f n,λ,λ ′ ∼ N (0, v j,f n,λ,λ ′ ), (9)
with

v j,f n,λ,λ ′ = h ex j,mn g ex j,km u ex j,lk w ex j,f l h ft j,m ′ n g ft j,k ′ m ′ u ft j,l ′ k ′ w ft j,f l ′ . ( 10
)
The mixing equation for the source sub-components can be rewritten as

x f n = A + f c f n + ǫ f , (11)
where A + f indicates the extended mixing matrix where the elements a i,j,f of A f are repeated Λj times.

Assuming the following factorization of the variational approximation of the posterior

q(c, θ) = " Y f,n q(c f n ) "" Y f q(A f )) " " Y j,m,n q(h ex j,mn ) • • • Y j,f,l ′ q(w ft j,f l ) " , (12)
we obtain the optimal factors of c f n , A f and w ex j,f l as in the following. Due to space limitations, we give the update equations of w ex j,f l among all of the NMF components. However the rest of the factors follow a similar shape:

q * (c f n ) ∼ N (µ c,f n , R cc,f n ), (13)
q * (A f ) ∼ N (µ A,f , R AA,f ), (14)
q * (w ex j,f l ) ∼ Inv-Gamma(α ex w,j,f l , β ex w,j,f l ),

where A f denotes the column vector obtained by concatenating the rows of the mixing matrix A f and

R cc,f n = " B c,f n + (σ 2 b I) -1 X i ([µ A,f µ H A,f + R AA,f ]ii) T " -1 , (16)
µ c,f n = R cc,f n (µ + A,f ) H (σ 2 b I) -1 x f n , (17)
R AA,f = " 1 σ 2 b X n diag (R T s,f n , . . . , R T s,f n ) | {z } I times " -1 , (18)
µ A,f = Σ A,f ( 1 σ 2 b X n R xs,f n ), (19)
α ex w,j,f l = N K ex j M ex j Λ ft j , (20)
β ex w,j,f l = N 5 F 6 (K ex j ) 5 (M ex j ) 5 (L ex j ) 6 (Λ ft j ) 5 X n,km,λ ′ (µ c,f n ) 2 (j,kml,λ ′ ) + (R cc,f n ) (j,kml,λ ′ ),(j,kml,λ ′ ) β ex h,j,mn β ex g,j,km β ex u,j,lk β ft h,j,m ′ n β ft g,j,k ′ m ′ β ft u,j,l ′ k ′ β ft w,j,f l ′ . (21)
In ( 16), [•]ii denotes the diagonal J × J block corresponding to channel i. In the same equation, the matrix B c,f n writes B c,f n = diag( N 6 F 6 (Λ ex j ) 6 β ex w,j,f l β ex u,j,lk β ex g,j,km β ex h,j,mn

(Λ ft j ) 6 β ft w,j,f l ′ β ft u,j,l ′ k ′ β ft g,j,k ′ m ′ β ft h,j,m ′ n ) J,Λ ex j ,Λ ft j j,kml,k ′ m ′ l ′ . (22)
In ( 19), R xs,f n is given by

R xs,f n = [x 1,f n µ H s,f n , . . . , x I,f n µ H s,f n ] T . (23)
Finally, the posterior distribution of the source coefficients is calculated by summing up the corresponding elements of the mean and the covariance of the source sub-components as follows:

µ s,j,f n = X λ X λ ′ µ c,j,f n,λ,λ ′ , ( 24
) (R ss,f n ) j,j ′ = X κ X λ ′ (R cc,f n ) (j,λ,λ ′ ),(j ′ ,λ,λ ′ ) . ( 25
)
and the second raw moment of the source coefficients R s,f n is given by

R s,f n = µ s,nf µ H s,nf + R ss,f n . ( 26
)
The VB update equations of the scale parameters of the NMF components {β ex w,j,f l , . . . , β ft h,j,m ′ n } turn out to be identical to the state-of-the-art expectation maximization (EM) updates of the NMF components { ŵex j,f l , . . . , ĥft j,m ′ n } [START_REF] Ozerov | Multichannel nonnegative matrix factorization in convolutive mixtures for audio source separation[END_REF] up to a variable change. For β ex w,j,f l this variable change writes

ŵex j,f l = β ex w,j,f l N K ex j M ex j K ft j M ft j L ft j . (27)
Therefore, we replaced these EM updates by the multiplicative updates in [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF], which converge much faster.

Equations ( 16) to (27) depend on each other. After proper initialization, we cycle through these equations by replacing the dependent values with their new estimates.

UNCERTAINTY PROPAGATION

We now present a moment matching method for propagating the uncertainty over the source components to the source images and in a following step to the MFCC features.

Uncertainty Propagation for the Source Images

Due to phase and scale indeterminacies in the source estimates s j,f n , we use the spatial source images y j,f n = A j,f s j,f n instead for our experiments, which do not suffer from such indeterminacies [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF].

Once the posterior distribution of the source coefficients s f n has been computed, the posterior distribution of the source images is calculated by propagating the first two moments of the sources to the source images as follows:

µ y,j,f n = µ post A,f µ post s,f n , (28)
R yy,j,f n = " X rr ′ ((R AA,f ) (ir,i ′ r ′ ) + (µ A,f ) (ir) (µ A,f ) (i ′ r ′ ) ) ((R ss,f n ) (r,r ′ ) + (µ s,f n ) (r) (µ s,f n ) H (r ′ ) ) " II ii ′ -µ y,j,f n µ H y,j,f n . (29)

Uncertainty Propagation for Feature Extraction

We calculate the expectation of the MFCCs for each source as

µ MFCC jn = Z MFCC(yj1n)P (yj1n) dyj1n (30)
where yj1n = [y j,1f n ] f =1...F are the STFT coefficients of the first channel of source image j in time frame n. Deterministic calculation without the use of the uncertainty model simply yields MFCC(yj1n) = 20 D log 10 (M|yj1n|). In this formulation, D is the DCT matrix and M is the matrix containing the mel filter coefficients. Note that we chose the scaling so that the MFCCs are expressed in decibels (dB).

In order to obtain the mean estimate of the MFCCs µ MFCC jn , we first compute the mean and variance of |yj1n| using formulae for the Rice distribution in [START_REF] Kolossa | Independent component analysis and time-frequency masking for speech recognition in multitalker conditions[END_REF] and then propagate them through the logarithm using the moment matching formulae in [START_REF] Gales | Model-Based Techniques for Noise Robust Speech Recognition[END_REF]. Note that this also makes it possible to estimate the variance of the MFCCs and exploit it for e.g. classification tasks [START_REF] Kolossa | Independent component analysis and time-frequency masking for speech recognition in multitalker conditions[END_REF]. In the following, due to space limitations, we simply consider the mean.

EXPERIMENTAL EVALUATION

Data and Algorithmic Settings

We considered the development dataset of the 2008 Signal Separation Evaluation Campaign (SiSEC) 1 . This dataset contains synthetic and live recorded convolutive, under-determined, stereo mixtures. There are 32 mixtures of 3 sources and 24 mixtures of 4 sources. Each mixture has a duration of 10 s.

We performed experiments with eight different sets of constraints over the parameters as considered in the experiments section of [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF]. These scenarios consist of all combinations of the following possibilities:

• Rank: Each source is either a single point source (1) or a subspace spanned by two point sources [START_REF] Delcroix | Static and dynamic variance compensation for recognition of reverberant speech with dereverberation preprocessing[END_REF].

• Spectral Structure: The narrowband spectral patterns w ex j,f l are either unconstrained (un) or fixed (co) to harmonic and noise-like patterns.

• Temporal Structure: The time-localized patterns h ex j,mn are either unconstrained (un) or fixed (co) to decreasing exponential patterns.

All other parameters are free. We initialize the mixing matrix Aj using the direction of arrival (DOA) estimation algorithm proposed in [START_REF] Blandin | Multi-source TDOA estimation using SNR-based angular spectra[END_REF]. The noise variance σ 2 b is initialized to 10 -6 . Finally, we performed 200 iterations for convergence.

Evaluation Criteria

In order to assess the impact of source separation on feature extraction, we evaluate the proposed algorithm according to both tasks. Source separation quality is evaluated in terms of the Signal-to-Distortion Ratio (SDR) in [START_REF] Vincent | Performance measures in blind audio source separation[END_REF] between the mean of the estimated source images µ y,j,f n and the true source images.

Feature extraction accuracy is evaluated in terms of the RMS error [START_REF] Adiloglu | An uncertainty estimation approach for the extraction of source features in multisource recordings[END_REF] between the estimated µ MFCC jn and the true MFCCs.We ignore the first MFCC coefficient and consider the MFCCs 2 to 20 only.

1-un-un 2-un-un 1-co-un 2-co-un 1-un-co 2-un-co 1-co-co

Results

Table 1 shows the source separation performance of VB compared to that of the state-of-the-art ML method in [START_REF] Ozerov | A general flexible framework for the handling of prior information in audio source separation[END_REF]. ML and VB perform similarly in almost all of the eight configurations. However VB is 0.2 dB better than the state-of-the-art ML method and yields the best performance for the 1 point source, spectrally and temporally constrained model (1-co-co) with 2.5 dB SDR. The baseline binary masking method [START_REF] Blandin | Multi-source TDOA estimation using SNR-based angular spectra[END_REF] yields 0.95 dB SDR.

Table 2 shows the total RMS error in dB over the MFCCs obtained either by deterministic computation or by moment matching for both VB and ML algorithms. As one can see, VB based MFCC estimation performs 0.05 dB better than the ML based estimation. Besides, the moment matching method outperforms deterministic MFCC estimation in all configurations with around 0.9 dB. Again, VB performs best for the RMS error by 6.54 dB for the 1 point source, spectrally unconstrained, temporally constrained model (1un-co). The baseline binary masking method [START_REF] Blandin | Multi-source TDOA estimation using SNR-based angular spectra[END_REF] performs significantly worse than both VB and ML algorithms and yields 17.25 dB RMS error.

CONCLUSION

In this paper, we presented a general, fully Bayesian source separation algorithm and an uncertainty propagation algorithm for the computation of the expectation of the MFCCs of individual sources in multisource recordings.

This algorithm provides a fundamental breakthrough towards mathematically rigorous estimation of uncertainty for robust feature extraction. The resulting MFCC coefficients are slightly more accurate than those obtained via the standard ML method, proving that the ML method provides a reasonable approximation, but that it is possible to obtain more accurate estimates.

In the future, we will seek to improve the tightness of the variational bound (and thereby the accuracy of the resulting features) and exploit the variance of the estimated features for uncertainty decoding.

  This work was supported by the Quaero Program funded by Oseo.

Fig. 1 .

 1 Fig. 1. Flow of the proposed Bayesian source separation and feature extraction approach.

Table 1 .

 1 SDR in dB achieved by ML or VB source separation over all mixtures.

	2-co-co

Table 2 .

 2 Total RMS error in dB for the MFCCs 2-20 obtained by ML-or VB-based source separation followed by deterministic (det) or moment matching (mm) feature extraction over all mixtures.

http://sisec2008.wiki.irisa.fr/tiki-index.php?page=Under-determined+speech+and+music+mixtures


Patrick Richard

Sean Mcnamara

Merline Tankeo

On the relevance of numerical simulations to booming sand

Keywords: numbers: 45.70.Ht, 46.40.-f, 91.60.Lj

We have performed a simulation study of 3D cohesionless granular flows down an inclined chute. We find that the oscillations observed in [L. E. Silbert, Phys. Rev. Lett., 94 098002 (2005)] near the angle of repose are harmonic vibrations of the lowest normal mode. Their frequencies depend on the contact stiffness as well as on the depth of the flow. Could these oscillations account for the phenomena of "booming sand"? We estimate an effective contact stiffness from the Hertz law, but this leads to frequencies several times higher than observed. However, the Hertz law also predicts interpenetrations of a few nanometers, indicating that the oscillations frequencies are governed by the surface stiffness, which can be much lower than the bulk one. This is in agreement with previous studies ascribing the ability to sing to the presence of a soft coating on the grain surface.

A "booming" or "singing" dune is a sand dune that emits a loud sound when an avalanche occurs on its slip face. The sound can be very loud -audible up to 10 km away -and has a well defined frequency of order one hundred Hertz. The physical origins of this phenomenon are still matter of debate in spite of much experimental and theoretical work [1][2][3][4][5][START_REF] Mills | Powders and Grains[END_REF]. A successful theory of booming sand must explain why the sound is not a mix of a wide range of frequencies, and therefore several frequency selection mechanisms have been proposed. We now summarize them non-exhaustively. According to Andreotti [5], for a given grain size, the frequency is set by the shear rate of the shear band separating the avalanche to the static part of the dune. This first explanation is not compatible with the one proposed by Douady et al. [3], in which the frequency is set by a resonance of the flowing layer. A third explanation has been proposed by Vriend et al. [4]: The frequency is not selected by the properties of the avalanche but by the acoustical resonance induced by the stratification of the dune, explaining why the frequency may vary with time.

More recently [2], Andreotti and Bonneau show that the shear band separating an avalanche from the static part of the dune induces an amplification of guided elastic waves, leading to a linear instability. The frequency is then set by the maximum value of the instability growth rate. Finally Mills and Chevoir [START_REF] Mills | Powders and Grains[END_REF] made the interesting remark that something like booming sand had already been observed in numerical simulations [START_REF] Silbert | [END_REF] that exhibited spontaneous oscillations in granular flows near the angle of repose. These oscillations were interpreted as signs of intermittency near the jamming transition.

Here, we examine Mills and Chevoir's re-interpretation of these oscillations as an acoustical phenomena -possibly "booming sand". We extend the simulations performed in [START_REF] Silbert | [END_REF] by carrying out an extensive study of the influence of the contact normal stiffness and of the height of the flow. After a description of the method, we will show that the observed oscillations are multi-body harmonic oscillations. Their frequency is thus governed by the stiffness of the springs that model the repulsive inter-granular forces. We estimate a reasonable value for this quantity and discuss the relevance of the simulations.

Method. We use our own implementation [8] of the classical "Discrete Element Method" method where Newton's equations of motion for a system of N "soft" grains are integrated. This requires giving an explicit expression for the forces that act between grains. Such a technique is able to reproduce successfully experimental results for gravity driven flows [9][10][11][12][13], sheared systems [14], granular materials close to jamming [15], silos [16] or rotating drums [17][18][19][20]. The Discrete Element Method is well known and can be found in many papers [9][10][11][12][13][14][15][16][17][18][19][20][21]. Therefore, we just present here the forces used in this work (and also in [START_REF] Silbert | [END_REF]). For the normal force between two overlapping spheres we use a standard linear spring-dashpot interaction model [22]:

f n = k n δ n -γ n v n
, where δ n is the normal overlap, k n is the spring constant, γ n the damping coefficient and v n the normal relative velocity. The damping models the dissipation characteristic of granular materials. Likewise we model the tangential force as a linear elastic and linear dissipative force in the tangential direction: f t = -k t δ t -γ t v t , where k t is the tangential spring constant, δ t the tangential overlap, γ t the tangential damping and v t the tangential velocity at the contact point. The magnitude of δ t is truncated as necessary to satisfy Coulomb law: |f t | ≤ µ |f n |, where µ is the grain-grain friction coefficient.

Numerical set-up. As in [START_REF] Silbert | [END_REF], we simulate gravity-driven chute flow. In every way the parameters are the same as in [START_REF] Silbert | [END_REF] except we vary the normal stiffness between grains and the height of the flow. The grains are monosized (diameter d, mass m). Unless otherwise specified the number of grains is N = 8000. The chute consists in a 3D cell whose base is flat and rectangular with size 20d × 10d. It can be inclined relative to the horizontal by an angle θ (angle between the horizontal and the long axis of the base) and is periodic in the directions tangent to the base. The bottom of the cell is obtained by pouring under gravity g a large number of grains in the cell (θ = 0 • ) and by fixing those that are in contact with the base. This disordered layer of fixed grains is sufficient to prevent crystallization throughout the system. Before any measurements are taken, the inclination is increased to θ ≈ 30 • , causing a rapid flow that erases any influence of the initial state. The angle θ is then set to a final value that is cited in the captions or text below. The following values of the parameters are used: 2 × 10 4 ≤ k n d/mg ≤ 2 × 10 6 , k t = 2k n /7, γ t = 0 and µ = 0.5. The value of γ n is adjusted to obtain a normal restitution coefficient e n = 0.88 [9]. We use dimensionless quantities by measuring distances, times, and elastic constants respectively, in units of d, d/g and mg/d.

Numerical Results. In accord with previous results [START_REF] Silbert | [END_REF], we observe an oscillation near the angle of repose which is defined, for a given flow height, by the angle below which the flow stops. This motion can be identified through measurements of the total kinetic energy. Note that in all our simulations the rotational kinetic energy is much lower than the translational one. As reported in Fig. 1a, the kinetic energy displays regular oscillations (characteristic frequency ≈ g/d) and irregular fluctuations with lower frequency (≈ 0.1 to 0.2 g/d). We next turn our attention to the spring stiffness k n by studying its effect on the amplitude and on the frequency of the oscillations. As reported in Fig. 1b, the oscillations indeed depend on the spring stiffness. As k n increases, their frequency increases while their amplitude decreases. Fig. 2 shows the frequency and amplitude of these oscillations versus spring stiffness. Frequency scales with √ k n , whereas amplitude scales with 1/k n . In the limit k n → ∞, the oscillations disappear. However, the avalanche does not disappear in this limit, indicating that the avalanche and the oscillations are two separate processes. For this reason, we do not think that the oscillations should be interpreted as arising from the jamming transition, or as avalanche precursors. Now let us consider the oscillation frequency. It scales as 1/t c where t c is the two-body collision time. Assuming that a two particle collision is a demi-cycle of a damped harmonic oscillator leads to t c = π

/[2k n /m -(γ n ) 2 /m 2 ] 1/2 ≈ π/(2k n /m) 1/2 ∝ 1/ √ k n .
This scaling of the frequency suggests that the oscillations are harmonic vibrations of a normal mode. Let us estimate the frequency of the vertical normal modes. In the simulations, there are a certain number n L of layers of grains resting on the bottom of the chute. We model the granular bed as a one dimensional chain of n L masses connected by linear springs. This model is the linear harmonic chain, used in solid state physics as a very elementary model of phonons [START_REF] Ashcroft | Solid State Physics[END_REF]. The damping added to the particle interactions affects short wavelength vibrations most strongly, but long wavelength ones only weakly [START_REF] Mcnamara | Discrete-element Modeling of Granular Materials[END_REF]. Thus the motion is dominated by the longest possible wavelength. If the bottom of the chute is considered to be fixed, and the top surface is free, the longest wavelength is four times the depth of the layer. This leads to a frequency f nL / g/d ≈ (1/4n L ) k n d/mg. With n L = 40 and k n d/mg = 2 × 10 6 , one obtains f 40 ≈ 9 g/d. This value is twice as large as the one shown in Fig. 2. The difference between the prediction and the simulation is probably caused by the proximity of the stability threshold. The macroscopic stiffness of a granular packing decreases as a yield condition is approached. Here, the yield condition (i.e the avalanche) is reached when θ equals the angle of repose. To support this idea, simulations were carried out using exactly the same protocol but with a final angle of θ = 0 • . The system comes rapidly to rest, with the longest-lived motions being persistent oscillations at frequencies closer to the predicted values. But the central point is that the model reproduces the correct scaling with stiffness f nL ∼ √ k n . Our results also predict that frequency should diminish as the granular layer is made deeper: f nL ∝ 1/n L . To check our analysis we carried out numerical simulations with k n d/mg = 2 × 10 6 , θ = 20 • and for several number of grains (8000 < N < 32000). The results shown in Fig. 3 confirm the predicted scaling. We conclude, therefore, that the oscillations observed in the simulations are simply harmonic oscillations of the lowest normal mode.

Relevance of the numerical results. Up to this point, we have discussed our results solely within the framework of the very idealized model where a relatively small number of perfect spheres flow down a fixed inclined plane, interacting via linear spring forces. The linear force law allowed us to reproduce the oscillations under exactly the same conditions as Silbert [START_REF] Silbert | [END_REF], and also facilitated the analysis of the observed frequencies. The idealized context of our work is emphasized by the exclusive use of dimensionless units. We did not try to give physical values to any parameter such as the particle diameter d. For the remainder of the paper, we will investigate the suggestion of Mills et al. [START_REF] Mills | Powders and Grains[END_REF] and discuss the relevance of our findings to "booming sand". This means that we must assign physical values to all quantities. The most troublesome (and the most important) parameter is the spring stiffness k n . Real sand grains have non-linear force-displacement laws that cannot be characterized by a single spring constant. We will, nevertheless, try to straddle the difference between the model and physical system by choosing a single value of k n relevant to the oscillations. This approach may be problematic: Modifying k n affects the rheology of the flow [13], even for the highest values of k n used in this work [12]. To determine the appropriate value of k n , we will consider a slightly less idealized model where the grains are spheres made out of an isotropic elastic material with Young modulus E and Poisson ratio ν. The contact force between two such spheres is given by the Hertz law:

F n = E √ 2d 3(1 -ν 2 ) δ 3/2 . (1)
We will take k n to be the stiffness seen by the grains when they oscillate at low amplitude about their equilibrium positions. We have thus

k n = dF n /dδ = (3/2)F n /δ.
Rearranging this equation, we obtain

F n = 2k n δ/3. (2)
In this equation, k n is not a constant: k n ∝ δ 1/2 . We seek a typical value of k n that will be determined by the 4) and ( 5) for a variety of diameters d (in meters) and number nL of layers of particles. typical contact force. We suppose that this force is of the order of the weight of a column of n L grains:

F n = n L mg = πn L ρgd 3 /6, (3)
where ρ is the density of the material making up the grains. Eqs. ( 1), (2), and (3) involve the unknowns F n , k n , and δ. Combining them yields to

k n d mg = 3 2 n L d δ = 3 2 2 √ 2n L E π(1 -ν 2 )ρdg 2/3 . (4)
We will consider modifying n L and d, while sweeping all the other constants into a single parameter that depends on the material:

k n d mg = K n 1/3 L d 2/3 , K = 3 2 2 √ 2E π(1 -ν 2 )ρg 2/3 . (5)
Putting in values of appropriate for glass: g = 10m/s 2 , ρ = 2.4 × 10 3 kg/m 3 , ν = 0.2, and E = 50 GPa, one obtains K = 2.3 × 10 4 m 2/3 . Various values of the dimensionless stiffness k n d/mg, are shown in Fig. 4. Note that at d ≈ 1 mm, even the softest contact (those that support the weight of only one grain) have a stiffness k n d/mg ≈ 10 6 , i.e., equal to the half of the highest value in Fig. 1. Diameters typical of "booming dunes" -around 200 µm -lead to even stiffer contacts. Now let us see whether the oscillation frequencies in the simulation correspond to those of "singing sand". The unit of frequency used here is g/d ≈ 220 Hz for d = 200 µm (which is, mentioned before, the typical diameter for grains of booming dunes). The observed frequencies in the field and in laboratory are around 90 Hz, i.e. 0.4 g/d (see Fig. 3 of [4] or table I of [3]). According to laboratory experiments [3], flow heights are of order of several centimeters. For d = 200 µm, this corresponds to H/d ≈ 500. Extrapolating the data in Fig 3 to H/d = 500 leads to frequencies tantalizingly close to those observed. But Fig. 3 was obtained for a fixed stiffness k n d/mg = 2×10 6 , much lower than the value k n d/mg = 5.3 × 10 7 predicted by Eq. ( 5) for n L = 500 and d = 200 µm. If the oscillations are indeed the origin of booming sand, the effective stiffness of the grains must be much lower than predicted by Eq. ( 5). To explain this disagreement, let us determine the overlap δ between two contacting grains from Hertz theory which assumes those grains as perfect spheres, without any surface asperities. Eqs. ( 1) and (3) yield to

δ = 3d 5/3 n 2/3 L 2K . (6)
For d = 200 µm and n L = 500, we obtain δ ≈ 3 nm. This is a maximum value of δ, concerning contacts that support the weight of 500 grains. Contacts near the free surface of the flow (n L < 5) have δ < 2 Å. For such overlaps, the contact between two real grains will be dominated by surface properties that might be quite different from the bulk properties considered in the Hertz model. In particular asperities or a layer of silica gel (as proposed in [START_REF] Goldsack | [END_REF]) could significantly reduce the stiffness seen by acoustical waves. These speculations, however, can only be confirmed by examining the sand grains themselves. These surface effects could also account for the rapid increase in sound speed with depth [4]. The contact stiffness near the dune surface would be anomalously soft and fixed by grain surface properties, whereas deeper in the dune, at higher contact forces, the stiffness would be much higher and dominated by the bulk properties.

Conclusion.

We have examined the oscillations observed [START_REF] Silbert | [END_REF] in numerical simulations of granular beds near the angle of repose. These oscillations are harmonic vibrations of the lowest vertical mode of the bed, and their frequency obeys the expected dependency on particle stiffness and bed depth. They are not part of the avalanche motion but may be connected to "booming dunes", if the effective contact stiffness is about 20 times smaller than expected from the Hertz contact law. Such a reduction of stiffness is possible because the Hertz law predicts extremely small grain overlaps, indicating that the stiffness is dominated by surface properties instead of bulk ones. If the oscillations are indeed related to "booming sand", it would mean that the sound originates from a resonance inside a flowing layer, similar to the explanation presented in [3]. However, several issues remain open. For example, what excites the oscillations? One clue is given by the effect of polydispersity: we performed simulations with uniform size distribution of width 2dσ, with σ = 0.05, 0.1 or 0.2, and observe that the oscillations disappear for σ = 0.2, consistent with the observation [26] that booming sand has a narrow grain size distribution. This suggests that the resonance is excited by quasi-periodic collisions in the shearing layer, consistent with most explanations that have been presented.

6 FIG. 1 :

 61 FIG. 1: (color online)(a) Average kinetic energy per particle for an angle of inclination θ = 19.5 • and k n d/mg = 2 × 10 6 . Fluctuations (main panel) at large time scales and oscillations at small time scales can be observed (inset). (b) Average kinetic energy per particle for θ = 20 • (the angle of repose is 19.1 • ) and for different values of kn.

FIG. 2 :

 2 FIG.2:(color online) (a) Frequency (cycles per time unit) of the kinetic energy oscillations (•), corresponding 1/2 power fit (dashed line), amplitude of the corresponding Fourier component ( ) and corresponding (k n ) -1 fit (full line). The angle of inclination is θ = 20 • . (b) same curves in a log-log plot.

2 )FIG. 3 :

 23 FIG. 3: (color online) Frequency (cycles per time unit) of the kinetic energy oscillations (•) versus the height of the flow and corresponding fit (dashed line). The angle of inclination is θ = 20 • and the normalized stiffness k n d/mg = 2 × 10 6 .

5 FIG. 4 :

 54 FIG.4: (color online) The dimensionless stiffness of the contact as estimated by Eqs. (4) and (5) for a variety of diameters d (in meters) and number nL of layers of particles.

Acknowledgments. We thank J.T. Jenkins, D. Tan for stimulating discussions. This work is supported by the ANR project STABINGRAM No. 2010-BLAN-0927-01 and by the région Bretagne (CREATE SAMPLEO). M. T. is supported by the région Bretagne (ARED grant).