Athanase Papadopoulos
email: athanase.papadopoulos@math.unistra.fr

Weixu Su
email: suweixu@gmail.com

ON THE FINSLER STUCTURE OF THE TEICHM ÜLLER METRIC AND THURSTON'S ASYMMETRIC METRIC

Keywords: AMS Mathematics Subject Classification: 32G15, 30F60, 57M50, 57N05 Teichmüller space, Thurston's asymmetric metric, Teichmüller metric, Finsler norm

We highlight several analogies between the Finsler (infinitesimal) properties of Teichmüller's metric and Thurston's asymmetric metric on Teichmüller space. Thurston defined his asymmetric metric in analogy with Teichmüllers' metric, as a solution to an extremal problem, which consists, in the case of the asymmetric metric, of finding the best Lipschitz maps in the hoomotopy class of homeomorphisms between two hyperbolic surface. (In the Teichmüller metric case, one searches for the best quasiconformal map between two conformal surfaces.) It turns out also that some properties of Thurston's asymmetric metric can be used to get new insight into Teichmüller's metric. In this direction, in analogy with Thurston's formula for the Finsler norm of a vector for the asymmetric metric that uses the hyperbolic length function, we give a new formula for the Finsler norm of a vector for the Teichmüller metric that uses the extremal length function. We also describe an embedding of projective measured foliation space in the cotangent space to Teichmüller space whose image is the boundary of the dual of the unit ball in the tangent space representing vectors of norm one for the Finsler structure associated to the Teichmüller metric.

Introduction

Let S = S g,n be an oriented surface of genus g ≥ 0 with n ≥ 0 punctures. We assume that 3g -3+n > 0. The Teichmüller space T(S) of S is the space of complex structures (or, equivalently, of complete and finite-area hyperbolic structures) on S up to equivalence, where two such structures X and Y are considered as equivalent if there is a conformal map (respectively, an isometry) h : (S, X) → (S, Y ) which is homotopic to the identity map of S.

There are several natural metrics on T(S); some of them are Riemannian and others are only Finsler. Some metrics on Teichmüller space (e.g. the Weil-Petersson metric, and the Teichmüller metric) are strongly related to the complex structure of that space. A Finsler metric is a length metrics where the distance between two points is defined by minimizing lengths of peicewise C 1 paths joining them and where the length of a path is computed by integrating norms of tangent vectors. The fact that a metric is Finsler but not Riemannian means that the norm function on each tangent space is not associated to a scalar product. In this paper, we consider the Teichmüller metric and Thurston's asymmetric metric. Both metrics are Finsler, and in the case of Thurston's asymmetric metric, the Finsler norm not only is non-Riemannian but it is even not symmetric. These metrics are defined in terms of distances between geometric structures on the surface (conformal structures, in the case of the Techmüller metric, and hyperbolic structures, in the case of Thurston's asymmetric metric). The question of comparing the various metrics on T(S) is a natural one, and it was also suggested by Thurston in [START_REF] Thurston | Minimal stretch maps between hyperbolic surfaces[END_REF]; see also [START_REF] Papadopoulos | Shift coordinates, stretch lines and polyhedral structures for Teichmüller space[END_REF].

We shall recall below the definition of Thurston's asymmetric metric. Thurston formulated several ideas concerning this asymmetric metric, and we shall review some of them below. He proved some of the major results and he outlined some other results in his paper [START_REF] Thurston | Minimal stretch maps between hyperbolic surfaces[END_REF]. After Thurston's paper was circulated, questions on the limiting behavior of stretch lines and of anti-stretch lines were considered, see [START_REF] Papadopoulos | On Thurston's boundary of Teichmüller space and the extension of earthquakes[END_REF], [START_REF] Papadopoulos | Shift coordinates, stretch lines and polyhedral structures for Teichmüller space[END_REF], [START_REF] Théret | On the negative convergence of Thurston's stretch lines towards the boundary of Teichmüller space[END_REF] and [START_REF] Théret | On elementary antistretch lines[END_REF]. Stretch lines are some special geodesics of that metric, and anti-stretch lines are stretch lines traversed in the opposite directions; we note that since the metric is not symmetric, a geodesic traversed in the opposite direction is not necessarily a geodesic. There were also generalizations of Thurston's asymmetric metric to the case of surfaces with boundary, see [START_REF] Liu | On length spectrum metrics and weak metrics on Teichmüller spaces of surfaces with boundary[END_REF] and [START_REF] Liu | Length spectra and the Teichmüller metric for surfaces with boundary[END_REF]. A renewal of interest in this metric has emerged recently, see e.g. the papers [START_REF] Walsh | The horoboundary and isometry groups of Thurston's Lipschitz metric[END_REF], [START_REF] Lenzhen | Bounded combinatorics and the Lipschitz metric on Teichmüller space[END_REF] and [START_REF] Liu | On the classification of mapping class actions on Thurston's asymmetric metric[END_REF]. It became also clear that among the known metrics on Teichmüller space, the Thurston metric is the one that has an interesting analogue on outer space, see [START_REF] Culler | Moduli of graphs and automorphisms of free groups[END_REF]. Thus, there are very good reasons to study Thurston's asymmetric metric. Some analogies and some differences between Thurston's asymmetric metric and Teichmüller's metric have described in the paper [START_REF] Papadopoulos | Shift coordinates, stretch lines and polyhedral structures for Teichmüller space[END_REF]. In the present paper, we point out new analogies that concern the Finsler character of these metrics.

Given a set X, a function d : X → X is said to be a weak metric on X if it satisfies all the axioms of a distance function except the symmetry axiom. The weak metric d is said to be asymmetric if it is strictly weak, that is if there exist two points x and y in X such that d(x, y) = d(y, x).

A related notion is that of a weak norm on a vector space E. This is a function p : E → R that is nonnegative, convex and positively homegeneous. In other words, p satisfies the following:

(1) p(x) ≥ 0 for all x ∈ E;

(2) p(λx) = λp(x) for all x ∈ E and for all λ ≥ 0;

(3) p(x 1 + x 2 ) ≤ p(x 1 ) + p(x 2 ) for all x 1 , x 2 ∈ E. In this paper, to simplify terminology, we shall use in general the term norm to denote a weak norm a,d the terms metric to denote a weak metric.

Let X and Y be now two complete hyperbolic metrics on S. In the paper [START_REF] Thurston | Minimal stretch maps between hyperbolic surfaces[END_REF], Thurston defined an asymmetric metric d L on T(S) by setting

(1) d L (X, Y ) = inf f log L f (X, Y ),
where the infimum is taken is over all homeomorphsims f : X → Y homotopic to the indentity map of S, and where

L f (X, Y ) is the Lipschitz constant of f , that is, Lip(f ) = sup x =y∈S d Y f (x), f (y) d X x, y .
We shall call this weak metric Thurston's asymmetric metric or, for short, Thurston's metric.

In the same paper, Thurston proved that there is a (non-necessarily unique) extremal Lipschitz homeomorphsim that realizes the infimum in [START_REF] Ahlfors | Conformal invariants: topics in geometric function theory[END_REF]. He also proved that we have the following formula for the asymmetric metric:

(2) d L (X, Y ) = log sup γ ℓ Y (γ) ℓ X (γ) ,
where ℓ X (γ) denotes the hyperbolic length of γ with respect to the metric X and γ ranges over all essential simple closed curves on S. Furthermore, Thurston proved that the asymmetric metric defined in (1) is Finsler, that is, it is a length metric which is defined by integrating a weak norm on the tangent bundle of T(S) along paths in Teichmüller space, and taking the minimum lengths over all peicewise C 1 -paths. Thurston gave an explicit formula for the weak norm of a tangent vector V at a point X in T(S), namely, (3)

V L = sup λ∈ML dℓ λ (V ) ℓ λ (X) .
Here, ML = ML(S) is the space of measured laminations on S, ℓ λ : T(S) → R is the length function on Teichmüller space associated to the measured lamination λ and dℓ λ is the differential of the function ℓ λ (X) at the point X ∈ T(S).

For X in T(S) and λ in ML, we shall use the notation ℓ λ (X) or ℓ X (λ) to denote the X-length of λ, depending on whether we consider the length function as a function on Teichmüller space or on measured lamination space.

We shall present an analogue of Formula (3) for the Teichmüller metric that is expressed in terms of extremal length, namely, we show that the Finsler weak norm associated to the Teichmüller distance is given by the following formula:

(4) V E = sup λ∈MF dExt 1/2 λ (V ) Ext 1/2 λ (X)
.

Here, X is a conformal structure on S, considered as a point in Teichmüller space, V is a tangent vector at the point X, MF = MF(S) is the space of measured laminations on S, Ext λ : T → R is the extremal length function associated to the measured foliation λ and dExt λ is the differential of that function at X.

There is a more geometric version of (4) in which the tangent vector V is interpreted as a Beltrami differential (see Corollary 3.7 below).

The following table summarizes some analogies between notions and results associated to Thurston's asymmetric metric and those associated to the Teichmüller metric.

(iii) d L (X, Y ) = log inf f L(f ) d T (X, Y ) = log inf f K(f ) (iv) d L (X, Y ) = log sup γ ℓ γ (Y ) ℓ γ (X) d T (X, Y ) = log sup γ Ext 1/2 γ (Y ) Ext 1/2 γ (X) (v) V L = sup λ∈ML dℓ λ (V ) ℓ λ (X) V T = sup λ∈MF dExt 1/2 λ (V ) Ext 1/2 λ (X) (vi) Thurston cataclysm coordinates The homeomorphism X ∈ T(S) → F µ (X) ∈ MF(µ) X ∈ T(S) → F h (Φ F (X)) ∈ MF(F ) (vii) ℓ λ (X) = i(λ, F λ * (X)) Ext λ (X) = i(λ, F h (Φ λ (X))) (viii) dℓ γ (µ) = 2 π Re < Θ α , µ > dExt λ (µ) = -2Re < Φ λ , µ > (ix)
The horofuction boundary is The horofuction boundary is Thurston's boundary Gardiner-Masur's boundary

(x) L X (λ) = ℓ X (λ) L X for λ ∈ ML E X (λ) = Ext 1/2 X (λ) K 1/2 X for λ ∈ ML
Most of the entries in this table are well know, some of them are known but need explanation, and some of the them are new and proved in this paper. We now make a few comments on all the entries.

(i) Stretch maps arise from the extremal problem of finding the best Lipschitz constant of maps homotopic to the identity between two hyperbolic structures on a surface S, in much the same way as Teichmüller maps arise from the extremal problem of finding the best quasiconformal constant of maps homotopic to the identity between two complex structures on S.

(ii) Stretch lines are geodesics for Thurston's metric. A stretch line is determined by a pair (µ, F ) where µ is a complete lamination (not necessarily measured) and F a measured foliation transverse to µ. Thurston proved that any two points in Teichmüller space can be joined by a geodesic which is a concatenation of stretch lines, but in general such a geodesic is not unique. Furthermore, there exist geodesics for Thurston's metric that are not concatenations of stretch lines. This contrasts with Teichmüller's theorem establishing that existence and uniqueness of geodesics joining any two distinct points, cf. [START_REF] Teichmüller | Extremale quasikonforme Abbildungen und quadratische Differentiale[END_REF], [START_REF] Teichmüller | Bestimmung der extremalen quasikonformen Abbildungen bei geschlossenen orientierten Riemannschen Flächen[END_REF] and Ahlfor's survery [START_REF] Ahlfors | Quasiconformal mappings, Teichmüller spaces, and Kleinian groups[END_REF].

(iii) The left hand side is Thurston's definition of Thurston's metric, the infimum is over all homeomorphisms f homotopic to the identity and L(f ) is the Lispschitz constant of such a homeomorphism. The right hand side is the definition of the Teichmüller metric, the infimum is over all quasiconformal homeomorphisms f homotopic to the identity and K(f ) is the dilatation of such a homeomorphism.

(iv) The left hand side is another expression (also due to Thurston) of Thurston's metric, and the right hand side is Kerckhoff's formula for the Teichmüller metric.

(v) Here, V is a tangent vector to Teichmüller space at a point X. The left hand side formula is due to Thurston [START_REF] Thurston | Minimal stretch maps between hyperbolic surfaces[END_REF]; it is the infinitesimal form of Thurston's metric. The right hand side is an infinitesimal form of the Teichmüller metric, and it is proved below (Theorem 3.6).

(vi) On the left hand side, µ is a complete geodesic lamination and the range of the map, MF(µ), is the subspace of MF of equivalence classes of measured foliations that are transverse to µ. On the right hand side, F is a measured foliation, and the range of the map, MF(F ), is the subspace of MF consisting of equivalence classes of measured foliations that are transverse to F . The left hand side map is a homeomorphism defined in Thurston's paper [START_REF] Thurston | Minimal stretch maps between hyperbolic surfaces[END_REF]. The right hand side map is a homeomorphism that arises from the fact that a pair F 1 , F 2 of measured foliations determines a unique point X ∈ T(S) and a unique quadratic differential Φ on X such that F 1 ad F 2 are measure equivalent to the vertical and horizontal foliations of Φ. Here, F h (Φ F (X)) is the (equivalence class of the) horizontal foliation of the quadratic differential on the Riemann surface X having F as vertical foliation.

(vii) The left hand side is an expression of the length of a complete lamination λ as the geometric intersection with the horocyclic foliation associated to a completion λ * of λ. The formula is proved in [START_REF] Papadopoulos | On Thurston's boundary of Teichmüller space and the extension of earthquakes[END_REF] for the case where λ is complete. The case where λ is not complete follows easily from the geometric arguments used in that proof. The right hand side formula is due to Kerckhoff, see Lemma 3.3 below.

(viii) In the left hand side formula, µ is a Beltrami differential. The formula was given by Gardiner [START_REF] Gardiner | Schiffer's interior variation and quasiconformal mapping[END_REF], and in this formula α is (the homotopy class of) a simple closed curve and Θ α is the Poincaré series of α, that is, a quadratic differential on the hyperbolic surface X defined by

Θ α = B∈<A>\Γ B * ( dz z ) 2
where Γ ⊂ PSL(2, R) is a Fuchsian group associated to X acting on the upper halfplane H 2 and A(z) = e ℓα z is the deck transformation corresponding to the simple closed geodesic α. The right hand side formula is called Gardiner's extremal length variational formula (see Lemma 3.2 below).

(ix) Walsh showed in [START_REF] Walsh | The horoboundary and isometry groups of Thurston's Lipschitz metric[END_REF] that the horofuction boundary of Thurston's metric is canonically identified with Thurston's boundary. Liu and Su showed in [START_REF] Liu | The horofunction compactification of Teichmüller space[END_REF] that the horofuction boundary of the Teichmüller metric is canonically identified with Gardiner-Masur's boundary.

(x) Here, we choose a basepoint X 0 ∈ T(S) and we fix a complete geodesic lamination µ on S. For each X ∈ T(S), L X is the Lipschitz constant of the extremal Lipschitz map between X 0 and X and K X is the quasiconformal dilatation of the extremal quasiconformal map between X 0 and X.

We illustrate the use of the functions L X and E X . We use Thurston's homeomorphism φ µ : T(S) → M F (µ) that maps each X ∈ T(S) to the horocylic foliation F µ (X). The definition of this homeomorphism is recalled in §2 below. We denote the projective class of a measured foliation F by [F ]. By a result in [START_REF] Papadopoulos | On Thurston's boundary of Teichmüller space and the extension of earthquakes[END_REF], a sequence (X k ) k≥1 in T(S) coverges to a limit [F ] ∈ PMF if and only if F µ (X k ) tends to infinity and [F µ (X k )] coverges to [F ]. We prove the following: [START_REF] Fathi | Travaux de Thurston sur les surfaces[END_REF], there exists a sequence (c k ) k≥1 of positive numbers such that c k → 0 and ( 5)

Proposition 1.1. If a sequence (X k ) k≥1 in T(S) coverges to a limit [F ] ∈ PMF, then L X k (•) converges to [F ] in the following sense: (⋆) up to a subsequence, L X k (•) converges to a positive multiple of i(F, •) uniformly on any compact subset of ML. Proof. Assume that a sequence (X k ) k≥1 coverges to [F ] ∈ PMF. Let L X k = L k and F µ (X k ) = F k . From Thurston's theory
c k i(F k , λ) → i(F, λ)
for each λ ∈ ML.

By the Fundamental Lemma in [START_REF] Papadopoulos | On Thurston's boundary of Teichmüller space and the extension of earthquakes[END_REF], there is a uniform constant C > 0 such that

(6) i(F k , λ) ≤ ℓ λ (X k ) ≤ i(F k , λ) + C. Note that c k ℓ λ (X k ) ≤ c k L k ℓ λ (X 0 )
and for each L k , there is a measured lamination

λ k with ℓ λ k (X 0 ) = 1 such that L k = ℓ λ k (X k ).
As a result, we have

(7) c k ℓ λ (X k ) ℓ λ (X 0 ) ≤ c k L k ≤ c k ℓ λ k (X k ).
From ( 5), ( 6), [START_REF] Gardiner | Schiffer's interior variation and quasiconformal mapping[END_REF], it follows that c k L k is uniformly bounded from above and uniformly bounded below away from zero. It follows that, up to a subsequence, ℓ λ (X k )

L k converges to a positive multiple of i(F, λ) for each λ ∈ ML. Since pointed convergence of ℓ λ (X k ) L k , λ ∈ ML is equivalent to uniform convergence on compact subsets of ML, property (⋆) holds.

Concerning the right hand side formula in (x), a similar result was obtaind by Miyachi in [START_REF] Miyachi | Teichmüller rays and the Gardiner-Masur boundary of Teichmüller space[END_REF] : Proposition 1.2. If a sequence (X k ) k≥1 in T(S) coverges to a limit P in Gardiner-Masur's boundary, then E X k (•) converges to some function E P (•) in the following sense: Up to a subsequence, E X k (•) converges to a positive multiple of E P (•) uniformly on any compact subset of ML.

Lipschitz Norm

This section contains results of Thurston from his paper [START_REF] Thurston | Minimal stretch maps between hyperbolic surfaces[END_REF] that we will use later in this paper. We have provided proofs because at times Thurston's proofs in [START_REF] Thurston | Minimal stretch maps between hyperbolic surfaces[END_REF] are considered as sketchy.

A geodesic lamination µ on a hyperbolic surface X is said to be complete if its complementary regions are all isometric to ideal triangles. (We note that we are dealing with laminations µ that are not necessarily measured, except if specified.) Associated with (X, µ) is a measured foliation F µ (X), called the horocyclic foliation, satisfying the following three properties:

(i) F µ (X) intersects µ transversely, and in each cusp of an ideal triangle in the complement of µ, the leaves of the foliation are pieces of horocycles that make right angles with the boundary of the triangle; (ii) on the leaves of µ, the transverse measure for F µ (X) agrees with arclength; (iii) there is a nonfoliated region at the centre of each ideal triangle of S \ µ whose boundary consists of three pieces of horocycles that are pairwise tangent (see Figure 1). We denote by MF(µ) the space of measured foliations that are transverse to µ.

Thurston [START_REF] Thurston | Minimal stretch maps between hyperbolic surfaces[END_REF] proved the following fundamental result.

Theorem 2.1. The map φ µ : T(S) → MF(µ) defined by X → F µ (X) is a homeomorphism.

The stretch line directed by µ and passing through X ∈ T(S) is the curve

R ∋ t → X t = φ -1 µ (e t F µ (X)
). We call a segment of a stretch line a stretch path.

Suppose that µ is the support of a measured geodesic lamination λ. Then, for any two points X s , X t , s ≤ t on the stretch line, their Lipschitz distance d L (X s , X t ) is equal to ts, and this distance is realized by log ℓ λ (X t ) ℓ λ (X s ) .

We denote by ML the space of measured geodesic laminations on X and we let ML 1 = {λ ∈ ML | ℓ λ (X) = 1}. We may identify ML 1 with PL, the space of projective measured laminations.

Thurston [START_REF] Thurston | Minimal stretch maps between hyperbolic surfaces[END_REF] introduced a Finsler structure on T(S) by defining the Finsler norm of a tangent vector V ∈ T X T(S) by the following formula :

(8) V L = sup λ∈ML dℓ λ (V ) ℓ λ (X) .
Note that we may write

V L = sup λ∈ML1 dℓ λ (V ).
Lemma 2.2. Let Γ µ (t) be a stretch line through X with Γ µ (0) = X and Γµ (0) = V 0 . For any measured lamination λ supported by µ (λ may be not unique), we have

V 0 L = dℓ λ (V 0 ) ℓ λ (X) .
Proof. By multiplying the transverse measure by a constant, we may assume that ℓ λ (X) = 1. The stretch map from X to Γ λ (t) is a Lipschitz map with Lipschitz constant e t , and the hyperbolic length of λ in Γ λ (t) is equal to e t . As a result, for any measured lamination γ ∈ ML 1 ,

ℓ γ (Γ λ (t)) ≤ e t .
It follows that sup

γ∈ML1 dℓ γ (V 0 ) = dℓ λ0 (V 0 ).
Theorem 2.3 (Thurston [START_REF] Thurston | Minimal stretch maps between hyperbolic surfaces[END_REF] p. 20). • L is the infinitesimal norm of Thurston's asymmetric distance d L .

Proof. Recall that the Finsler norm • L induces an (asymmetric) distance:

d(X, Y ) = inf Γ 1 0 Γ L ,
where the infimum is taken over all piecewise C 1 curves Γ : [0, 1] → T in Teichmüller space joining X to Y . By compactness of PML, there exists an element λ ∈ ML such that the distance d L (X, Y ) is attained at λ, that is,

d L (X, Y ) = log ℓ λ (Y ) ℓ λ (X) .
For any piecewise C 1 path Γ : [0, 1] → T(S) satisfying Γ(0) = X and Γ(1) = Y , we have

d L (X, Y ) = log ℓ λ (Y ) ℓ λ (X) = 1 0 d log ℓ λ ( Γ(t)) ≤ 1 0 Γ(t) L . Therefore, d L (X, Y ) ≤ d(X, Y ).
To prove the reverse inequality, we use a result proved by Thurston in [START_REF] Thurston | Minimal stretch maps between hyperbolic surfaces[END_REF], namely, that any two points in Teichmüller space can be connected by a finite concatenation of stretch paths, each of which stretches along some common measured geodesic lamination λ.

As a result, for any distinct X and Y in T(S), we can assume that there exists a path Γ connecting X and Y , with Γ being a concatenation Γ 1 * • • • * Γ n of stretch paths. Up to reparametrization, we may assume that each Γ i is defined on [0, 1]. Furthermore, by a result of Thurston, we may assume that each stretch path Γ i is directed by a complete geodesic lamination µ i that contains a common measured lamination λ, which is the maximally stretched lamintation from X to Y . This is a consequence of Theorem 8.2 in Thurston [START_REF] Thurston | Minimal stretch maps between hyperbolic surfaces[END_REF], in which Thurston shows that there is a unique maximal ratio-maximizing chain recurrent lamination which contains all other ratio-maximizing chain recurrent lamitations. For our purposes, we take λ to be the maximal (with respect to inclusion) measured lamination contained in µ(X, Y ).

Along each Γ i , by Lemma 2.2, we have

sup ν∈ML dℓ ν ( Γi ) ℓ ν = dℓ λ ( Γi ) ℓ λ .

Now we have

Γ Γ L = n i=1 Γi Γ L = n i=1 Γi dℓ λ ( Γi ) ℓ λ = n i=1 log ℓ λ (Γ i (1)) ℓ λ (Γ i (0)) = log ℓ λ (Y ) ℓ λ (X) .
For the last equality, note that we have Γ i (1) = Γ i+1 (0) for each i and Γ

1 (0) = X, Γ n (1) = Y . It follows that d(X, Y ) ≤ d L (X, Y ).

Teichmüller Norm

Extremal length is an tool in study of the Teichmüller metric. The notion is due to Ahlfors and Beurling [START_REF] Ahlfors | Conformal invariants and function-theoretic null-sets[END_REF], see also [START_REF] Ahlfors | Conformal invariants: topics in geometric function theory[END_REF]. We briefly recall the definition. Given a Riemann surface X, a conformal metric σ on X is a metric that is locally of the form σ(z)|dz| where z is a local holomorphic parameter and σ(z) ≥ 0 is a Borel measurable function in the local chart. We define the σ-area of X by

A(σ) = X σ 2 (z)|dz| 2 .
Given a homotopy class of simple closed curves α, its σ-length is defined by

L σ (α) = inf α ′ α ′ σ(z)|dz|,
where the infimum is taken over all essential simple closed curves α ′ in the homotopy class α. (Note that it follows from the invariance property of the expression σ(z)|dz| that the two integrals A(σ) and L σ (α) are well defined, that is, they can be computed in the holomorphic local coordinates and the value obtained does not depend on the choice of the coordinates, see [START_REF] Ahlfors | Conformal invariants: topics in geometric function theory[END_REF].)

With the above notation, we define the extremal length of α on X by

Ext α (X) = sup σ L 2 σ (α) A(σ) ,
where σ(z)|dz| ranges over all conformal metrics on X satisfying 0 < A(σ) < ∞. By a result of Kerckhoff [START_REF] Kerckhoff | The asymptotic geometry of Teichmüller space[END_REF], there is a unique continuous extension of the extremal length function to the space of measured foliations MF, with Ext aγ (X) = a 2 Ext γ (X) for a > 0 and γ a homotopy class of simple closed curves on X. Recall that a tangent vector to Teichmüller space at a point X can be represented by a Beltrami differential µ = µ(z) dz dz [START_REF] Imayoshi | An Introduction to Teichmüller Spaces[END_REF]. (Note that in this section we use the letter µ to denote a Beltrami differential since this is the traditional notation, although in other sections the same letter is used to denote a lamination. Hopefully, there will be no confusion.) Using the extremal length function, we define a norm on the tangent space at X by setting: ( 9)

µ E = sup λ∈MF dExt 1/2 λ (µ) Ext 1/2 λ (X)
.

The following theorem is due to Hubbard-Masur [START_REF] Hubbard | Quadratic differentials and foliations[END_REF]; we refer to Kerckhoff [START_REF] Kerckhoff | The asymptotic geometry of Teichmüller space[END_REF] for a short proof. Theorem 3.1. For any Riemann surface X and for any foliation on X, there is exactly one quadratic differential, denoted by Φ λ , whose vertical measured foliation is measure-equivalent to λ. Lemma 3.2 (Gardiner [START_REF] Gardiner | Measured foliations and the minimal norm property for quadratic differentials[END_REF]). The extremal length function Ext λ is differentiable and we have the following formula, called the "first variational formula":

(10) dExt λ (µ) = -2Re < Φ λ , µ >,
where < Φ λ , µ > is the natural pairing

< Φ λ , µ >= X Φ λ (z)µ(z)dxdy.
The following observation is due to Kerckhoff [START_REF] Kerckhoff | The asymptotic geometry of Teichmüller space[END_REF]:

Lemma 3.3. With the above notation, (11) Ext λ
(X) = X |Φ λ (z)|dxdy.
Proof. The proof we give here is due to Ivanov [START_REF] Ivanov | Isometries of Teichmüller spaces from the point of view of Mostow rigidity[END_REF]. We include it for completeness. By continuity and the density of weighted simple closed curve in MF, it suffices to prove [START_REF] Imayoshi | An Introduction to Teichmüller Spaces[END_REF] for the case where µ = aγ ∈ MF, where γ is (the homotopy class of) a simple closed curve and a > 0.

Let Φ be the one-cylinder Strebel differential on X determined by aγ. The complement of the vertical critical leaves of Φ is a cylinder foliated by circles isotopic to γ. Let us also set ρ = |Φ| 1/2 |dz|. Then ρ is a flat metric on S, with a finite number of singular points, which are conical singularities. Measured in the flat metric ρ, the circumference and height of the cylinder are equal to L ρ (γ) and a respectively. By a theorem of Jenkins-Strebel [START_REF] Strebel | Quadratic differentials[END_REF], the extremal length Ext X (γ) of γ is equal to

Ext γ (X) = L ρ (γ) a , where ρ = |Φ| 1/2 |dz|. The area A(ρ) = X |Φ λ (z)|dxdy of the cylinder is equal to aL ρ (γ). As a result, Ext aγ (X) = a 2 Ext γ (X) = aL ρ (γ) = A(ρ).
The next result follows from [START_REF] Ivanov | Isometries of Teichmüller spaces from the point of view of Mostow rigidity[END_REF] and [START_REF] Imayoshi | An Introduction to Teichmüller Spaces[END_REF].

Proposition 3.4. The norm µ E satisfies

(12) µ E = sup λ∈MF -Re < Φ λ , µ > X |Φ λ (z)|dxdy = sup Φ =1 Re < Φ, µ >,
where Φ varies over all holomorphic quadratic differentials Φ.

Note that we have Φ = X |Φ(z)|dxdy.

Now we consider the Teichmüller metric on T(S). We recall that it is defined by

d T (X, Y ) := 1 2 inf f log K(f )
where f : X → Y is a quasi-conformal map homotopic to the identity map of S and

K(f ) = x∈X K x (f ) 1
is the quasi-conformal dilatation of f (the sup here denotes essential supremum), with

K x (f ) = |f z (x)| + |f z (x)| |f z (x)| -|f z (x)
| being the pointwise quasiconformal dilatation at the point x ∈ X with local conformal coordinate z.

Teichmüller's theorem states that given any X, Y ∈ T(S), there exists a unique quasi-conformal map f : X → Y , called the Teichmüller map, such that

d T (X, Y ) = 1 2 log K(f ). The Beltrami coefficient µ := ∂f ∂f is of the form µ = k Φ |Φ|
for some quadratic differential Φ on X and some constant k with 0 ≤ k < 1. In some natural coordinates given by Φ on X and for some associated quadratic differential Φ ′ on Y , the Te-

ichmüller map f is given by f (x+iy) = K 1/2 x+iK -1/2 y, where K = K(f ) = 1 + k 1 -k .
It is known that the Teichmüller metric is a Finsler metric and that between any two points in T(S) there is exactly one geodesic. A geodesic ray with initial point X is given by the one-parameter family of Riemann surfaces {X t } t≥0 , where there is a holomorphic quadratic differential Φ on X and a family of Teichmüller maps f t : X → X t , with initial Beltrami differential µ(f t ) = e 2t -1 e 2t + 1 Φ |Φ| . Here µ(f t ) is chosen such that the geodesic ray has unit speed, that is, d T (X s , X t ) = ts for all s ≤ t. We also recall the following formula due to Kerckhoff [START_REF] Kerckhoff | The asymptotic geometry of Teichmüller space[END_REF].

Theorem 3.5. Let X, Y be any two points in T(S). Then

d T (X, Y ) = 1 2 log sup λ Ext λ (Y ) Ext λ (X) ,
where λ ranges over elements in MF.

Now we can prove the following:

Theorem 3.6. The metric induced by the norm • E defined in ( 9) is the Teichmüller metric.

Proof. Denote by d E the length metric on T(S) induced by the norm µ E . For any X, Y in T(S), by Kerckhoff's formula, the Teichmüller distance is realized by

d T (X, Y ) = 1 2 log Ext λ (Y ) Ext λ (X)
for some measured foliation λ. An argument similar to the one in the first part of the proof of Theorem 2.3 shows that d T (X, Y ) ≤ d(X, Y ). For the converse, let Γ(t) be a Teichmüller geodesic connecting X and Y . We parametrize Γ with unit speed, and define it on the interval [0, T ] with T = d T (X, Y ) with Γ(0) = X and Γ(T ) = Y .

For each t in [0, T ], the Teichmüller map between X = Γ(0) and Γ(t) has quasiconformal dilatation e 2t . It follows from the geometric definition of a quasiconformal map that for any λ in MF, we have Ext λ (Γ(t)) ≤ e 2t Ext λ (Γ(0)).

Taking square roots, we get Ext

1/2 λ (Γ(t)) ≤ t 1/2 λ (Γ(0)) and then lim t→0 + Ext 1/2 λ (Γ(t)) -Ext 1/2 λ (Γ(0)) tExt 1/2 λ (Γ(t)) ≤ lim t→0 + e t Ext 1/2 λ (Γ(0)) -Ext 1/2 λ (Γ(0)) tExt 1/2 λ (Γ(t)) = lim t→0 + e t -1 t = 1.
As a result, d log Ext

d log Ext 1/2 λ ( Γ(t)) ≤ 1 for all t ∈ R.
Considering the integral of the norm • E along Γ(t), we have

T 0 Γ(t) E ≤ T = d T (X, Y ).
This proves that d(X, Y ) ≤ d T (X, Y ).

Combining Proposition 3.4 and Theorem 3.6 , we have the following corollary.

Corollary 3.7. The Teichmüller Finsler norm, denoted by µ T , is given by

µ T = sup Φ =1 Re < Φ, µ >,
where Φ varies over all holomorphic quadratic differentials Φ.

The above result was already known and can be obtained by using the famous Reich-Strebel inequality [START_REF] Reich | Extremal quasiconformal mappings with given boundary values, Contributions to Analysis[END_REF]. The result means that the Teichmüller norm

• T is dual to the L 1 -norm Φ = X |Φ(z)|dxdy
on the cotangent space.

Convex embedding of measure foliation space using extremal length

Following a usual trend, we shall call the boundary of a convex body a convex sphere. Note that the boundary of a convex body is always homeomorphic to a sphere.

Thurston [START_REF] Thurston | Minimal stretch maps between hyperbolic surfaces[END_REF] proved that for any X ∈ T(S), the function dℓ : ML 1 → T * X T(S) defined by λ → dℓ λ embeds ML 1 as a convex sphere in the cotangent space of Teichmüller space at X containing the origin. This embedding is the dual of the boundary of the unit ball in the tangent space representing vectors of norm one for the Finsler structure associated to Thurston's metric. We prove an analogous result for the extremal length function.

4.1. Given a point X in T(S), let PF = MF 1 = {λ ∈ MF | Ext λ (X) = 1}.
Then, the function PF → T * X T(S) defined by λ → dExt λ embeds PF as a convex sphere in T * X T(S) containing the origin. This sphere is the boundary of the dual of the ball the tangent space representing vectors of norm one for the Finsler structure associated to the Teichmüller metric.

Proof. We fix a point X in T(S). By Gardiner's formula, dExt λ (µ) = -2Re < Φ λ , µ > (see Theorem 3.1 for the notation). Kerckhoff proved that if a sequence λ n in PF converges to λ then Φ λn -Φ λ → 0 (see Page 34-35 of Kerckhoff's paper [START_REF] Kerckhoff | The asymptotic geometry of Teichmüller space[END_REF], where he proved the result for the case where each λ n is a simple closed curve. Since the subset of weighted simple closed curves is dense in MF, the result we need follows). It follows that the map dExt is continuous.

Moreover, since the pairing between quadratic differentials and harmonic Beltrami differentials is nondegenerate [START_REF] Imayoshi | An Introduction to Teichmüller Spaces[END_REF], if dExt λ = dExt λ ′ , then Φ λ = Φ λ ′ . This means that λ = λ ′ . Therefore, the map dExt is injective.

Let dExt(PF) be the image of PF under the map dExt, that is,

dExt(PF) = {dExt λ ∈ T * X T(S) | λ ∈ PF}.
Since PF is homeomorphic to a sphere, by invariance of domain, dExt(PF) is open and dExt is a homeomorphism between PF and dExt(PF).

For any λ ∈ PF, consider the Beltrami differential

V λ = - Φ λ 2|Φ λ | .
It is a tangent vector in T X T(S). By Gardiner's formula,

dExt λ (V λ ) = -2Re < Φ λ , - Φ λ 2|Φ λ | > = X |Φ λ |dxdy = Ext λ (X) = 1.
As a result, the derivative dExt λ in the direction V λ is 1. Using Hölder's inequality, for any γ ∈ PF, we have

dExt γ (V λ ) = Φ γ Φ λ |Φ λ | dxdy ≤ X |Φ γ |dxdy = Ext γ (X) = 1.
Moreover, equality holds if and only if γ = λ. Therefore, the derivative dExt γ (V λ ) of any other γ ∈ PF is strictly less than 1. Denote by C(dExt(PF)) the convex hull of dExt(PF) in T * X T(S). Thus, V λ defines a non-constant linear functional on the cotangent space which attains its maximal value on C(dExt(PF)) at dExt λ . As a result, dExt λ is an extreme point of the convex set C(dExt(PF)).

We now use Thurston's argument in his proof of Theorem 5.1 of [START_REF] Thurston | Minimal stretch maps between hyperbolic surfaces[END_REF]. The set of extreme points of any convex set is a sphere of some dimension. Since the dimension of dExt(PF) is one less than the dimension of T * X T(S), the set of extremal points of C(dExt(PF)) coincides with dExt(PF). As a result, dExt(PF) is a convex sphere. To see that the convex set C(dExt(PF)) contains the origin in its interior, note that since this convex hull has a nonempty interior, there must be at least one line through the origin of T * X T(S) which intersects dExt(PF) in at least two points. For each of these points, there is a linear functional which attains its positive maximum value there. It follows that 0 must separate the two points, so 0 is a convex combination of them.

Comparison between the Lipschitz and Teichmüller norms

From an inequality called Wolpert's inequality [START_REF] Wolpert | The length spectra as moduli for compact Riemann surfaces[END_REF], we have d L ≤ 2d T . (In fact, the inequality is contained in Sorvali's paper [START_REF] Sorvali | On the dilatation of isomorphisms between covering groups[END_REF], and it was rediscovered by Wolpert. ) Choi and Rafi [START_REF] Choi | Comparison between Teichmüller and Lipschitz metrics[END_REF] proved that the two metrics are quasi-isometric in any thick part of the Teichmüller space, but that are (X n ), (Y n ) in the thin part, with d L (X n , Y n ) → 0.

It is interesting to compare the Lipschitz norm and the Teichmüller norm. In particular, we ask the following: Question 5.1. Determine a function C(ǫ) such that

V L ≤ V T ≤ C(ǫ) V L
for any V ∈ T X T(S) and for all X ∈ T ǫ (S), the ǫ-thick part of T(S).

Since V L (respectively V T ) is given by the logarithmic derivative of the hyperbolic (respectively extremal) length function of measured laminations, a further study of the relation between hyperbolic and extremal length may give the answer to the above question. Question 5.2. Consider the length-spectrum metric d ls on T(S). This can be defined as a symmetrization of Thurston's metric 1 , by the formula

d ls (X, Y ) = max{log sup γ ℓ Y (γ) ℓ X (γ) , log sup γ ℓ X (γ) ℓ Y (γ) } = max{d L (X, Y ), d L (Y, X)}.
Is the length-spectrum metric d ls a Finsler metric ? If yes, give a formula for the infinitesimal norm of a vector on Teichmüller space with respect to this Finsler structure.

There are comparisons between the length-spectrum metric and the Teichmüller metric. It was shown by Liu and Su [START_REF] Liu | Almost-isometry between Teichmüller metric and length-spectra metric on moduli space[END_REF] that the divergence of d ls and d T is only caused by the action of the mapping class group. In fact, they showed that d ls and d T are "almost" isometric on the moduli space and that the two metrics on the moduli space determine the same asymptotic cone. Question 5.3. One can also seek for results analogous to those presented here for the weak metric L * dual to Thurston's metric, that is, the weak metric on Teichmüller space defined by L * (X, Y ) = L(Y, X), as well as for the length spectrum metric d ls on the same space.

Figure 1 .

 1 Figure 1. The horocyclic foliation of an ideal triangle.

1 / 2 λ

 12 ( Γ(0)) ≤ 1 and then

Theorem


Anaïs Crestetto

Philippe Helluy

Multi-Water-Bag Model And Method Of Moments For The Vlasov Equation

Keywords: Vlasov-Poisson, water-bag approximation, Burgers equation, multivalued solution MSC2010: 35Q83, 44A60, 65M08

The kinetic Vlasov-Poisson model is very expensive to solve numerically. It can be approximated by a multi-water-bag model in order to reduce the complexity. This model amounts to solve a set of Burgers equations, which can be done easily by finite volume methods. However, the solution is naturally multivalued (filamentation). The multivalued solution can be computed by a moment method. We present here several numerical experiments.

Introduction

Kinetic equations are used in several domains, such as plasma physics or bubble flows in gases or liquids. The distribution function depends on space and time but also on an additional velocity variable. Computations are thus very expensive. It is of great interest to reduce the complexity of the resolution by using fluid methods.

We consider a plasma containing ions of positive charge and electrons of negative charge. Ions are much heavier than electrons so that we can neglect their displacement and assume that their density n 0 is constant. The electrons move following the system of Vlasov-Poisson in a periodic domain in x:

∂ t f (x, v, t) + v • ∇ x f (x, v, t) + q m E (x, t) • ∇ v f (x, v, t) = 0, ( 1
) div E (x, t) = ρ (x, t) = q m f (x, v, t) dv -n 0 , (2)
E (x, t) dx = 0, (3)
f (x, v, 0) = f 0 (x, v) ≈ n 0 . (4)
The unknowns are the distribution function of electrons f , and the electric field E. The electric field depends on space and time x, t, while the distribution function depends also on an additional velocity variable v. The charge and the mass of one electron are noted q < 0 and m > 0 respectively. Without loss of generality, we can take q m = -1 and n 0 = 1. In one dimension of space, this system becomes:

∂ t f (x, v, t) + v∂ x f (x, v, t) -E (x, t) ∂ v f (x, v, t) = 0, (5)
∂ x E (x, t) = ρ (x, t) = 1 -f (x, v, t) dv, (6)
E (x, t) dx = 0, (7)
f (x, v, 0) = f 0 (x, v) ≈ 1. (8) * Calvi, INRIA Nancy -Grand Est & IRMA, Université de Strasbourg † IRMA, Université de Strasbourg 1
The distribution function f is initially a perturbation of the equilibrium n 0 . After simple calculations, Equation (6) can also be written

∂ t E (x, t) = vf (x, v, t) dv -vf (0, v, t) dv. (9)
In higher dimensions, Equation ( 6) would be replaced by a Poisson equation, assuming that E = -∇Φ, where Φ is the electric potential. The solution can be stable (for example in the case of Landau damping). But since there is no dissipation, the solution can become unstable and filamentation can appear.

In order to reduce the complexity of our model, we propose to approximate the kinetic equation by fluid models. We consider two possibilities:

• the multi-water-bag model, • the method of moments.

Then we propose numerical approximations of these models by simple finite volume schemes. The numerical results will be compared to those obtained with a full resolution of the kinetic model by the popular Particle-In-Cell (PIC) method (described for example in [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF]).

Multi-water-bag model

The multi-water-bag model, detailed in [START_REF] Besse | The multi-water-bag equations for collisionless kinetic modeling[END_REF], generalizes the water-bag model (see for example [START_REF] Bertrand | Non linear electron plasma oscillation: the water bag model[END_REF]). It consists of replacing f by a piecewise constant approximation in the velocity variable. Each piece is called a "water-bag". It is possible to compute only the boundaries of the water-bags.

Presentation of the model

Let N be an integer, v + j (x, t) and v - j (x, t) velocities, A j constants, j = 1, . . . , N , such that we can write:

f (x, v, t) = N j=1 A j H v + j (x, t) -v -H v - j (x, t) -v , ( 10
)
where H is the Heaviside function:

H (u) = 0 if u < 0, 1 if u > 0.
Injecting this expression in the Vlasov equation, and assuming that there is no two v ± j (x, t) equal, we obtain the following system:

∂ t v ± j + v ± j ∂ x v ± j + E = 0, ∀j = 1, . . . , N, (11)
∂ t E = - N j=1 A j v + 2 j (x, 0) -v -2 j (x, 0) 2 - v + 2 j (x, t) -v -2 j (x, t) 2 . ( 12
)
The velocities follow a Burgers equation with a source term. Instead of evolving the distribution function f , we evolve these velocities. The natural solution can become multivalued (filamentation). The weak entropy solution is only an approximation in this context, when shocks appear.

Numerical scheme

We discretize our domain: x i = x 0 + i∆x, with ∆x being the spacial step, and consider a time step ∆t such that t n = n∆t. We evolve each velocity independently by using, for example, the Godunov scheme:

v ±,n+1 j,i = v ±,n j,i - ∆t ∆x    v ±,⋆ j,i+ 1 2 2 2 - v ±,⋆ j,i-1 2 2 2    -∆tE n i , (13) 2
where v ±,n j,i ≃ v ± j (x i , t n ) and E n i ≃ E (x i , t n ). We compute the v j ±,⋆ with an exact Riemann solver.

We compute the electric field with the following scheme:

E n+1 i = E n i -∆t N j=1 A j    v +,0 j,i 2 -v -,0 j,i 2 2 - v +,n j,i 2 -v -,n j,i 2 2    . (14)
In higher dimension, this step should be replaced by the numerical resolution of a Poisson equation.

Remarks

Before the shock, the weak solution and the multivalued solution coincide. After the shock, the natural solution is multivalued, filaments or branches appear, and the weak solution becomes discontinuous. In Section 4, we compare numerically the two kinds of solutions. It is also possible to compute several branches by the method of moments [START_REF] Brenier | A kinetic formulation for multi-branch entropy solutions of scalar conservation laws[END_REF].

Method of moments

The method of moments, presented in [START_REF] Desjardins | A quadrature-based moment method for dilute fluid-particle flows[END_REF][START_REF] Fox | Numerical simulation of spray coalescence in an Eulerian framework: direct quadrature method of moments and multi-fluid method[END_REF], consists in taking the first moments of the equation that we have to solve in order to reduce the number of variables. The system is closed by assuming that the distribution function is made of water-bags.

Presentation of the method

Definition 1: The moment M k of order k of f is defined by:

M k (x, t) = +∞ -∞ v k f (x, v, t) dv. (15)
Taking the 2N first moments of the Vlasov equation:

v k ∂ t f (x, v, t) dv + v k+1 ∂ x f (x, v, t) dv -E (x, t) v k ∂ v f (x, v, t) dv = 0, (16)
we obtain the following system of moments, coupled with the equation for the electric field:

∂ t M 0 + ∂ x M 1 = 0 (17) ∂ t M k + ∂ x M k+1 + kEM k-1 = 0, for k = 1, . . . , 2N -1, ( 18
)
∂ t E (x, t) = M 1 (x, t) -M 1 (0, t) . ( 19
)
We now have a system of 2N + 1 equations, in which the velocity no longer appears. There are 2N + 2 unknowns: the moments M k for k = 0, . . . , 2N and E. We have to close this system by finding an expression for M 2N .

Closure relation

We represent f by water-bags for closing the system:

f (x, v, t) = N j=1 A j H v + j (x, t) -v -H v - j (x, t) -v . ( 20
)
We obtain:

M k (x, t) = N j=1 A j v + k+1 j (x, t) -v -k+1 j (x, t) k + 1 , ∀ k = 0, . . . , 2N, (21)
and thus have an expression of M 2N , assuming that we know the v ± j .

Numerical scheme

At time t n , we know the 2N first moments of f and the electric field. We solve the system:

M k (x, t) = N j=1 A j v + k+1 j (x, t) -v -k+1 j (x, t) k + 1 , ∀ k = 0, . . . , 2N -1, (22)
to obtain v ± j , at time t n . In general, the system may have several solutions. Uniqueness can be recovered through an entropy argument [START_REF] Brenier | A kinetic formulation for multi-branch entropy solutions of scalar conservation laws[END_REF].

To solve this system, we can use the Newton method. But we have no rigorous result of convergence.

When the problem is such that f can be written:

f (x, v, t) = 2N j=1 (-1) j-1 H a + j (x, t) -v , (23)
where -∞ < a 2N ≤ • • • ≤ a 1 < +∞, we can use the algorithm described in [START_REF] Brenier | A kinetic formulation for multi-branch entropy solutions of scalar conservation laws[END_REF][START_REF] Gosse | Resolution of the finite Markov moment problem[END_REF], which solves rigorously such a system. It is then possible to catch numerically the multivalued solutions. For approximating the system of moments, we use a natural kinetic scheme. Formally, we write :

∂ t f (x, v, t) + v + ∂ x f (x, v, t) + v -∂ x f (x, v, t) -E (x, t) ∂ v f (x, v, t) = 0, (24)
where v + = max (0, v) and v -= min (0, v). Denoting ∆x as the space step, ∆t as the time step, and f n i as the approximation of f (x i , v, t n ), we use the following upwind discretization for the Vlasov equation:

f n+1 i -f n i ∆t + 1 ∆x v + f n i -f n i-1 + v -f n i+1 -f n i -E n i ∂ v f n i = 0. (25)
We multiply it by v k and integrate it in v:

1 ∆t v k f n+1 i -f n i dv + 1 ∆x v + k+1 f n i -f n i-1 + v -k+1 f n i+1 -f n i dv -E n i v k ∂ v f n i dv = 0. ( 26
)
We obtain a finite volume scheme:

M n+1 k,i -M n k,i ∆t + 1 ∆x F M n k,i , M n k,i+1 -F M n k,i-1 , M n k,i -kE n i M n k-1,i = 0. ( 27
)
The scheme for the electric field is:

E n+1 i -E n i ∆t = M n 1,i -M n 1,0 . (28)

Numerical results

We validate our models on classical test cases: Landau damping and two stream instability, and obtain good decrease rates for the electric potential energy, when N is big enough. An example for a Landau damping test case is given in Figure 1, with N = 5.

We are now interested in solutions that can initially exactly be depicted by the multi-waterbag model, and that are unstable. We compare the three methods in Figure 2, for N = 1:

• the method of moments with an approximation by water-bags,

• the multi-water-bag model with the scheme of Godunov,

• the Particle-In-Cell (PIC) method, considered as the reference.

Before the shock, the two fluid methods describe precisely the solution. After the shock, they describe the main part of the solution, but cannot catch the filaments. More test cases will be presented at the conference, with higher N .

Figure 1 :

 1 Figure 1: Landau damping. Multi-water-bag model, method of moments and PIC method compared to the exact solution.

Figure 2 :

 2 Figure 2: Test case for N = 1 in the phase space at times T = 2 and T = 5. Multi-water-bag model (black circles), method of moments (empty squares) and PIC method (dots). 6


Ahmed Yureidini

Erwan Kerrien

Stéphane Cotin

Robust RANSAC-based blood vessel segmentation

come

DESCRIPTION OF PURPOSE

The segmentation of vascular structures is particularly valuable for diagnosis assistance, treatment and surgery planning. A wide range of applications may benefit from an improved blood vessel segmentation process: quantitative studies of pathologies (e.g. stenoses), automated vascular navigation, accurate blood flow computation and computer-based simulations for surgeons education. In such cases, a mere surface depiction is not enough and the extraction of the vessel medial axes is also required.

Vascular segmentation has resulted in a vast literature. [START_REF] Lesage | A review of 3D vessel lumen segmentation techniques: Models, features and extraction schemes[END_REF] Many previous works extract the vessel centerline tree, and claim to be robust to the kissing vessel issue: [START_REF] Wong | Probabilistic vessel axis tracing and its application to vessel segmentation with stream surfaces and minimum cost paths[END_REF][START_REF] Li | 3D multi-branch tubular surface and centerline extraction with 4D iterative key points[END_REF][START_REF] Friman | Multiple hypothesis template tracking of small 3D vessel structures[END_REF] two vessels may happen to be locally tangent, or a vessel may run along a dense structure, e.g. an aneurism or bone. Such works address this issue by using a shape prior. But the downside of this assumption is that the vessel cross-section is circular, which is not always true, especially for large vessels. Thus, a trade-off is usually made between accuracy and robustness. Our work aims at preserving this desired robustness against noise and topology issues while not compromising accuracy on extraction of the vessel surface and its centerline. This paper presents a tracking procedure which builds a vessel tree through successive robust fitting of cylinders to the image. Thereby, our dedicated tracking process: delineates the centerline of the vessel, supplies a local estimation of its radius and robustly extracts dense points on the vessel surface.

The reminder of this paper is organized as follows. Our RANSAC-Based Tracking (RBT) algorithm is introduced in Section 2.1. In Section 2.2, we describe the methodology of comparison between our proposal and the Multiple Hypothesis Tracking (MHT) algorithm. Section 3 discusses the results produced by both procedures on a data set of 10 patients, and conclusions are presented in Section 4.

METHOD

Vessel tracking algorithm

Efficient solutions for tracking blood vessels rely on the prior assumption that the shape is locally a cylinder. [START_REF] Wong | Probabilistic vessel axis tracing and its application to vessel segmentation with stream surfaces and minimum cost paths[END_REF][START_REF] Friman | Multiple hypothesis template tracking of small 3D vessel structures[END_REF] Nevertheless, previous works have hard times dealing with a large range of vessel diameters. [START_REF] Friman | Multiple hypothesis template tracking of small 3D vessel structures[END_REF] Moreover, the prior is often imposed too strongly, leading to vessels with circular cross-sections. The algorithm introduced in this section suggests rather than imposes a cylinder shape as a guide to select points, in a robust manner, on the local surface of the blood vessels. See also Figure 1 for a depiction of the main steps. . One fitting cylinder is found using RANSAC, per axis. (d) Best cylinder found, together with its consensus set (green dots) and the final points used to set the cylinder height.

RANSAC-Based Cylinder Estimation

The elementary tracking step consists in fitting a cylinder to the local vessel shape. We first suppose that this step was carried out and an initial cylinder of center C 0 , axis direction d 0 , radius r 0 and height h 0 was fitted. Hereafter, we describe how the next cylinder is found.

The center C and the direction d of the cylinder are first predicted: C = C 0 + (h 0 /4) d 0 , and d = d 0 . The radius r is bounded to be in [r 0 /2, 3r 0 /2]. An initial set of points are extracted using ray casting: N r rays, with a length of 4r 0 , equidistributed on a sphere, are casted from C; then, the point of minimum of the directional gradient is computed along each ray * . This extraction scheme is very similar to taking the points of minimal gradient along the columns of a Bounded Spherical Projection (BSP) image 2 (see Figure 1b).

Then, a cylinder is fitted to the points, using RANSAC. [START_REF] Chaperon | Extracting cylinders in full 3D data using a random sampling method and the Gaussian image[END_REF] Chaperon et al. [START_REF] Chaperon | Extracting cylinders in full 3D data using a random sampling method and the Gaussian image[END_REF] estimate the cylinder axis direction from the normals at the input points. This may fail at curved sections of vessels for example. Also, in our case, normals would be estimated as the image gradient at the extracted points. Such an input can be noisy at some places, making it hard to get a trustful cylinder direction. We therefore followed a strategy similar to Friman et al. [START_REF] Friman | Multiple hypothesis template tracking of small 3D vessel structures[END_REF] for robustness reasons (see Fig. 1c): N d directions are tested, equidistributed on the half-unit sphere oriented along d. Each direction is associated by RANSAC to a cylinder respecting the constraint on r, using a direct formula giving the cylinder center and radius knowing its direction. [START_REF] Beder | Direct solutions for computing cylinders from minimal sets of points[END_REF] Pruning the candidate cylinder directions drastically speeds up the process. The N d directions are reviewed by increasing angle with d 0 . As soon as a cylinder presents with an inlier rate above p inl , the algorithm stops and returns the cylinder. In case of failure, the algorithm returns the cylinder with the maximum inlier rate if it is at least p inl /2 or otherwise, the algorithm returns that no valid cylinder could be found. Finally, the center C of the returned cylinder is refined as the median point of its consensus set, and its parameters (center, radius and direction) are refined with Powell's optimization algorithm, to get a smoother sequence of cylinder directions. A new set of points is extracted by ray casting from C and the subset P of the inlier points is returned with the cylinder. Its height is set so as to encompass 75% of P (see Fig. 1d).

Tracking

The tracking is initialized manually by a point C 0 . Most vessels irrigating the brain are vertical: d 0 is the upward direction. A rough estimation of the radius r 0 is also supplied by the user as well as an inlier rate p inl . The tracking is stopped when no valid cylinder can be found or when it turns back. The other branches are added in a similar way. A point is placed manually on the new vessel, close to the bifurcation. The vector linking this point to the nearest point on the parent vessel centerline provides the direction for the new vessel. The vessel tree is built up by assembling the tracked branches. As a result, only one point per vessel is necessary, hence according to Schaap et al. [START_REF] Schaap | Standardized evaluation methodology and reference database for evaluating coronary artery centerline extraction algorithms[END_REF] classification, our method belongs to the category of algorithms with minimal user-interaction (category 2). For comparison, the MHT algorithm is assigned to category 3: interactive extraction, since a seed point, an initial tracking direction and an estimate of the vessel radius are required.

The algorithm output is thereafter threefold: the centerline is provided by the centers of the successive cylinders, each cylinder gives a robust estimation of the local vessel radius, and a dense sampling of the vessel surface is provided by the cylinder inliers.

Validation: Material and methods

Image Data

A set of 10 patient data was used for validation. Each patient data set consisted of a 3DRA acquired on a vascular C-arm (Innova 4100, GE Healthcare) during the intra-arterial injection of the internal carotid artery. Each 3DRA volume presents as a 512 3 isotropic voxel cube, between 0.18-0.22 mm voxel size. A total of 744 vessels and between 55-96 vessels per patient were examined.

The complex morphology of the vasculature makes any comparison of segmentation outcomes difficult. In the case of multibranch tracking algorithms, it is necessary to set branches into correspondence at each bifurcation, and to handle possible missing branches, which is difficult to automate. In our case of single branch tracking algorithms, the algrithms to compare may follow different branches at bifurcations, which is very hard to distinguish from tracking errors. As a result, only comparisons of vessel sections in-between bifurcations. For the sake of comparison, the algorithms were both initialized with the same four points, placed manually on each vessel section: one seed point C 0 was placed at the proximal end of the vessel; the vessel direction d 0 was estimated by placing a second point around the vessel axis, a little distal from C 0 ; and the vessel radius r 0 was estimated as the distance between two diametrically opposed points placed on the vessel surface, in a cut plane showing the vessel cross-section around C 0 .

The criterion to include a vessel in the validation set was that it could be tracked visually through successive cut planes in the 3DRA volume data. Only the internal carotid artery (the main vessel) was excluded because its radius is too large for the MHT algorithm used as the reference algorithm (see below). Note that our tracking algorithm gives visually excellent result on this large artery with the same parameters.

Parametrization of our algorithm

The RBT procedure was run for each vessel with the following parameters: the inlier rate threshold p inl was 70%, a relative threshold equal to 10% the local cylinder radius was taken on the distance to the fitted cylinder for the outlier rejection test; the number of rays thrown to extract points was N r = 162, and the number of axis directions tested was N d = 81. A minimum of 220 tests and a maximum of 500 tests were set for RANSAC. Note that the same parameters were used for all vessels, let they be rather large or small. Though our segmentation algorithm only needs one point per vessel, C 0 , d 0 and r 0 , as estimated above, were used to initialized our tracking procedure in order to be in the same conditions as the reference method.

Reference method

For comparison, the vessels were segmented with the Multiple Hypothesis Tracking (MHT) Algorithm, [START_REF] Friman | Multiple hypothesis template tracking of small 3D vessel structures[END_REF] available on MevisLab software † under the TubularTracking module. The manual seed C 0 and direction d 0 were used to track each branch. The single hypothesis tracking option was set (one single vessel to track). The tracking step was set to 0.1 the radius of current tracking radius in the MHT, the maximum search angle angle between 2 tracking steps to 85 • and the number of search angles was 5. A maximum of 2000 iterations were allowed and all other parameters initialized to default values.

The MHT algorithm proved to be quite sensitive to the allowable bounds set onto the radius. To have MHT working on the most possible vessels, the radius was allowed to vary between 0.25 to 1.25 times the estimated radius r 0 at the seed point.

Besides, since the reference MHT method has no efficient stopping criteria, the resulting centerline was cut manually where the vessel visual tracking was not reliable anymore.

Validation protocol

Since we noticed a certain amount of tracking failures for MHT, we first visually assessed the tracking success. A tracking was considered as successful if the resulting centerline remained inside the vessel of interest. Otherwise, the result was considered as a failure. This first evaluation resulted in a sucess rate in percents (see Table 1).

However, the first examination was tough for MHT. In some cases, rated as failures, the MHT successfully tracked the proximal section of the vessel, but leakage was commonly observed where neighboring dense structures misled the tracking. Since the successfully tracked portion was sometimes quite long, a second experiment was set up. Note that no leakage was noticed for RBT on the entire data set.

In this second experiment, these MHT tracking results, considered as failures in a first place, were visually cut before the problem occurred. According to the length of the detected centerline, results, both for MHT and RBT, were classified as: short, medium, long if the centerline was respectively correct along 1-33%, 33-66%, more than 66% of the visually estimated length of the targeted vessel. A second evaluation was therefore performed, recomputing the success rate, while only taking into account the medium and long trackings. Moreover, successful tracking lengths were computed and compared (see Table 2).

A third evaluation was carried out, focusing on a quantitative assessment of the centerline detection accuracy. The Average Symmetric Surface Distance (ASSD) 8 was used. Only vessels labelled as medium or long were considered, since we aimed at reporting figures on significant trackings. The ASSD measure translates the mean error between two centerlines. Its computation required to find the centerline sections that were potentially common to both. Each point on the centerline of the ground truth (MHT) was selected as a measure point if the test centerline (RBT) came enough in its vicinity. Each measure point {P i } i=1,...,N was associated with its closest point {Q i } i=1,...,N on the test centerline. Symmetrically, measure points {Q j } j=1,...,M were selected on the test centerline, and their corresponding point {P j } j=1,...,M were computed on the ground truth. The ASSD was thereafter computed as follows:

ASSD = 1 N + M    ∑ i d(P i , Q i ) + ∑ j d(P j , Q j )    (1)
where d(P, Q) is the Euclidean distance between points P and Q.

Two ASSD measures were carried out with different distance thresholds: first, a fixed value of 3 voxels and second, the mean radius of the targeted vessel was taken as an adaptive threshold. Indeed, the former is suitable for large vessels and the latter is a more appropriate choice to consider the scale of vessels and avoid taking into account erroneous centerlines at leaking portions on small vessels. Moreover comparing both measures enables to assess a potential influence of the distance threshold on the result. For the latter measure, the estimation of the mean radius was given by our algorithm, since it appeared more trustful than MHT's estimation (see Fig. 4).

RESULTS AND DISCUSSIONS

Table 1 reports on the first experiment. A total number of 744 vessel sections were tracked on the 10 patient data included in our experiments. For the visual assessment, the overall performance of RBT algorithm was 94%, confirming its efficiency and strength. In constrast, the MHT algorithm reached a success rate of 65%. This rather poor performance for MHT is the resultant of a lack of stopping criterion, which made MHT keep on tracking, even wrong vessels, when the target vessel became hardly visible. Leakages, due to neighboring structures, were the cause of many failures. On the contrary, the RBT procedure dealt well with this kissing vessel problem. Fig. 2 illustrates the difficulty for MHT to sometimes handle the kissing vessel problem and the robustness of our method to such topological difficulties. 1: Success rate of MHT and RBT on a data set of 10 patients. A total of 744 test vessels were used to evaluate both algorithms. The average success rate for MHT and for our algorithm is respectively 65% and 94%. The second evaluation took into account the actual tracked length, to assess the success rate. In particular, short trackings that were counted as successes in the first evaluation, were here discarded. Table 2 summarizes its results. As expected, MHT increased its performance up to 69%, since some trackings were previously rated as failures due to leakage issues, though they were correct on the proximal vessel section. On the contrary, RBT performance slightly dropped to 89%. This abatement could be mostly explained by RBT stopping prematurely due to drops in image intensity, as illustrated in Fig. 3.

When both algorithms succeeded, on average, RBT went further (725.5 mm) than MHT (688.6 mm). Nonetheless, MHT definitevely delivered a lesser number of successful trackings which impacted negatively the average tracked length per patient up to 778.6 mm against 1446.1 mm for RBT. Table 2 also reports quantitative ASSD measures. The average ASSD between both centerlines, with both 3 voxels and the mean radius of the targeted vessel, was below one voxel. This showed that the proposed method produced similar results as the MHT.

On one hand, our RBT algorithm showed the same behaviour over all cases. The most important observation is that it handled nicely the kissing vessel problem. Vessels of varied diameters could also be tracked with the same set of parameters. Nevertheless, our proposal reached its limits when struggling with drops in image intensity. On the other hand, the MHT algorithm coped especially well with gaps but it was very sensitive to both bad initialization and radius variation along blood vessels. A propos this limitation, the user-defined radius estimation bound could lead to a saturation of the radius estimation and ineluctably to a wrong estimation of the centerline as pointed out in Fig. 4. In contrast, RBT used the last estimated radius to restrain the estimate. Thereby, a good behaviour was noticed when contending with variations of density along vessels while the MHT centerline estimation was heavily disturbed (see Fig. 5).

One of the reasons we identified as potentially impeding MHT performances is the fixed tracking step that it uses. The increase in resolution in 3DRA implies that we try to track very small vessels (0.5 mm in diameter or less). Such vessels present with a very high tortuosity and potentially acute bends. RBT defines the tracking step as the quarter of the height of the last fitted cylinder. But since the height is adapted to the longitudinal extension of the inlier point set, fitted cylinders are shorter in twisted and bent portions of the blood vessels, thus naturally adapting the tracking step. 2: Success rate of MHT and RBT algorithms with tracking results classified as medium and long (M & L). The mean success rates for MHT and for our RBT algorithm are respectively 69% and 89%. Furthermore, the mean ASSD, with a distance threshold of 3 voxels (µ <3 ) and the mean radius of the targeted vessel (µ <R ), between the extracted centerlines by both methods is below one voxel (0.17-0.18 mm). The comparison length (CL) corresponds to the total length of centerline used to compute ASSD on vessels successfully tracked by both methods. On average the total tracked lengths (TL) for MHT and RBT are respectively 778.6 mm and 1446.1 mm.

CONCLUSION

The general context of our work is interventional neuroradiology. While most previous works present segmentation results on CTA or MRA data, we only addressed the segmentation of 3D Rotational Angiography (3DRA) volume. Indeed, 3DRA is the modality of choice for interventional radiologists that routinely acquire 3DRA data before, during and after the treatment. While 3DRA is arguably of a much higher quality than CTA and MRA, which would supposedly ease the segmentation task, noise and artefacts are still present. Moreover, even in 3DRA, and as in MRA and CTA, our aim remains the segmentation of vessels whose radius is about the voxel size. A state of the art MHT vessel segmentation algorithm 4 was used as a reference for validation. The performance of MHT were below our expectations, when applied to 3DRA, which supports our claim that 3DRA presents with specific issues regarding blood vessel segmentation: improved blood vessel visibility implies more kissing vessels, and vessels with a much stronger tortuosity.

The proposed algorithm for blood vessel segmentation, called RBT, builds a vessel tree structure, without imposing the vessels to be of circular cross-section. Furthermore, robustness is incorporated in the overall tracking procedure by using RANSAC to handle outliers in candidate points at the vessel surface. The results demonstated the strength of the proposed method as well as its ability to accurately detect centerlines. In particular, RBT proved to be able to capture very complex vascular topologies and to be very robust to the kissing vessel issue. In addition, it tracked very tortuous tiny blood vessels. As such, it improves upon previous blood vessel tracking algorithms.

Besides, RBT is also able to continue vessel tracking past bifurcations. After tracking a vessel, a child vessel is segmented by providing a seed point on the parent vessel. This point is connected to its closest point on the centerline of the parent vessel, which provide a initial direction for the tracking. The initial radius is the same as that of the parent vessel at the connecting point. This rough initialization is possible due to RBT being very robust. It did not work with MHT, that required tighter initial parameter values. Furthermore, RBT was even robust enough for us to use the same set of parameters (aside from center, direction and radius) for all patients and vessel. These characteristics of our method minimize user interaction and allow for segmentation of vessels of interest to the user. This counterbalances the fact that RBT does not automatically handles bifurcations. We would even argue that it is a clinical advantage since the radiologist is often interested in only a part of the vasculature (the one related to the pathology), that he/she would probably like to choose and control. Anyway, user interaction would be necessary, even with a so-called fully automated procedure, to correct for potential errors. Finally, the proposed segmentation algorithm seems to be easily extended to other modalities. Preliminary outcomes showed a promising behaviour of our algorithm on CTA data.

A path to improve RBT relies on a recognized strength of MHT. MHT allows the tracking to locally degrade the fitting score, as long as a raise occurs. As a result, MHT resists to local drops in the density along a vessel. However, this strategy prevents MHT for stopping per se as it too greedly follows hints of vesselness on the data. Our method stops as soon as no cylinder is found with the minimum percentage of inliers and may occasionnaly stops prematurely. In order to increase the robustness of our proposal, we need to integrate an exploration phase during the tracking procedure. A difficult compromise would have to be found between our stopping criteria and a greedy MHT exploration.

A final word should be said on the segmentation characteristic of RBT. It is a tracking algorithm since the vessel centerlines are extracted, building a vessel tree with little manual interaction. But, as MHT, the local vessel radius is also estimated. However, RBT goes further since it also provides a set of points that reliably lie on the local vessel surface. These points provide a rather dense sampling of the vessel surface. Our current and future work aim at providing a complete vessel implicit surface representation that interpolate these points.

Figure 1 :

 1 Figure 1: Outline of RBT vessel tracking algorithm. (a) Cut plane through the original 3DRA data; (b) Parameter update: Starting from C 0 (green dot) and d 0 (white arrow), the new center C (red dot) is found. d = d 0 . Candidate points (white dots) are extracted at the vessel surface by casting rays in the volume, from C. Only the points around the displayed cut plane are shown. (c) N d directions are tested for the cylinder axes (yellow arrows). One fitting cylinder is found using RANSAC, per axis. (d) Best cylinder found, together with its consensus set (green dots) and the final points used to set the cylinder height.

  (a) MHT result (b) RBT result (c) Full result (32 arteries)

Figure 2 :

 2 Figure 2: Comparison of MHT (left) and our algorithm (middle) with regard to a kissing vessel case on Patient 1: a blood vessel runs along the aneurism sack. Our RBT algorithm successfully handles this case, while MHT gets lost in the aneurism. (right) 32 arteries segmented on this patient data with our algorithm.

Figure 3 :

 3 Figure 3: MHT deals well with drops in image intensity. (left) MHT result; (middle) RBT result: the tracking stops too early. (right) a cut plane along the MHT centerline shows the drop in image intensity in the vessel. The tracking was performed on a branch of the Middle Cerebral Artery of Patient 10.

Figure 4 :

 4 Figure 4: Influence of user-defined radius estimation bound on MHT tracking where underestimation of the vessel radius leads to an innacurate outcome on Patient 5: (left) the upper bound is reached leading to leakage issues and (middle) computation with larger values avoid this problem. (right) Our RBT algorithm successfully tracks the vessel without resorting to this information.

  (a) Cut plane along the centerline of Fig.4a (b) Cut plane along the centerline of Fig.4b

  (c) Cut plane along the centerline of Fig.4c

Figure 5 :

 5 Figure 5: Plane views of Fig.4: the detection of the centerline for the MHT on a fuzzy portion of 0.8 mm of diameter, with saturation of the radius (left) and without saturation (middle), is disturbed whereas our algorithm detects correctly the centerline.

† http://www.mevislab.de

Acknowlegments

The authors would like to express their gratitude to Pr. René Anxionnat who provided the data for evaluation. Figures 2a,2b, 3 (left), 3 (middle), and 4, were produced with MeVisLab software (http://www.mevislab.de). Figures 2c,3 (right), and 5 were generated with Mayavi software (http://mayavi.sourceforge.net).