Mohammed Lemou
email: mohammed.lemou@univ-rennes1.fr

AND Ana Maria Luz

Florian Méhats
email: florian.mehats@univ-rennes1.fr

NONLINEAR STABILITY CRITERIA FOR THE HMF MODEL

published or not. The documents may come

Nonlinear stability criteria for the HMF Model

Mohammed Lemou, Ana Maria Luz, Florian

Introduction and main result

1.1. The HMF model. In this paper, we are interested in the nonlinear stability of a class of inhomogeneous steady state solutions to the Hamiltonian mean-field (HMF) model [START_REF] Messer | Statistical Mechanics of the Isothermal Lane-Emden Equation[END_REF][START_REF] Antoni | Clustering and relaxation in Hamiltonian long-range dynamics[END_REF]. The HMF system is a kinetic model describing particles moving on a unit circle interacting via an infinite range attractive cosine potential. This model has been used as a toy-model of the Vlasov-Poisson system in the physical community, for the study of non equilibrium phase transitions [START_REF] Chavanis | Lynden-Bell and Tsallis distributions for the HMF model[END_REF][START_REF] Staniscia | Out-of-equilibrium phase transitions in the HMF model : a closer look[END_REF][START_REF] Antoniazzi | Nonequilibrium tricritical point in a system with long-range interactions[END_REF][START_REF] Ogawa | Precise determination of the nonequilibrium tricritical point based on Lynden-Bell theory in the Hamiltonian mean-field model[END_REF], of travelling clusters [START_REF] Barré | Small traveling clusters in attractive and repulsive Hamiltonian mean-field models[END_REF][START_REF] Yamaguchi | Construction of traveling clusters in the Hamiltonian mean-field model by nonequilibrium statistical mechanics and Bernstein-Greene-Kruskal waves[END_REF] or of relaxation processes [START_REF] Yamaguchi | Stability criteria of the Vlasov equation and quasi-stationary states of the HMF model[END_REF][START_REF] Barré | The Vlasov equation and the Hamiltonian mean-field model[END_REF][START_REF] Chavanis | Dynamics and thermodynamics of a simple model similar to self-gravitating systems : the HMF model[END_REF]. The dynamics of perturbations of inhomogeneous steady states of the HMF model has been investigated in [START_REF] Barré | Dynamics of perturbations around inhomogeneous backgrounds in the HMF model[END_REF][START_REF] Barré | Algebraic damping in the one-dimensional Vlasov equation[END_REF] and the formal linear stability of steady states has been studied in [START_REF] Campa | Inhomogeneous Tsallis distributions in the HMF model[END_REF][START_REF] Ogawa | Spectral and formal stability criteria of spatially inhomogeneous solutions to the Vlasov equation for the Hamiltonian mean-field model[END_REF][START_REF] Barré | On the neighborhood of an inhomogeneous stable stationary solution of the Vlasov equation -Case of an attractive cosine potential[END_REF]. In particular, a simple criterion of linear stability has been derived in [START_REF] Ogawa | Spectral and formal stability criteria of spatially inhomogeneous solutions to the Vlasov equation for the Hamiltonian mean-field model[END_REF]. Our aim here is to prove the nonlinear stability of inhomogeneous steady states under the same criterion, by adapting the techniques developed in [START_REF] Lemou | Orbital stability of spherical galactic models[END_REF][START_REF] Lemou | A new variational approach to the stability of gravitational systems[END_REF] for the 3D Vlasov-Poisson system. However, we emphasize that the steady state solutions to the Vlasov-Poisson system studied in [START_REF] Lemou | Orbital stability of spherical galactic models[END_REF] are compactly supported. Here this assumption is not needed and a weaker assumption is made in the case of the HMF model, as we will see later on. Note finally that the long-time validity of the N-particle approximation for the HMF model has been investigated in [START_REF] Caglioti | Long time estimates in the mean field limit[END_REF][START_REF] Caglioti | Quasi-stationary states for particle systems in the mean-field limit[END_REF] and the Landau-damping phenomenon near spatially homogeneous state has been studied recently in [START_REF] Faou | Landau damping in Sobolev spaces for the Vlasov-HMF model[END_REF].

In the HMF model, the distribution function of particles f pt, θ, vq solves the initialvalued problem

B t f `vB θ f ´Bθ φ f B v f " 0, pt, θ, vq P R `ˆT ˆR, (1.1)
f p0, θ, vq " f init pθ, vq ě 0,

The authors acknowledge support by the ANR project Moonrise ANR-14-CE23-0007-01. The work of A. M. Luz was supported by the Brazilian National Council for Scientific and Technological Development (CNPq) under the program "Science without Borders" 249279/2013-4.

where T is the flat torus r0, 2πs and where the self-consistent potential φ f associated to a distribution function f is defined by

φ f pθq " ´ż 2π 0 ρ f pθ 1 q cospθ ´θ1 qdθ 1 , ρ f pθq " ż R f pθ, vqdv. (1.
2)

The so-called magnetization is the two-dimensional vector defined by

M f " ż 2π 0 ρ f pθqupθqdθ,
with upθq " pcos θ, sin θq T (1.3) and we have φ f pθq " ´Mf ¨upθq.

(1.4)

The Cauchy problem for (1.1) is much simpler than the one for the Vlasov-Poisson system, since the interaction kernel is smooth, and it can be shown that the HMF model is well-posed in the natural energy space. The following quantities are invariant during the evolution:

-the Casimir functions ij Gpf pθ, vqqdθdv

for any function G P C 1 pR `q such that Gp0q " 0; -the nonlinear energy Hpf q " 1 2

ij v 2 f pθ, vqdθdv `1 2 ż ρ f pθqφ f pθqdθ " 1 2 ij v 2 f pθ, vqdθdv ´1 2 M f ¨ż ρ f pθqupθqdθ " 1 2 ij v 2 f pθ, vqdθdv ´1 2 |M f | 2 ;
(1.6)

-the total momentum ij vf pθ, vqdθdv.

(1.7)

Moreover, the HMF system enjoys the Galilean invariance, that is, if f pt, θ, vq is a solution, then so is f pt, θ `v0 t, v `v0 q, for v 0 P R.

1.2. Statement of the main result. We consider a stationary state of the form f 0 pθ, vq " F pe 0 pθ, vqq, with e 0 pθ, vq " v 2 2 `φ0 pθq, (1.8) and where the potential associated to f 0 according to (1.2) takes the form φ 0 pθq " ´m0 cos θ, with m 0 ą 0.

Here F is a given function satisfying the following assumption.

Assumption 1.1. The function F is a C 0 function on R satisfying the following properties. It is a C 1 function on p´8, e ˚q, for some e ˚P R Y t`8u, with F 1 ă 0 on this interval. We also assume that F peq " 0 for e ě e ˚when e ˚is finite, and that lim eÑ`8 F peq " 0 if e ˚" `8. We denote by F ´1 its inverse function, which is a C 1 function defined from p0, sup F q onto p´8, e ˚q. The function f 0 given by (1.8) is supposed to belong to the energy space L 1 pp1 `|v| 2 qdθdvq. Moreover, in the case e ˚ă `8 and m 0 " e ˚, we assume further that ş m 0 ´m0 logpm 0 ´eqF 1 peqde ă `8.

Examples. All the following typical examples that can be found in the literature fulfill our Assumption 1.1:

(i) Maxwell-Boltzmann distributions [START_REF] Chavanis | Dynamics and thermodynamics of a simple model similar to self-gravitating systems : the HMF model[END_REF], F peq " A expp´βeq.

(ii) Polytropic distributions with compact support [START_REF] Campa | Inhomogeneous Tsallis distributions in the HMF model[END_REF], F peq " Ape ˚´eq 1 q´1

`with q ą 1. We used the usual notation x `" maxp0, xq.

(iii) Polytropic distributions with non compact support [START_REF] Campa | Inhomogeneous Tsallis distributions in the HMF model[END_REF], F peq " Ape 0 `eq

1 q´1 ẁith 1 3 ă q ă 1. (iv) Lynden-Bell distributions [11], F peq " A 1`B exppβeq . Remark 1.2. Note that Assumption 1.1 implies in particular that f 0 P L 8 since }f 0 } L 8 ď F p´m 0 q.
It is also clear that e ˚is finite if and only if f 0 is compactly supported. We finally note that we must have e ˚ą ´m0 , otherwise f 0 " 0 and this contradicts the assumption m 0 ą 0.

Our aim is to prove the orbital stability of such steady state under the following criterion.

Assumption 1.3 (Nonlinear stability criterion). We will assume that f 0 satisfies the following criterion

κ 0 ă 1, with κ 0 " ż 2π 0 ż `8 ´8 ˇˇF 1 pe 0 pθ, vqq ˇˇ¨ż D pcos θ ´cos θ 1 qpe 0 pθ, vq ´φ0 pθ 1 qq ´1{2 dθ 1 ż D pe 0 pθ, vq ´φ0 pθ 1 qq ´1{2 dθ 1 ‹ ‹ ' 2 dθdv,
(1.9) where D " θ 1 P T : φ 0 pθ 1 q ă e 0 pθ, vq ( .

Remark 1.4. Direct computations show that our criterion κ 0 ă 1 is the same as the one derived in [START_REF] Ogawa | Spectral and formal stability criteria of spatially inhomogeneous solutions to the Vlasov equation for the Hamiltonian mean-field model[END_REF], that is

0 ă 1 `ij F 1 pe 0 pθ, vqq cos 2 θdθdv ´4 ? m 0 ż m 0 ´m0 Kpkpeqq ˆ2Epkpeqq Kpkpeqq ´1˙2 F 1 peqde ´4 ? m 0 ż `8 m 0 Kp1{kpeqq kpeq ˆ2kpeq 2 Ep1{kpeqq Kp1{kpeqq `1 ´2kpeq 2 ˙2 F 1 peqde,
with kpeq " ´e`m 0 2m 0

¯1{2

and where Kpkq and Epkq are respectively the complete elliptic integrals of first and second kinds, see e.g. [START_REF] Barré | Dynamics of perturbations around inhomogeneous backgrounds in the HMF model[END_REF].

Before stating our main result, we first recall the usual notion of rearrangement which we adapt here to functions defined on the domain T ˆR. For any nonnegative function f P L 1 pT ˆRq, we define its distribution function as µ f psq " ~tpθ, vq P T ˆR : f pθ, vq ą su~, for all s ě 0, (1.10) where |A| denotes the Lebesgue measure of a set A. Note that µ f p0q may be infinite, but µ f psq is finite for s ą 0. Let f 7 be the pseudo-inverse of the function µ f , defined by f 7 psq " inf tt ě 0, µ f ptq ď su " sup tt ě 0, µ f ptq ą su , for all s ě 0 with, in particular, f 7 p0q " }f } L 8 P R Y t`8u and f 7 p`8q " 0. It is well known that µ f is right-continuous and that for all s ě 0, t ě 0,

f 7 psq ą t ðñ µ f ptq ą s. (1.11)
Next, we define the rearrangement f ˚of f by

f ˚pθ, vq " f 7 ´ffl ffl fflBp0, a θ 2 `v2 q X T ˆRffl ffl ffl ¯,
where Bp0, Rq denotes the open ball in R 2 centered at 0 with radius R.

Our main result is the following theorem.

Theorem 1.5. Let f 0 be a steady state of the form (1.8) satisfying Assumptions 1.1 and 1.3. There exists δ ą 0 such that, for all

f P L 1 `p1 `|v| 2 qdθdv ˘satisfying |M f ´Mf 0 p¨´θ f q | ă δ, we have }f ´f0 p¨´θ f q} 2 L 1 ď C ˆż Hpf q ´Hpf 0 q `Cp1 `}f } L 1 q}f ˚´f 0 } L 1 `C ż `8 0 s 2 ´f 7 0 psq ´f 7 psq ¯`ds `C ż `8 0 µ f 0 psq 2 β f ˚,f 0 psqds ˙, (1.12)
where β f ˚,f 0 psq " ~tpθ, vq P T ˆR : f ˚pθ, vq ď s ă f 0 pθ, vqu~, for all s ě 0, and where C is a positive constant depending only on f 0 . The parameter θ f is defined by M f " |M f |pcos θ f , sin θ f q T , where M f is given by (1.3). In particular, if f 0 is a compactly supported steady state, then (1.12) can be replaced by

}f ´f0 p¨´θ f q} 2 L 1 ď C ˆż Hpf q ´Hpf 0 q `Cp1 `}f } L 1 q}f ˚´f 0 } L 1 ˙. (1.13)
The proof of this theorem is given in Section 5 and uses several steps which are developed in the following sections. In Section 2, we introduce the generalized rearrangements with respect to the microscopic energy, which enable to define a reduced energy function depending on the magnetization vector only. In Section 3 we show that, under the stability criterion κ 0 ă 1, the magnetization of the steady state is a strict local minimizer of this reduced energy function and, in Section 4, we use a result in [START_REF] Lemou | Extended rearrangement inequalities and applications to some quantitative stability results[END_REF] to establish a functional inequality that enables the control of f ´f0 . We finally end the proof of Theorem 1.5 in section 5.

1.3. Proof of the orbital stability of f 0 . In this subsection, we show how to derive a stability result for the HMF model directly from our main Theorem 1.5.

Corollary 1.6. Let f 0 be a steady state of the form (1.8) satisfying Assumptions 1.1 and 1.3. Then f 0 is orbitally stable in the energy space, i.e., for all ε ą 0, there exists η ą 0 such that the following holds. For all solution f ptq to the HMF model with initial data f init , that preserves the mass and the energy, we have: if }p1 `v2 qpf init ´f0 q} L 1 ď η , then }p1 `v2 qpf p¨´θ f q ´f0 q} L 1 ď ε, where θ f is defined by M f " |M f |pcos θ f , sin θ f q T and M f is given by (1.3).

Proof. We distinguish the two cases: e ˚ă `8 and e ˚" `8.

Case 1: e ˚ă `8 . In this case f 0 is compactly supported and we can apply (1.13), that is we have

}f ´f0 p¨´θ f q} 2 L 1 ď C ˆż Hpf q ´Hpf 0 q `Cp ż 1 `}f } L 1 q}f ˚´f 0 } ˙(1.14)
for all f satisfying |M f ´Mf 0 p¨´θ f q | ă δ. Let f init P L 1 `p1 `|v| 2 qdvdθ ˘be any initial data for the HMF equation (1.1) such that

}p1 `|v| 2 qpf init ´f0 q} L 1 ă η,
where 0 ă η ă minp1, δ{2q will be made precise later on. This implies in particular that

|M f init ´Mf 0 | ď }f init ´f0 } L 1 ă η ă δ{2, (1.15)
and then

|Hpf init q ´Hpf 0 q| " ˇˇˇij ˆv2 2 `φ0 pθq ˙pf init ´f0 qdθdv ´1 2 |M f init ´Mf 0 | 2 ˇˇď pm 0 `1qη.
Now the contractivity property of the rearrangement implies that }f init ´f 0 } L 1 ď }f init ´f0 } L 1 ă η and then

|Hpf init q ´Hpf 0 q| `Cp ż 1 `}f init } L 1 q}f init ´f 0 } ď rm 0 `1 `C p2 `}f 0 } L 1 qs η.
We then choose η such that

η ă min ´1, δ{2, rm 0 `1 `C p2 `}f 0 } L 1 qs ´1 δ 2 {p4Cq ¯. (1.16)
Let now f ptq be a solution to the HMF model with initial data f init . From the conservation properties of this model, to wit Hpf ptqq " Hpf init q and f ptq ˚" f init , and from (1.14) we then get

}f ptq ´f0 p¨´θ f ptq q} 2 L 1 ď C ˆHpf ptqq ´Hpf 0 q `Cp ż 1 `}f ptq} L 1 q}f ptq ˚´f 0 } ă δ 2 {4,
as long as |M f ptq ´Mf 0 p¨´θ f ptq q | ă δ. In fact we shall prove that we have |M f ptq ´Mf 0 p¨´θ f ptq q | ă δ, @t ě 0.

(1.17) Indeed, at t " 0 we have |M f p0q ´Mf 0 p¨´θ f p0q q| ă δ{2 by assumption on f init (see (1.15)). If at some time t we have |M f ptq ´Mf 0 p¨´θ f ptq q| ě δ, then by continuity in time there is some time t 0 such that |M f pt 0 q ´Mf 0 p¨´θ f pt 0 q q| " 2δ{3 ă δ. We thus get }f pt 0 q ´f0 p¨´θ f pt 0 q q} L 1 ă δ{2.

But this implies

2δ{3 " |M f pt 0 q ´Mf 0 p¨´θ f pt 0 q q | ď }f pt 0 q ´f0 p¨´θ f pt 0 q q} L 1 ă δ{2, which is a contradiction, and claim (1.17) is proved. We conclude from Theorem 1.5 that

}f ptq ´f0 p¨´θ f ptq q} 2 L 1 ď C ˆ|Hpf init q ´Hpf 0 q| `Cp ż 1 `}f init } L 1 q}f init ´f 0 } L 1 ˙, (1.18
) for all t ě 0. The orbital stability of the solution f ptq is then proved in the L 1 norm in a quantitative way, since the right-hand side of (1.18) goes to zero as }p1 v2 qpf init ´f0 q} L 1 goes to zero, as a consequence of the usual contractivity property of the rearrangement }f init ´f 0 } L 1 ď }f init ´f0 } L 1 . It remains to prove this stability in the whole energy norm. We argue by contradiction. Assume that there exists ε ą 0 and a sequence f n init such that }p1 `v2 qpf n init ´f0 q} L 1 Ñ 0 as n Ñ `8 and, for some

t n ą 0, we have inf θPr0,2πs }v 2 pg n p¨´θq ´f0 q} L 1 ą ε,
where g n " f n pt n q and f n is a solution of the HMF model associated with the initial data f n init . We have already shown the L 1 stability, which means that we have }g n p¨´θ gn q ´f0 } L 1 Ñ 0 as n Ñ `8. In particular, up to a subsequence, we have v 2 g n p¨´θ gn q Ñ v 2 f 0 almost everywhere as n Ñ 8. Now from the Brézis-Lieb lemma, we have

}v 2 g n p¨´θ gn q ´v2 f 0 } L 1 ´}v 2 g n p¨´θ gn q} L 1 `}v 2 f 0 } L 1 Ñ 0 (1.19)
as n Ñ `8. From the conservation of the energy and the convergence }p1`v 2 qpf n init f0 q} L 1 Ñ 0, we have Hpg n p¨´θ g n qq " Hpg n q " Hpf n init q Ñ Hpf 0 q.

(1.20)

Note that we have used the convergence of the magnetization vector

ˇˇM f n init ´Mf 0 ˇˇď }f n init ´f0 } L 1 Ñ 0.
We apply this inequality to g n ˇˇM g n p¨´θ g n q ´Mf 0 ˇˇď }g n p¨´θ g n q ´f0 } L 1 Ñ 0, and obtain from (1.20)

}v 2 g n p¨´θ g n q} L 1 ´}v 2 f 0 } L 1 Ñ 0.
Using (1.19), this implies that

}v 2 g n p¨´θ g n q ´v2 f 0 } L 1 Ñ 0,
and yields a contradiction.

Case 2: e ˚" `8 . In this case f 0 is not compactly supported and we shall use inequality (1.12) of Theorem 1.5. The quantity µ f 0 psq involved in (1.12) is no longer bounded and presents a singularity at s " 0. Therefore we shall need to prove the following claim: if

}f n init ´f0 } L 1 Ñ 0 then ż `8 0 s 2 ´f 7 0 psq ´f n7 init psq ¯`ds Ñ 0 and ż `8 0 µ f 0 psq 2 β f ni nit ,f 0 psqds Ñ 0, (1.21)
as n Ñ `8 up to the extraction of a subsequence. Once this claim is proved, the rest of the stability proof is exactly the same as in the case of a compactly supported steady state f 0 . Let us then prove claim (1.21). We start by proving the first limit of this claim. Assume that

}f n init ´f0 } L 1 Ñ 0. Since }f n7 init ´f 7 0 } L 1 pR `q ď }f n init ´f0 } L 1 Ñ 0, we deduce that s 2 ´f 7
0 psq ´f n7 init psq ¯`Ñ 0 as n Ñ `8 for almost every s ě 0, up to an extraction of a subsequence. But we have

s 2 ´f 7 0 psq ´f n7 init psq ¯`ď s 2 f 7 0 psq, and
ż `8 0 s 2 f 7 0 psq ă `8,
(see (3.4) for the second inequality). Therefore, by dominated convergence we can pass to the limit inside the integral and get the first convergence in claim (1.21). Now we prove the second limit of claim (1.21). Assume again that

}f n init ´f0 } L 1 Ñ 0, then ż `8 0 β f ni nit ,f 0 psqds " ij pf 0 ´f ni nit q `dθdv ď }f ni nit ´f 0 } L 1 ď }f n init ´f0 } L 1 Ñ 0 as n Ñ `8.
This means that β f ni nit ,f 0 psq Ñ 0 for almost every s ě 0, up to an extraction of a subsquence. We then conclude that the quantity µ f 0 psq 2 β f ni nit ,f 0 psq arising in (1.21) converges to 0 for almost every s ě 0 (up to an extraction). Therefore, to end the proof of the second limit in claim (1.21), it is sufficient to dominate this quantity by an L 1 function in s P p0, }f 0 } L 8 q uniformly in n. To this purpose, we observe that β f ni nit ,f 0 psq ď µ f 0 psq and then

µ f 0 psq 2 β f ni nit ,f 0 psq ď µ f 0 psq 3 , @ s ą 0.
To prove that the rhs of this inequality is integrable on R `, we write

ż `8 0 s 2 f 7 0 psqds " ż `8 0 s 2 ˜ż f 7 0 psq 0 dt ¸ds " ż `8 0 ˜żf 7 0 psqąt
s 2 ds ¸dt and using (1.11) we get

ż `8 0 s 2 f 7 0 psqds " ż `8 0 ˜ż0ďsăµ f 0 ptq s 2 ds ¸dt " 1 3 ż `8 0 µ f 0 ptq 3 dt.
Since from (3.4) we have ş `8 0 s 2 f 7 0 psqds ă `8, the proof of claim (1.21) is complete. This proves the orbital L 1 stability. To get the stability in the energy space, we proceed as in the case e ˚ă `8. This ends the proof of the orbital stability in all cases.

The reduced energy functional

The aim of this section is to introduce a reduced energy functional J p|M f |q which depends only on the modulus of the magnetization and which is such that J p|M f |q ´J pm 0 q (recall that M f 0 " pm 0 , 0q T with m 0 ě 0) is controlled by the relative nonlinear energy Hpf q ´Hpf 0 q, up to conserved quantities.

2.1. Generalized rearrangements with respect to the microscopic energy. Our purpose now is to define a generalized symmetric nonincreasing rearrangement with respect to the microscopic energy e " v 2 2 `φpθq, where the potential φ is a given C 8 function on T. We introduce the quantity a φ peq " It is readily seen that a φ is continuous on R, vanishes on p´8, min φs and is strictly increasing from rmin φ, `8q to r0, `8q. This enables to define its inverse a ´1 φ on r0, `8q. Note that, for all e P R,

ffl
4π ? 2pe ´max φq 1{2 `ď a φ peq ď 4π ? 2pe ´min φq 1{2 `, (2.2)
which implies, for all s P R `,

s 2 32π 2 `min φ ď a ´1 φ psq ď s 2 32π 2 `max φ. (2.3)
We now introduce the generalized rearrangement with respect to the microscopic energy.

Lemma (2.4)

The function α 1 peq and its derivative α 1 1 peq are represented on Figure 1. The proof of the following lemma is deffered to the Appendix. This function satisfies the following properties: (i) α 1 is a continuous nondecreasing function from R to R `and α 1 peq " 0 for e ď ´1.

(ii) α 1 is a strictly increasing and strictly convex C 1 function on r´1, 1q. Its derivative for e P p´1, 1q is given by

α 1 1 peq " 2 ? 2 ż θ 1 peq 0 pe `cos θq ´1{2 dθ. (2.6)
and its right-derivative at e " ´1 is equal to 2π.

(iii) α 1 is a strictly increasing and strictly concave C 1 function on p1, `8q. Its derivative for e P p1, 8q is still given by (2.6) and we have

α 1 peq " 4π ? 2e, α 1
1 peq " 2π a 2{e as e Ñ `8.

(iv) We have α 1 p1q " 16 and

α 1 1 peq " ´2 log |e ´1| as e Ñ 1. (2.7)
(v) The inverse α ´1 1 of the function α 1 : r´1, `8q Þ Ñ r0, `8q is a strictly increasing C 1 function, defined on r0, `8q, satisfying pα ´1 1 q 1 psq "

1 α 1 1 ˝α´1 1 psq
for s P R `zt0, 16u, with α 1 1 given by (2.6), and

pα ´1 1 q 1 p0q " 1 2π , pα ´1
1 q 1 p16q " 0.

2.2.

Reduction to a functional of the magnetization vector. In this subsection, we prove the following result.

Proposition 2.3. For all f P L 1 pp1 `|v| 2 qdθdvq, we have

J p|M f |q ´J p|M f 0 |q ď Hpf q ´Hpf 0 q `3}f } L 1 }f ˚´f 0 } L 1 (2.8) `1 32π 2 ż `8 0 s 2 ´f 7 0 psq ´f 7 psq ¯`ds
where, for all m P R `, J pmq "

m 2 2 `ż `8 ´8 ż 2π 0 ˆv2 2 `φ˙f ˚φ 0 dθdv (2.9)
with φpθq " ´m cos θ.

Proof. Writing the difference Hpf q ´Hpf 0 q between the nonlinear energies as

Hpf q ´Hpf 0 q " ij ˆv2 2 `φf ˙pf ´f0 q dθdv ´1 2 `|M f | 2 ´|M f 0 | 2 ˘´ij φ f pf ´f0 qdθdv " ij ˆv2 2 `φf ˙´f ´f ˚φf ¯dθdv `ij ˆv2 2 `φf ˙´f ˚φf ´f ˚φf 0 ¯dθdv `ij ˆv2 2 `φf ˙´f ˚φf 0 ´f0 ¯dθdv `1 2 |M f ´Mf 0 | 2 "I 1 `I2 `I3 `1 2 |M f ´Mf 0 | 2 ,
we organize the proof in three steps.

Step 1: Identification of J p|M f |q ´J p|M f 0 |q.

Let us first prove that

f 0 " f ˚φ0 0 , (2.10)
which amounts to proving that F peq " f 7 0 ˝aφ 0 peq, @e ě min φ 0 .

(2.11)

Recall that a φ 0 is invertible from rmin φ 0 , `8q to r0, `8q and denote G " F ˝a´1 φ 0 on r0, `8q. Recall also that F is assumed to be continuously decreasing. Hence, so is the function G and then it is standard that G 7 " G, see for instance [START_REF] Rakotoson | Réarrangement relatif, Un instrument d'estimations dans les problèmes aux limites[END_REF]. Now, for all t ě 0,

µ
v 2 2 `φ0 pθq ă a ´1 φ 0 pµ G ptqq *ffl ffl ffl ffl " a φ 0 ˝a´1 φ 0 pµ G ptqq " µ G ptq.
From this, we deduce that f 7 0 " G 7 " G, which gives (2.11) and ends the proof of (2.10). Consequently,

I 3 `1 2 |M f ´Mf 0 | 2 " ij ˆv2 2 `φf ˙´f ˚φf 0 ´f0 ¯dθdv `1 2 |M f ´Mf 0 | 2 " ij ˆv2 2 `φf ˙´f ˚φf 0 ´f ˚φ0 0 ¯dθdv `1 2 |M f ´Mf 0 | 2 " ij ˆv2 2 `φf ˙f ˚φf 0 dθdv ´ij ˆv2 2 `φ0 ˙f ˚φ0 0 dθdv `ż pφ 0 ´φf q ρ f 0 dθ `1 2 |M f ´Mf 0 | 2 " ij ˆv2 2 `φf ˙f ˚φf 0 dθdv `1 2 |M f | 2 ´ij ˆv2 2 `φ0 ˙f ˚φ0 0 dθdv ´1 2 |M f 0 | 2 .
We observe now that φ f can be written as φ f pθq " ´|M f | cospθ ´θM q for some θ M P T. Hence, by periodicity, we have

ij ˆv2 2 `φf ˙f ˚φf 0 dθdv `1 2 |M f | 2 " J p|M f |q,
where J is defined by (2.9), and the same holds for φ 0 . We thus have

I 3 `1 2 |M f ´Mf 0 | 2 " J p|M f |q ´J p|M f 0 |q.
Step 2: Positivity of I 1 .

We have, using Fubini,

I 1 " ij ˆv2 2 `φf ˙´f ´f ˚φf ¯dθdv " ij ˆv2 2 `φf ˙˜ż f 0 dt ´ż f ˚φf 0 dt ¸dθdv " ż `8 0 ¨ij f ąt ˆv2 2 `φf ˙dθdv ´ij f ˚φf ąt ˆv2 2 `φf ˙dθdv ‹ 'dt " ż `8 0 ¨ij Aptq ˆv2 2 `φf ˙dθdv ´ij Bptq ˆv2 2 `φf ˙dθdv ‹ 'dt
where, for all t ě 0, we have denoted

Aptq " ! pθ, vq P T ˆR : f ˚φf pθ, vq ď t ă f pθ, vq ) , Bptq " ! pθ, vq P T ˆR : f pθ, vq ď t ă f ˚φf pθ, vq
) .

Since f ˚φf is a decreasing function of v 2 2 `φf , we clearly have @pθ, vq P Aptq, @pθ 1 , v

1 q P Bptq, v 2 2 `φf pθq ą v 12 2 `φf pθ 1 q.
Moreover, from the equimeasurability of f and f ˚φf , we have |Aptq| " |Bptq|. Consequently, we obtain I 1 ě 0.

Step 3: Control of |I 2 | by }f ˚´f 0 } L 1 .

Let us first state an elementary result.

Lemma 2.4. Let φpθq " ´m cospθ ´θ0 q for pm, θ 0 q P R `ˆT.

´a´1 φ f psq ´min φ f ¯´f 7 psq ´f 7 0 psq ¯ds ´}φ f } L 8 }f 7 ´f 7 0 } L 1 ě ´ż `8 0 a ´1 φ f psq ´f 7 0 psq ´f 7 psq ¯`ds ´2}φ f } L 8 }f 7 ´f 7 0 } L 1 .
Using (2.3), we deduce that

I 2 ě ´1 32π 2 ż `8 0 s 2 ´f 7 0 psq ´f 7 psq ¯`ds ´3}φ f } L 8 }f 7 ´f 7 0 } L 1 .
We now conclude by observing that, for all θ P T, we have

|φ f pθq| ď |M f | |upθq| " |M f | ď }f } L 1 .
(2.13)

3. Study of the functional J .

In this section, we study the function J pmq defined for m P R `by (2.9), with φpθq " ´m cos θ. For e P R and m P R ˚, we recall that

a φ peq " α m peq " ? m α 1 ´e m ¯,
where α 1 was defined by (2.5). Clearly, (2.9) and (2.12) yield, for m ą 0,

J pmq " m 2 2 `ż `8 ´8 ż 2π 0 ˆv2 2 ´m cos θ ˙f 7 0 ˝αm ˆv2 2 ´m cos θ ˙dθdv " m 2 2 `m ż `8 0 f 7 0 psq α ´1 1 ˆs ? m ˙ds.
Proposition 3.1. The function J defined by (2.9) is a C 2 function on R `. Denoting φpθq " ´m cos θ, we have

J 1 pmq " m ´ij f ˚φ 0 pθ, vq cos θ dθdv (3.1)
and J 2 pmq " 1`ij pf 7 0 ˝aφ q 1 pepθ, vqq ¨cos θ ´ż θmpepθ,vqq

0 cos θ 1 pepθ, vq `m cos θ 1 q ´1{2 dθ 1 ż θmpepθ,vqq 0 pepθ, vq `m cos θ 1 q ´1{2 dθ 1 ‹ ‹ ‹ ' 2 dθdv, (3.2
) where epθ, vq " v 2 2 `φpθq and θ m is defined by (2.4).

From this Proposition and from (2.10), it is immediate to deduce the Corollary 3.2. Under Assumption 1.3, the magnetization m 0 of the stationary state f 0 is a strict local minimizer of J : one has J 1 pm 0 q " 0 and J 2 pm 0 q " 1 ´κ0 ą 0.

Proof of Proposition 3.1. e `m cos θ dθde.

By passing to the limit in this formula, we also get that J is differentiable at m " 0, with J 1 p0q " 0. Finally, coming back to the variables pθ, vq, we obtain (3.1).

In order to compute the second derivative of J , let us transform this expression into a more suitable one, using an integration by parts in e. We denote r e ˚" a ´1 φ ˝aφ 0 pe ˚q, where e ˚is defined in Assumption 1.1. By (2.10), we have f 7 0 ˝aφ " F ˝a´1 φ 0 aφ , this function being continuous on r´m, `8q, of class C 1 on r´m, `8qztm, r e ˚u, nonincreasing, and vanishes on rr e ˚, `8q. Therefore, in the case e ˚ă `8, one can directly integrate by parts to obtain ż `8 ´m f 7 0 ˝aφ peqb 1 φ peq de " ´ż `8 ´m pf 7 0 q 1 ˝aφ peqa 1 φ peqb φ peq de.

Now we deal with the case r e ˚" e ˚" `8. Since f 7 is a nonincreasing function on R `and belongs to L 1 pR `q, we deduce that f 7 psq Ñ 0 when e Ñ `8. Therefore, according to (3.7), we have f 7 0 ˝aφ peqb φ peq Ñ 0 when e Ñ `8, and the integration by parts giving (3.11) is also valid in the case e ˚" `8.

Consequently, we have

J 1 pmq " m `ż `8 ´m pf 7 0 q 1 ˝aφ peqa 1 φ peqb φ peqde " m `ż `8 0 pf 7 0 q 1 psqb φ ˝a´1 φ psqds " m `?m ż `8 0 pf 7 0 q 1 psq β 1 ˝α´1 1 ˆs ? m ˙ds.
Consider the function hpm, sq " pf 7 0 q 1 psq β 1 ˝α´1 1 ˆs ? m ˙.

Using again Lemma 2.2, we get that h is continuously differentiable with respect to m P R ˚for all m P R ˚zts 2 {32u, with

Bh Bm pm, sq " ´spf 7 0 q 1 psq 2m 3{2 α 1 1 ˝α´1 1 p s ? m q β 1 1 ˝α´1 1 ˆs ? m ˙.
Since |b 1 φ peq| ď a 1 φ peq, we deduce that

ˇˇˇB h

Bm pm, sq ˇˇˇÀ ´spf 7 0 q 1 psq, for all m P rm 1 , m 2 s, 0 ă m 1 ă m 2 .

We now claim that the fonction s Þ Ñ spf 7 0 q 1 psq belongs to L 1 pR `q.

(3.12) Indeed, since f 7 0 is decreasing, we have

ż r 0 s 2 f 7 0 psqds ě f 7 0 prq ż r 0 s 2 ds " r 3 3 f 7 0 prq.
Hence, using (3.4), we get

f 7 0 psq À 1 s 3 , @s ą 0.
In particular sf 7 0 psq Ñ 0 when s Ñ `8. On the other hand, the function f 7 0 " F ˝a´1 φ 0 is continuous on R `, of class C 1 and decreasing on r0, a φ 0 pe ˚qq, vanishing on ra φ 0 pe ˚q, `8q (with possibly a φ 0 pe ˚q " `8). Therefore we can perform the following integration by parts ´ż `8 0 spf 7 0 q 1 psqds "

ż `8 0 f 7 0 psqds ă `8.
This ends the proof of claim (3.12) and enables to conclude by dominated convergence that J 1 is continuously differentiable on R `and that

J 2 pmq " 1 `1 2 ? m ż `8 0 hpm, sqds `?m ż `8 0 Bh Bm pm, sqds " 1 `ż `8 0 pf 7 0 q 1 psq ¨1 2 ? m β 1 ˝α´1 1 ˆs ? m ˙´sβ 1 1 ˝α´1 1 ´s ? m 2m α 1 1 ˝α´1 1 p s ? m q 'ds " 1 `1 2m ż `8 ´m pf 7 0 ˝aφ q 1 peq a 1 φ peq `a1 φ peqb φ peq ´aφ peqb 1 φ peq ˘de.
Finally, observing from (3.9) and from

eb 1 φ peq `2m ? 2 ż θmpeq 0 pcos θq 2 pe `m cos θq ´1{2 dθ " 1 2 b φ peq that 1 2ma 1 φ peq `a1 φ peqb φ peq ´aφ peqb 1 φ peq " ´pb 1 φ peqq 2 a 1 φ peq `2? 2 ż θmpeq 0 pcos θq 2 pe `m cos θq ´1{2 dθ " 2 ? 2 ż θmpeq 0 ¨cos θ ´ż θmpeq 0 cos θ 1 pe `m cos θ 1 q ´1{2 dθ 1 ż θmpeq 0 pe `m cos θ 1 q ´1{2 dθ 1 ‹ ‹ ‹ ' 2 pe `m cos θq ´1{2 dθ,
we obtain (3.2) by coming back to the pθ, vq variables.

Control of f

Our previous analysis has allowed the control of the magnetization vector by the relative Hamiltonian and the relative rearrangements. It remains to control the whole distribution function f . To this aim we now write the relative energy in the following form:

Hpf q ´Hpf 0 q " ij ˆv2 2 `φf 0 ˙pf ´f0 q dθdv ´1 2 |M f ´Mf 0 | 2 . (4.1)
In particular, this means that the following quantity ij ˆv2 2 `φf 0 ˙pf ´f0 q dθdv is controlled and the problem is to show how this quantity controls f ´f0 . This task was achieved in the context of the gravitational Vlasov-Poisson system [START_REF] Lemou | Orbital stability of spherical galactic models[END_REF] using compactness arguments. Here we will rather use a functional inequality established in [START_REF] Lemou | Extended rearrangement inequalities and applications to some quantitative stability results[END_REF] to get a quantitative control of }f ´f0 } L 1 by this quantity, up to rearrangement terms depending only on f ˚and f 0 which are preserved by the flow. We emphasize that the steady states to Vlasov-Poisson system studied in [START_REF] Lemou | Orbital stability of spherical galactic models[END_REF] are compactly supported and this property was essential to successfully drive the stability analysis in this context. Here this assumption is not needed and a much weaker assumption is made in the case of the HMF model. More precisely, we have the following inequality:

Proposition 4.1. Let f 0 be given by (1.8) where F satisfies Assumption 1.1. Then, there exist a constant K 0 depending only on f 0 such that, for all f P L 1 pp1`|v| 2 qdvdθq we have

p}f ´f0 } L 1 `}f 0 } L 1 ´}f } L 1 q 2 ďK 0 ij ˆv2 2 `φf 0 ˙pf ´f0 q dθdv `m0 }f ˚´f 0 } L 1 `1 8π 2 ż `8 0 µ 0 psq 2 β f ˚,f 0 psqds, (4.2)
where β f ˚,f 0 psq " ~tpθ, vq P T ˆR : f ˚pθ, vq ď s ă f 0 pθ, vqu~, for all s ě 0.

Proof. We shall apply Theorem 1 in [START_REF] Lemou | Extended rearrangement inequalities and applications to some quantitative stability results[END_REF]. We use the rearrangement with respect to e 0 pθ, vq " v 2 2 `φf 0 and recall that the function a φ 0 is strictly increasing and a one-to-one function from rmin φ 0 , `8q to r0, `8q. Following [START_REF] Lemou | Extended rearrangement inequalities and applications to some quantitative stability results[END_REF], we introduce the functions

B 0 pµq " ż µ 0 a ´1 φ 0 psqds, @µ ě 0, (4.3)
and H 0 pµq " inf 0ăsďµ B 0 pµ `sq `B0 pµ ´sq ´2B 0 pµq s 2 .

Then from Theorem 1 in [START_REF] Lemou | Extended rearrangement inequalities and applications to some quantitative stability results[END_REF], we have

p}f ´f0 } L 1 `}f 0 } L 1 ´}f } L 1 q 2 ďKpf 0 q ij ˆv2 2 `φf 0 ˙pf ´f0 q dθdv `ż `8 0 a ´1 φ 0 p2µ f 0 psqqβ f ˚,f 0 psqds ´ż `8 0 a ´1 φ 0 pµ f 0 psqqβ f 0 ,f ˚psqds (4.4)
where

Kpf 0 q " 4 ż }f 0 } L 8 0 ds H 0 pµ f 0 psqq
, and (4.5) β f,g psq " meastpθ, vq P T ˆR; f pθ, vq ď s ă gpθ, vqu, @s ě 0. (4.6) Using the estimates (2.3) we then get from (4.4)

p}f ´f0 } L 1 `}f 0 } L 1 ´}f } L 1 q 2 ďKpf 0 q ij ˆv2 2 `φf 0 ˙pf ´f0 q dθdv `1 8π 2 ż `8 0 µ 0 psq 2 β f ˚,f 0 psqds `m0 ż `8 0 ´βf 0 ,f ˚psq `βf ˚,f 0 psq ¯ds (4.7)
Observing that ż `8 0 β f,g psqds " ij pg ´f q `dθdv, we get ż `8 0 ´βf 0 ,f ˚psq `βf ˚,f 0 psq ¯ds " }f ˚´f 0 } L 1 , and therefore

p}f ´f0 } L 1 `}f 0 } L 1 ´}f } L 1 q 2 ďKpf 0 q ij ˆv2 2 `φf 0 ˙pf ´f0 q dθdv `1 8π 2 ż `8 0 µ 0 psq 2 β f ˚,f 0 psqds `m0 }f ˚´f 0 } L 1 . (4.8)
To end the proof of inequality (4.2), it only remains to show that the quantity Kpf 0 q is finite. First we rewrite H 0 pµq as

H 0 pµq " inf 0ăsďµ B 0 pµ `sq `B0 pµ ´sq ´2B 0 pµq s 2 " inf 0ăsďµ ż 1 0
p1 ´λqppa ´1 φ 0 q 1 pµ `λsq `pa ´1 φ 0 q 1 pµ ´λsqqdλ.

Then, from the properties of a φ stated in Lemma 2.2, we claim that pa ´1 φ 0 q 1 pµ `λsq `pa ´1 φ 0 q 1 pµ ´λsq ě pa ´1 φ 0 q 1 pµq, for all 0 ď λ ď 1, 0 ă s ď µ. Indeed, if µ ď 16 ? m 0 then µ ´λs ď 16 ? m 0 , and since the function a ´1 φ 0 is a concave function on r0, 16 ? m 0 s, we have pa ´1 φ 0 q 1 pµ λsq ě pa ´1 φ 0 q 1 pµq. Therefore we have the desired claim in this case. Similarly, if µ ě 16 ? m 0 then µ `λs ě 16 ? m 0 , and since the function a ´1 φ 0 is a convex function on r16 ? m 0 , `8s, we have pa ´1 φ 0 q 1 pµ `λsq ě pa ´1 φ 0 q 1 pµq. Therefore the above claim holds in both cases. Using this claim, we then get

H 0 pµq ě 1 2a
Kpf 0 q ď 8 ż }f 0 } L 8 0 a 1 φ 0 `F ´1ptq ˘dt.
We then perform the change of variable e " F ´1ptq to get

Kpf 0 q ď 8 ż em 0 a 1 φ 0 peq|F 1 peq|de. (4.9)
Now we claim that the rhs integral in this inequality is finite. Indeed, assume first that e ˚ă `8. The only possible singularities in this integral are at e " m 0 and e " e ˚, since the function e Þ Ñ a 1 φ 0 peq|F 1 peq| is continuous on r´m 0 , `8qztm 0 , e ˚u. If we suppose that e ˚‰ m 0 , then we have a 1 φ 0 peq|F 1 peq| " a 1 φ 0 pe ˚q|F 1 peq| when e Ñ e ˚and, from Lemma 2.2 we have (for m 0 ă e ˚otherwise F vanishes in the neighborhood of m 0 ) a 1 φ 0 peq|F 1 peq| " C log |e´m 0 | when e Ñ m 0 . These two possible singularities are thus integrable (the first is integrable by assumption on F ).

If e ˚" m 0 then our Assumption 1.1 (i) ensures that ş em 0 a 1 φ 0 peq|F 1 peq|de is finite, since a 1 φ 0 peq " C log |e ´m0 | as e Ñ m 0 . Assume now that e ˚" `8. Using assertion (iii) of Lemma 2.2, we have a 1 φ 0 peq|F 1 peq| " C|F 1 peq|{ ? e as e Ñ `8, where C is a constant. But, as for (2.12), we have

ż ż F ˆv2 2 `φ0 pθq ˙dθdv " ż ż ¨ż F ´v2 2 `φ0 pθq 0 dt 'dθdv " ż `8 0 meas " F ˆv2 2 `φ0 pθq ˙ą t * dt " ż `8 0 a φ 0 pF ´1ptqqdt " ż `8 ´m0 a φ 0 peq|F 1 peq|de.
This implies that the integral ş `8 ´m0 a φ 0 peq|F 1 peq|de is convergent. From assertion (iii), we know that a 1 φ 0 peq|F 1 peq| " C 1 |F 1 peq|{ ? e ď C 1 |F 1 peq| ? e " C 2 a φ 0 peq|F 1 peq| for e large enough, where C 1 and C 2 are some positive constants. This proves the fact that the rhs integral of (4.9) is finite, and ends the proof of Proposition 4.1.

Proof of Theorem 1.5

We first insert identity (4.1) into inequality (4.2) and get

p}f ´f0 } L 1 `}f 0 } L 1 ´}f } L 1 q 2 ďK 0 " Hpf q ´Hpf 0 q `1 2 |M f ´Mf 0 | 2  `m0 }f ˚´f 0 } L 1 `1 8π 2 ż `8 0 µ 0 psq 2 β f ˚,f 0 psqds. (5.1)
We write M f " |M f |upθ f q where upθq " pcos θ, sin θq T , and denote by f 0 p¨´θ f q the fonction f 0 p¨´θ f qpθ, vq " f 0 pθ ´θf , vq. We then apply this inequality (5.1) to f 0 p¨´θ f q and get

p}f ´f0 p¨´θ f q} L 1 `}f 0 } L 1 ´}f } L 1 q 2 ď K 0 " Hpf q ´Hpf 0 q `1 2 |M f ´Mf 0 p¨´θ f q | 2  `m0 }f ˚´f 0 } L 1 `1 8π 2 ż `8 0 µ 0 psq 2 β f ˚,f 0 psqds.

Now we observe that

M f 0 p¨´θ f q " ż 2π 0 ρ f 0 pθ ´θf qupθqdθ " ż 2π 0 ρ f 0 pθqupθ `θf qdθ " pm 0 cospθ f q, m 0 sinpθ f qq T " m 0 upθ f q. Therefore p}f ´f0 p¨´θ f q} L 1 `}f 0 } L 1 ´}f } L 1 q 2 ď K 0 " Hpf q ´Hpf 0 q `1 2 p|M f | ´m0 q 2  `m0 }f ˚´f 0 } L 1 `1 8π 2 ż `8 0 µ 0 psq 2 β f ˚,f 0 psqds. (5.2)
Now we use Corollary 3.2 together with the fact that J is a C 2 function to conclude that there exist δ ą 0 and C ą 0 such that J pmq ´J pm 0 q ě Cpm ´m0 q 2 for all m P pm 0 ´δ, m 0 `δq.

Reporting this into estimate (5.2) yields

p}f ´f0 p¨´θ f q} L 1 `}f 0 } L 1 ´}f } L 1 q 2 ď K 0 " Hpf q ´Hpf 0 q `1 2C pJ p|M f |q ´J pm 0 qq  `m0 }f ˚´f 0 } L 1 `1 8π 2 ż `8 0 µ 0 psq 2 β f ˚,f 0 psqds.
for all f such that |M f | P pm 0 ´δ, m 0 `δq. Now using inequality (2.8), we get

p}f ´f0 p¨´θ f q} L 1 `}f 0 } L 1 ´}f } L 1 q 2 ď C " Hpf q ´Hpf 0 q `Cp1 `}f } L 1 q}f ˚´f 0 } L 1 ż `C ż `8 0 s 2 pf 7 0 psq ´f 7 psqq `ds `C ż `8 0 µ 0 psq 2 β f ˚,f 0 psqds  . (5.3)
for some positive constant C only depending on f 0 . To end the proof of Theorem 1.5, we observe that, from the inequality pa `bq 2 ě 1 2 a 2 ´b2 ,

p}f ´f0 p¨´θ f q} L 1 `}f 0 } L 1 ´}f } L 1 q 2 ě 1 2 }f ´f0 p¨´θ f q} 2 L 1 ´p}f 0 } L 1 ´}f } L 1 q 2 ě 1 2 }f ´f0 p¨´θ f q} 2 L 1 ´}f 0 ´f ˚}2 L 1 ě 1 2 }f ´f0 p¨´θ f q} 2 L 1 ´p}f 0 } L 1 `}f } L 1 q }f 0 ´f ˚}L 1 ě 1 2 }f ´f0 p¨´θ f q} 2 L 1 ´C p1 `}f } L 1 q }f 0 ´f ˚}L 1 with C " maxp1, }f 0 } L 1 q.
We then report this into (5.3) and get inequality (1.12) for all f such that |M f | P pm 0 ´δ, m 0 `δq.

Let us deduce (1.13) in the case where if f 0 is a compactly supported steady state. In this case, the support of f This enables to deduce (1.13) from (1.12) and this ends the proof of Theorem 1.5.

Appendix

Proof of Lemma 2.2. The proof of Item (i) is straightforward. Let us prove Item (ii).

It is already clear from (2.5) that α 1 is strictly increasing. In order to prove that α 1 1 peq is given by (2.6) for all e P p´1, 1q, we perform the change of variable u " cos θ in (2.5): and, for all e P p´1, 1q, ż 1

α 1 peq " ż 1 ´1 gpe,
´1 q e puq du "

2 ? 2 ? 1 ´e π.
Hence, by Brézis-Lieb's Lemma [START_REF] Lieb | Loss. Analysis[END_REF], we have q e Ñ q e 0 in L 1 p´1, 1q, for all e 0 P p´1, 1q, and using a generalized dominated convergence theorem as stated in [START_REF] Lemou | A new variational approach to the stability of gravitational systems[END_REF] (Appendix A), we deduce that α 1 is C 1 on p´1, 1q, with (5.4)

Performing again the change of variable u " cos θ in (5.4) yields (2.6). Now, we perform the change of variable t " u`e 1´u in (5.4) and get, for e P p´1, 1q, From this expression, we clearly see that α 1 1 is strictly increasing, which yields the convexity of α 1 on p´1, 1q. We also deduce that the right-derivative of α 1 at e " ´1 is finite and its value is given by 2 ż `8 0 dt p1 `tq ? t " 2π.

Item (iii) is an easy consequence of the following expression, valid for e ą 1:

α 1 peq " 2 ? 2 ż 2π 0 ?
e `cos θ dθ.

Let us now prove Item (iv). The value α 1 p1q " 16 is obtained by a direct calculation.

In order to prove the equivalent (2.7), we first consider the case e Ñ 1, e ă 1. The change of variable s " 1{t in (5.5) yields α 2 `?1 `eq 2 " ´2 logpe ´1q as e Ñ 1, e ą 1, we infer that α 1 1 peq " ´2 logpe ´1q as e Ñ 1, e ą 1, which end the proof of (iv). Finally, Item (v) is a straightforward consequence of Items (i), (ii), (iii), (iv), and the proof of Lemma 2.2 is complete.

Aknowledgements: M. Lemou and F. Méhats acknowledge support from the CNRS project: ANR Moonrise, ANR-14-CE23-0007. A. M. Luz declares that she has no conflict of interest.

Figure 1 .

 1 Figure 1. Function α 1 (left) and its derivative α 1 1 (right).

Lemma 2 . 2 ( 2 ż θ 1

 2221 Properties of the function α 1 ). Let α 1 peq " 4 ? peq 0 pe `cos θq 1{2 dθ for e P R. (2.5)

7 0 8 0s 2 2 ż ` 8 0 ´f 7 0 8 0µ f 0 2 ż ` 8 0 2 ij

 78228780282 is r0, |Suppf 0 |s, so ż `´f 7 0 psq ´f 7 psq ¯`ds ď |Suppf 0 | psq ´f 7 psq ¯`ds ď |Suppf 0 | 2 }f ˚´f 0 } L 1 .Furthermore, for all s ě 0, we haveµ f 0 psq ď |Suppf 0 |, hence ż `psq 2 β f ˚,f 0 psqds ď |Suppf 0 | β f ˚,f 0 psqds " |Suppf 0 | pf 0 ´f ˚q`d θdv ď |Suppf 0 | 2 }f ˚´f 0 } L 1

  Proof. Recall that from the right continuity of µ f , we have(1.11), for all s ě 0, t ě 0.Finally, we state a technical lemma dealing with the case of potentials which have the special form of potentials of the HMF model. For e P R, m P R ˚and φpθq " ´m cos θ we denote

	α m peq " a φ peq " 2 ?	2	ż 2π 0	a pe `m cos θq `dθ "	?	m α 1	´e m	ānd
	introduce the angle
	$ &	0,						if e ď ´m,
	θ m peq "	arccos p´e{mq P p0, πq,	if ´m ă e ă m,
	%	π,						if e ě m.
	Therefore, for all t ě 0,
	µ f ˚φ ptq "	ffl ffl ffl ffl " pθ, vq P T ˆR : f 7	ˆaφ	ˆv2 2	*ffl ffl ffl `φpθq ˙˙ą t ffl
	"	ffl ffl ffl ffl " pθ, vq P T ˆR : a φ	ˆv2 2	*ffl ffl ffl `φpθq ˙ă µ f ptq ffl
	"	ffl ffl ffl ffl " pθ, vq P T ˆR :	v 2 2	`φpθq ă a ´1 φ pµ f ptqq	*ffl ffl ffl ffl
		´a´1
	" a φ

2.1. Let φ P C 8 pTq and let a φ be the function defined by (2.1). Let f P L 1 pT ˆRq, nonnegative. Then the function

f ˚φpθ, vq " f 7 ˆaφ ˆv2 2 `φpθq ˙˙, pθ, vq P T ˆR is equimeasurable to f , that is µ f ˚φ " µ f

, where µ f is defined by (1.10). In the sequel, the function f ˚φ is called (decreasing) rearrangement with respect to the microscopic energy v 2 2 `φpθq. φ pµ f ptqq ¯" µ f ptq.

  f 0 ptq "

		ffl ffl ffl ffl " pθ, vq P T ˆR : F	ˆv2 2	*ffl ffl ffl `φ0 pθq ˙ą t ffl
	"	ffl ffl ffl ffl " pθ, vq P T ˆR : G ˝aφ 0	ˆv2 2	*ffl ffl ffl `φ0 pθq ˙ą t ffl
	"	ffl ffl ffl ffl " pθ, vq P T ˆR : G 7 ˝aφ 0	ˆv2 2	*ffl `φ0 pθq ˙ą t ffl ffl ffl .
	Hence, applying the (1.11) to the function
	µ G psq " |tt ě 0 : Gptq ą su| , for all s ě 0,
	and to its pseudo-inverse G 7 , we get
	µ f 0 ptq "	ffl ffl ffl ffl	" pθ, vq P T ˆR : a φ 0	ˆv2 2	*ffl ffl ffl `φ0 pθq ˙ă µ G ptq ffl
		ffl
		ffl
	"	ffl
		ffl

"

pθ, vq P T ˆR :

  uq du, where gpe, uq " 4

						? 2	a	pe `uq ?1 .
								´u2
	For u P p´1, 1q, we have
	0 ď	Bg Be	ď q e puq "	2 ? 1 ? 2 ´e	a	1 pe `uqp1 ´uq	1 ´eďuď1

  `eq 2 " ´2 logp1 ´eq as e Ñ 1, e ă 1, thus α 1 1 peq " ´2 logp1 ´eq as e Ñ 1, e ă 1. To deal with the case e Ñ 1, e ą 1, we perform for e ą 1 the change of variable t " 1´cos θ 1`cos θ in (2.6):

	Using again (5.6), we get
		ˇˇα 1 1 peq ´I2 peq ˇˇď 2 ?	2	ż `8 0	p1 `tq 3{2	a	dt tptpe ´1q `e `1q	ď C 0
	where we used that tpe ´1q `e `1 ě 2 and where we set
						I 2 peq " 2 ?	2	ż `8 0	p1 `tq	dt tpe ´1q `e `1q a	.
	Since	I 2 " ´2 log	p ?	e	´1
								1 1 peq " 2	? 2	ż `8 0	ds sp1 `sqp2 `p1 ´eqsq a	.
	Let								I 1 peq " 2 ?	2	ż `8 0	p1 `sq	ds 2 `p1 ´eqs a	.
	From	0 ď	1 ? s	´1 ? 1 `s "	a	sp1 `sqp	1 ?	s `?1 `sq	ď	1 p1 `sq ?	s	,	(5.6)
	we deduce that
	ˇˇα 1 1 peq ´I1 peq ˇˇď 2	?	2	ż `8 0	p1 `sq 3{2	ds sp2 `p1 ´eqsq a	ď 2	ż `8 0	ds p1 `sq 3{2 ?	s	" C 0 .
	A direct computation yields
		I 1 peq "	´2? ? 1 `e log 2 2 `?1 α 1 1 ´e p ? 1 peq " 2 ? 2 ż `8 0	dt tp1 `tqptpe ´1q `e `1q a	.


Kévin Gosse
email: kevin.gosse@etudiant.univ-reims.fr

Stéphanie Jehan-Besson
email: stephanie.jehan@ensicaen.fr

Su Ruan
email: su.ruan@univ-rouen.fr

Comparison of 2D and 3D Region-based Deformable Models and Random Walker Methods for PET Segmentation

Keywords: Biomedical imaging, Active contours, Deformable models, level sets, PET imaging, Random walk, Poisson law, segmentation

In this paper, we propose to compare different methods for tumor segmentation in positron emission tomography (PET) images. We first propose to tackle this problem under the umbrella of shape optimization and 3D deformable models. Indeed, 2D active contours have been widely investigated in the literature but these techniques do not take advantage of 3D informations. On the one hand, we use the well-known model of Chan and Vese. On the other hand we use a criterion based on parametric probabilities which allows us to test the assumption of Poisson distribution of the intensity in such images. Both will be compared to their 2D equivalent and to an improved randomwalker algorithm. For this comparison, we use a set of simulated, phantom and real sequences with a known ground-truth and compute the corresponding Dice Coefficients. We also give some examples of 2D and 3D segmentation results.

I. INTRODUCTION

The current work is devoted to the segmentation of tumors regions in Positron Emission Tomography (PET) images. Indeed such medical images are powerful for cancer tumors visualization by detecting photons from a radiotracer localized in abnormals cells. However the conception of automatic segmentation methods in these images is difficult due to different factors and notably a high level of noise, a low contrast and unclear boundaries. The patient breathing can lead to motion artifacts and noise and there is also a large variability in shape or even in texture for the tumors to be segmented. A wide variety of image segmentation techniques have been proposed up to now and a large review of these numerous methods used to segment PET images is available in [START_REF] Foster | A review on segmentation of positron emission tomography images[END_REF].

In this paper, due to the shape variability of tumors, we rather decide to give a new insight on region-based deformable models (active contours or surfaces) for the segmentation of tumors in PET images. Active contours have been introduced by Kass, Witkin and Terzopoulos [START_REF] Kass | Snakes: Active contour models[END_REF] and later by Caselles et al [START_REF] Caselles | Geodesic active contours[END_REF] and were originally boundary methods. Due to the low contrast and unclear boundaries of PET images, we focus on region-based models where region functionals are used in addition to boundary functionals. Such powerful methods have appeared later. Introduced by [START_REF] Cohen | Surface reconstruction using active contour models[END_REF] and [START_REF] Ronfard | Region-based strategies for active contour models[END_REF], they have been further developed in [START_REF] Zhu | Region competition: Unifying snakes, region growing, and Bayes/MDL for multiband image segmentation[END_REF], [START_REF] Chan | Active contours without edges[END_REF], [START_REF] Paragios | Geodesic active regions: A new paradigm to deal with frame partition problems in computer vision[END_REF]. In [START_REF] Martin | Influence of the noise model on level set active contour segmentation[END_REF], [START_REF] Galland | Multi-component image segmentation in homogeneous regions based on description length minimization: Application to speckle, Poisson and Bernoulli noise[END_REF], [START_REF] Lecellier | Statistical region-based active contours with exponential family observations[END_REF], [START_REF] Lecellier | Region-based active contours with exponential family observations[END_REF], [START_REF] Jehan-Besson | Biomedical diagnostics and clinical technologies : Applying high performance cluster and Grid Computing, chapter Medical image segmentation and tracking through the maximisation or the minimisation of divergences between pdfs[END_REF], the authors propose to take advantage of the minimization of statistical functionals that involve parametric pdfs of the exponential family in order to better modelize the distribution of the intensity in the regions to segment. In PET images, the distribution can be assimilated to a Poisson law because it is based on the detection of Photons in the tumoral zone. In order to verify this assumption, we propose to compute the evolution equation of a deformable model that will minimize the antilog-likelihood function of a Poisson distribution. We compare it to the classical Chan and Vese method [START_REF] Chan | Active contours without edges[END_REF] which appears to be usually very robust for segmentation. We also propose to take benefit of a 3D segmentation model based on the level set method [START_REF] Osher | Fronts propagating with curvaturedependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF]. We compare the 3D segmentation to the result obtained using 2D segmentations of each slice. Moreover, we also compare these methods to two random-walker models, the classical one from [START_REF] Grady | Random walks for image segmentation[END_REF] and the Locally Adaptive Random Walk (LARW) from [START_REF] Onoma | Segmentation of heterogeneous or small FDG PET positive tissue based on a 3D-locally adaptive random walk algorithm[END_REF]. The comparison of all these methods is performed on the basis of a strong database that includes simulated images, real phantom acquisitions and clinical data. The ground truth delineations provided with this database allow us to compute the Dice Coefficient as well as sensitivity and specificity parameters in order to objectively compare the different assumptions, implementations or methods.

This paper is organized as follows: we briefly remind the principle of region-based deformable models and we explain both Chan and Vese and Poisson model in Section 2. The random-walker model used as a comparison with active contours is presented in Section 3. Experimental results are discussed in Section 4. We finally conclude and give some perspectives.

II. ACTIVE CONTOURS

Let U be a class of domains (open, regular bounded sets, i.e. C 2 ) of R n (n = 2 or 3), and Ω i an element of U of boundary ∂Ω i . The image domain is denoted by Ω I . A region-based segmentation problem aims at finding a partition of Ω I in 2 regions {Ω in , Ω out } that minimizes the following criterion:

E(Ωin, Ωout, Γ) = Jr(Ωin) + Jr(Ωout) + αJreg(Γ), (1)
where J r (Ω i ) describes the homogeneity of the region Ω i and will be further developed in sections A and B according to the underlying distribution model. The term Γ denotes the common interface between the two regions. The energy term J reg is a regularization term, balanced with a positive real parameter α, which is chosen as J reg = Γ da where da represents an area element. In 2D, this term corresponds to the minimization of the curve length and in 3D to the minimization of the surface area. This term will help us to provide a smooth curve or surface to delineate the tumor.

In order to minimize such a criterion, we propose to classically perform a shape gradient descent by computing the evolution equation of an active contour using shape derivation tools. The interested reader may refer to [START_REF] Delfour | Shapes and Geometries[END_REF] [18] [START_REF] Aubert | Image segmentation using active contours : Calculus of variations or shape gradients[END_REF] for further details in shape derivation methods and active contours based on shape gradients. Formally, from the computed shape derivatives, we can deduce the evolution equation that will drive the active contour towards a minimum of the criterion in the following general form:

∂Γ(p, τ ) ∂τ = v(x, Ω)N (x) (2)
with v(x, Ω) the velocity of the active contour (in 2D) or surface (in 3D) deduced from the shape derivatives. The velocity is directed along the unit inward normal N of ∂Ω. The active contour evolves from an initial curve Γ(τ = 0) = Γ 0 . In the equation (2), we have x = Γ(p, τ ).

Let us now detail the two main criteria that we have tested for PET segmentation and the corresponding evolution equations.

A. THE CHAN-VESE MODEL

When looking for an homogeneous region, a classic and robust model was proposed in [START_REF] Chan | Active contours without edges[END_REF] by Chan and Vese, originally inspired by the Mumford-Shah functional [START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF]. This criterion implies the use of the mean value of the intensity inside the region Ω i , denoted by µ i . We then search for the partition of the image which minimizes the square error between the intensity and the mean value of the region inside each region by using the following region descriptor J r in criterion (1):

Jr(Ωi) = Ω i (I(x) -µi(Ωi)) 2 dx with i = in or out (3)
The shape derivative of this criterion can be found in [START_REF] Jehan-Besson | DREAM 2 S: Deformable regions driven by an Eulerian accurate minimization method for image and video segmentation[END_REF] [19] and leads to the following velocity term in equation (2):

v(x, Ωin, Ωout) = (I -µin) 2 -(I -µout) 2 + ακ, ( 4
)
where κ is the curvature of Γ. This model is implicitely based on the assumption of a Gaussian model with a fixed variance and a piece-wise constant image model.

B. THE POISSON MODEL

Besides the well-known model of Chan and Vese, we propose to investigate a model based on a Poisson law. Indeed, our assumption lies on the Poisson nature of photon-counting in PET images as first investigated by [START_REF] Chan | Level set method for positron emission tomography[END_REF] for segmentation. We focus here on the minimization of the anti-log likelihood classically used by [START_REF] Zhu | Region competition: Unifying snakes, region growing, and Bayes/MDL for multiband image segmentation[END_REF], [START_REF] Goudail | Statistical image processing techniques for noisy images[END_REF] and [START_REF] Lecellier | Statistical region-based active contours with exponential family observations[END_REF], [START_REF] Lecellier | Region-based active contours with exponential family observations[END_REF]. In that case, the region based term is the following :

Jr(Ωi) = - Ω i log p(I(x), λi)dx, (5)
with p the Poisson distribution of the intensity I inside the region Ω i indexed with the parameter λ i (i = in or out), such that:

p(I(x), λ i ) = λ I(x) i e -λi I(x)! .
We can calculate the maximum likelihood estimate of the parameter λ i as follows:

λi M V = 1 |Ωi| Ω i I(x)dx. (6)
The log-likelihood is then:

log p(I(x), λi M V = I(x) log λi M V -λi M V -log(I(x)!),
which leads to the following region-based term:

Jr(Ωi) = λi M V |Ωi| + Ω i log(I(x)!)dx -log λi M V Ω i I(x)dx. (7)
We then propose to compute the shape derivative of the previous region-based term which leads to the following theorem:

Theorem 1: The shape derivative in the direction V of the functional (5), with p the pdf of Poisson distribution and λ iM V the parameter of this distribution estimated by the maximum likelihood estimate, is

< J r (Ω), V >= - ∂Ω vi(x, Ω) < V, N > da where v i (x, Ω) = -I(x) log( λ iM V ) + λ iM V + log(I(x)!) .
The complete proof is given in Appendix. Note that a more general proof for shape derivatives of statistical region-based terms using pdfs of the exponential family can be found in [START_REF] Lecellier | Region-based active contours with exponential family observations[END_REF]. The proof for the Poisson law is here given for completeness. Using the fact that the two regions Ω in and Ω out share the same interface with opposite normals, we can get our velocity term v(x, Ω in , Ω out ) for equation (2) :

v(x, Ωin, Ωout) = -I log( λin M V ) + λin M V (8)
+ I log( λout M V ) -λout M V + ακ.

C. IMPLEMENTATION OF THE ACTIVE CONTOUR

As far as the numerical implementation is concerned, we use the level set method introduced in [START_REF] Osher | Fronts propagating with curvaturedependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF]. The key idea is to introduce an auxiliary function U (x, τ ) such that Γ(τ ) is the zero level set of U . The function U is often chosen to be the signed distance function of Γ(τ ). The evolution equation then becomes:

∂U (x, τ ) ∂τ = v(x, Ωin, Ωout)|∇U (x, τ )|. ( 9
)
This method is accurate and allows to automatically handle the topological changes of the initial curve. It can be used in 2D and 3D. For the 3D implementation, κ is chosen as the mean curvature. The 3D implementation is interesting in the case of PET sequences since it allows to introduce automatically a global 3D regularization term. In the experimental results, we show that the 3D segmentation allows a better segmentation of tumors than the 2D segmentation applied to each slice separately.

III. THE RANDOM-WALKER MODEL

In this section, we first explain the classical random walk algorithm [START_REF] Grady | Random walks for image segmentation[END_REF] and then its extension for PET segmentation known as the LARW (Locally Adaptive Random Walk) method [START_REF] Onoma | Segmentation of heterogeneous or small FDG PET positive tissue based on a 3D-locally adaptive random walk algorithm[END_REF].

A. CLASSICAL RANDOM WALK (RW)

In this model, an image is considered like a graph G = (V, E). Each voxel is associated with a node of the graph and edges link adjacent voxels. This method requires a small number, called seeds, of labeled voxels as an initial step. The key idea behind this method is to imagine a random walker starting from each unlabeled voxels. This random-walker has a certain probability to first arrive at a specific seed. The label is then chosen as the same label as the seed with the maximum probability of arrival of the random-walker. A common weight used for edges is:

ωij = exp -β(gi -gj) 2 , ( 10
)
where both g i and g j are the image intensity values at points i and j. The term β denotes a free parameter depending on the application. We define the Laplacian matrix of the graph by

Lij =      di if i = j, -ωij if vi et vj are adjacent nodes, 0 otherwise.
with d i = ω ij for each node. This matrix can be rearranged with labeled and unlabeled voxels as follows:

L = LM B B T LU ,
where M and U correspond to labeled and unlabeled pixels. It has been shown in [START_REF] Grady | Random walks for image segmentation[END_REF] that the solution of the random-walk for a specific label minimizes the following energy:

D[x] = 1 2 x T Lx = 1 2 e ij ∈E ωij(xi -xj) 2 ,
which leads to the linear system:

LU X = -B T M. (11)

B. LOCALLY ADAPTIVE RANDOM WALK (LARW)

An improved version of the classical random-walker algorithm was proposed in [START_REF] Onoma | Segmentation of heterogeneous or small FDG PET positive tissue based on a 3D-locally adaptive random walk algorithm[END_REF]. The proposed segmentation method initializes automatically seeds in tumor voxels using Fuzzy-C Means (FCM), and then delineates the tumor volume using the improved RW. They made two major changes to the algorithm described above. First, they changed the weight of edges to take care of the spatial distance between voxels. Equation ( 10) is then modified to:

ωij = exp -β(gi -gj) 2 hij ( 12
)
where the h ij term is the spatial distance between voxels i and j.

They also proposed to strengthen the grouping of voxels having similar intensity by adding the likelihood of probability to each class (tumor and non tumor). Supposing each class k having a Gaussian distribution, their observation likelihood function can be written as:

p k (I, µ k , σ k ) = 1 σ k √ 2π exp -(I -µ k ) 2 2σ 2 k ( 13
)
where µ k and σ k are respectively the mean and the variance of the class k. The matrix B T is then modified

B = p k .B T
which leads to the new linear system to solve

LU X = -B M (14)
the complete algorithm is detailed in [START_REF] Onoma | Segmentation of heterogeneous or small FDG PET positive tissue based on a 3D-locally adaptive random walk algorithm[END_REF].

IV. EXPERIMENTAL RESULTS

A. DATA SET AND METHOD OF EVALUATION

In order to test the different algorithms, we use the training data set provided by the MICCAI PETSEG Challenge [START_REF]PET segmentation challenge using a data management and processing infrastructure[END_REF]. Following recommendations by the TG2111 [START_REF] Hatt | Classification and evaluation strategies of auto-segmentation approaches for PET report of the american association of physicists in medicine task group 211[END_REF], this evaluation dataset consists of a mixture of numerical simulations with realistic uptake distributions and a variety of tumor shapes, physical phantoms acquisitions, and real clinical images [START_REF] Shepherd | Design of a benchmark platform for evaluating PET-based contouring accuracy in oncology applications[END_REF], [START_REF] Berthon | Towards a standard for the evaluation of PET autosegmentation methods: requirements and implementation[END_REF]. Simulated images consist of realistic tumors embedded within numerical phantoms, with corresponding SIMSET [START_REF] Aristophanous | The development and testing of a digital PET phantom for the evaluation of tumor volume segmentation techniques[END_REF] or GATE [START_REF] Lemaitre | Incorporating patient specific variability in the simulation of realistic whole body 18F-FDG distributions for oncology applications[END_REF], [START_REF] Papadimitroulas | Investigation of realistic PET simulations incorporating tumor patient's specificity using anthropomorphic models: creation of an oncology database[END_REF] simulated data to better account for the physics of PET acquisition. Voxel-by-voxel ground-truth is available from the simulations. The physical phantoms are zeolites with different sizes and shapes incorporated in an anthropomorphic phantom, with repeated acquisitions of the same zeolites [START_REF] Zito | The use of zeolites to generate PET phantoms for the validation of quantification strategies in oncology[END_REF]. The surrogate of truth for these images is the corresponding high resolution CT, thresholded as to obtain the exact known volume of the zeolite. Clinical images are lung [START_REF] Wanet | Gradient-based delineation of the primary GTV on FDG-PET in non-small cell lung cancer: a comparison with threshold-based approaches, ct and surgical specimens[END_REF], [START_REF] Lapuyade-Lahorgue | Speqtacle: An automated generalized fuzzy C-means algorithm for tumor delineation in PET[END_REF] and head and neck [START_REF] Geets | A gradient-based method for segmenting FDG-PET images: methodology and validation[END_REF] tumors, with either histopathology-based contours [START_REF] Wanet | Gradient-based delineation of the primary GTV on FDG-PET in non-small cell lung cancer: a comparison with threshold-based approaches, ct and surgical specimens[END_REF], [START_REF] Geets | A gradient-based method for segmenting FDG-PET images: methodology and validation[END_REF] or a consensus of several manual delineations by experts [START_REF] Lapuyade-Lahorgue | Speqtacle: An automated generalized fuzzy C-means algorithm for tumor delineation in PET[END_REF] obtained with STAPLE [START_REF] Warfield | Simultaneous truth and performance level estimation (STAPLE): an algorithm for the validation of image segmentation[END_REF] to serve as a surrogate of truth.

Image sequences are then divided into three parts. The different images resolutions are detailed in Table 1. The different segmentation results from the above algorithms were compared to the provided ground-truth using the Dice Similarity Coefficient which measures the extent of spatial overlap between two binary sequences X and Y and is expressed as follows:

DC = 2|X ∩ Y | |X| + |Y | ( 15
)
where |X| stands for the area of X.

We also use the sensitivity and the specificity coefficients. Sensitivity, also called the true positive rate (TPR), is the proportion of the tumor correctly identified by the algorithm:

B. 2D AND 3D COMPARISON

For this experiment, we perform 2D and 3D segmentations using the level set approach detailed in section II with both the Chan-Vese and the Poisson models. Then, we compute the Dice Similarity Coefficient for each slice of each sequence which allows us to compare 2D and 3D segmentations for a same model. In Fig. 1, we trace the obtained DC values for each segmented slice of the first sequence using the Chan-Vese model. The tumor is located between the slices 5 to 20 which corresponds to the numbers 1 to 16 in the Fig. 1 . Experimental results, not reported here for space reasons, are very similar for other sequences and Poisson model. It clearly appears that considering the whole 3D information improves segmentation results especially at the beginning and the end of the 3D sequences when tumors tend to be less visible and very small. In the example given in Fig. 2, we can see that the 2D segmentation shows a lot of false positives. Indeed, it takes into account the lower left of the image which is of similar intensity, but not part of the tumor. The effect is clearly less visible in 3D.

C. METHODS COMPARISON

For this second experiment, we propose to compare the different methods presented in this paper. We limit ourselves with 3D segmentation since we have seen in the previous section that 3D segmentation outperforms 2D segmentation in all cases. For the random-walker method, we make two versions of the algorithm. The first one is the classical randomwalker algorithm [START_REF] Grady | Random walks for image segmentation[END_REF] described briefly in section III. The second one is the Locally Adaptive Random-Walk version from [START_REF] Onoma | Segmentation of heterogeneous or small FDG PET positive tissue based on a 3D-locally adaptive random walk algorithm[END_REF] with minor changes (auto-initialization of seed's background like the one's from the tumor and unsystematic median-filtering).

First of all, in Fig. 3 and Fig. 4, we give two examples of the different segmentation results and corresponding contours obtained on two extracted slices (one from the simulated set and one from the clinical set). For each method, we compare the mask to the ground truth (expert segmentation) and we show in red the false positive points and in green the false negative points. We can observe that both methods are able to segment accurately the tumor from the simulated set (Fig. 3). In the case of more complicated clinical data with non homogeneous tumors (see Fig. 4), the different automatic methods are less accurate. In this case the Poisson model leads to less false negative but more false positive than the Chan and Vese methods. For such images, a pre-processing such as denoising or contrast enhancement could be valuable. Let us now consider the computation of the DC (Dice Coefficient) between the different methods and the ground truth for all the sequences of the training data-set (Table 2).

Considering these results it appears that the assumption of Poisson distribution in PET images is relatively true especially for simulated images (sequences 1 to 6). In the case of real phantom acquisitions (sequences 7 to 15), the method is less robust than the Chan and Vese method. Indeed the piece-wise constant model assumed by the Chan-Vese method seems more robust in this case due to the fact that the phantom has a constant intensity. As far as clinical data are concerned, the Chan-Vese model and Poisson model seem quite equivalent. The Chan and Vese method is less accurate in some cases (sequences 16 and 18) due to the low contrast and more accurate in sequences 17 and 19. In fact the real distribution model of the intensity in these images is probably a mixture between a Poisson law and a Gaussian model due to the reconstruction process and also to the artifacts encountered in acquisition methods (patient breathing, low resolution, camera noise model). Globally, for all the different sequences, the Chan and Vese method is more robust.

When computing the specificity and sensitivity coefficients (see Table 3), we can observe that the Poisson method offers very good sensitivity coefficients but lowest specificity values. On the contrary, the LARW method shows very good specificity values with lower sensitivity values. The Chan and Vese method offers a good compromise between sensitivity and specificity coefficients.

The 3D Chan-Vese and LARW methods are the most robust algorithms and can perform great segmentation in many kind of images. We can observe that Chan-Vese method tends to sub-segment tumors and it is worth noting that Chan-Vese segmentation is often a subset of Poisson segmentation. An interesting fact for RW and LARW is the importance of the initialization of seeds. The more the better is not necessary true here. If we have more information on where the background is exactly, and made no change in tumor's seeds, it will necessarily sub-segment the previous segmentation because non-labeled voxels have more probability to reach background then. In a general way, we have to choose correctly the amount of seeds for background and tumor. Here, we choose them in order to obtain the best results. More detailed experimental results can be found in [START_REF] Onoma | Segmentation of heterogeneous or small FDG PET positive tissue based on a 3D-locally adaptive random walk algorithm[END_REF].

V. CONCLUSION AND PERSPECTIVES

In this paper, we propose to investigate the use of 2D and 3D deformable models for PET images segmentation. These deformable models are currently based on the minimization of statistical criteria and are implemented using the level set method. In one of the tested criterion, we take into account the Poisson distribution of PET images. Our results show a great improvement of the segmentation from 2D to 3D and validate the implication of Poisson law in PET images. However, the Chan-Vese model (implicitely based on the assumption of a Gaussian model) still remains the most robust algorithm and gives comparable results to the LARW method [START_REF] Onoma | Segmentation of heterogeneous or small FDG PET positive tissue based on a 3D-locally adaptive random walk algorithm[END_REF]. Our ongoing research is now directed towards two main challenges. First of all, the use of an improved statistical model based on a mixture of a Poisson and a Gaussian law could be valuable. Secondly, an intelligent fusion of the different segmentation results may also be investigated in order to improve the segmentation results.

VI. ACKNOWLEDGEMENT

We would like to thank the organizing committee of the MICCAI PETSEG Challenge [START_REF]PET segmentation challenge using a data management and processing infrastructure[END_REF] for their permission to publish our evaluation results using their training database.

APPENDIX PROOF OF THEOREM 1

In the case of Poisson distribution, the region-based term that we want to derive according to the domain Ω is the following: We have to compute the shape derivative of J r (Ω). For this purpose, we remind here the following fundamental theorem which establishes a relation between the shape derivative of region-based terms of the form J r (Ω) = Ω k(x, Ω)dx:

Jr(Ω) = -
Theorem 2: The shape derivative < J r (Ω), V > of the functional J r (Ω) = Ω k(x, Ω)dx in the direction V is

< J r (Ω), V >= Ω k (x, Ω, V )dx - ∂Ω k(x, Ω) < V • N > da
where k is the shape derivative of k. The term N denotes the unit inward normal to ∂Ω and da its area element (in R 2 , we have da = ds where s stands for the arc length). Using this theorem and the equation of the Maximum Likelihood estimate for the Poisson parameter λ given in (6), we find:

1 )

 1 Sequences 1 to 6 : simulated images (GATE and SIM-SET). 2) Sequences 7 to 15 : real phantom acquisitions 3) Sequences 16 to 19 : patient clinical data

  Sensitivity =Number of true positives Size of the tumor[START_REF] Onoma | Segmentation of heterogeneous or small FDG PET positive tissue based on a 3D-locally adaptive random walk algorithm[END_REF] Specificity, also called the true negative rate (TNR), is the proportion of non-tumor correctly identified as such: Specificity = Number of true negatives Size of the non-tumor[START_REF] Delfour | Shapes and Geometries[END_REF]

Fig. 1 .

 1 Fig. 1. 2D and 3D segmentation comparison on a simulated image sequence (sequence 1). The graph represents the values of the DC coefficient computed for each 2D slice (from 1 to 16) between the expert contour and the resulting contour of the 2D segmentation method (in blue) and the 3D one (in red).

Fig. 2 .

 2 Fig. 2. a. PET image, b. Ground-Truth, c. Initial contour, d. Final contour in red for 2D Chan-Vese segmentation (smoothing coefficient α = 1), e. 2D segmentation comparison with ground-truth (false positive in red and false negative in green), f. 3D segmentation comparison with ground-truth (smoothing coefficient α = 1, false positive in red and false negative in green).

  Fig. 3. a. PET image from the simulated set with the different contours (in yellow the ground truth, in blue the 3D Chan and Vese segmentation, in red the 3D Poisson segmentation and in green the LARW segmentation). We also display the masks of the ground truth in b., of the Chan and Vese segmentation in c., of the Poisson segmentation in d. and of the LARW segmentation in e. For the different methods, the false positives are shown in red and the false negative in green.

  Fig. 4. a. PET image from the clinical set with the different contours (in yellow the ground truth, in blue the 3D Chan and Vese segmentation, in red the 3D Poisson segmentation and in green the LARW segmentation). We also display the masks of the ground truth in b., of the Chan and Vese segmentation in c., of the Poisson segmentation in d. and of the LARW segmentation in e. For the different methods, the false positives are shown in red and the false negative in green.

Ω

  log p(I(x), λ)dx= λMV |Ω| + Ω log(I(x)!)dx -log λMV Ω I(x)dx

+

  <J r (Ω), V > = -∂Ω log(I(x)!) < V • N > da + |Ω| < λMV , V > -λMV ∂Ω < V • N > da -< λ M V , V > |Ω| 1 -I |Ω| λMVwhereI = Ω I(x)dx and < λ M V , V >= -1 |Ω| ∂Ω (I(x)λ M V ) < V • N > da. Since 1 -I |Ω| λ M V= 0, we can then conclude that :< J r (Ω), V >= ∂Ω log p(I(x), λMV ) < V • N > dawhere p(I(x), λ M V ) is the Poisson pdf model.

Table 1 .

 1 Resolutions of the different sequences of the PETSEG data set

	Sequence	Name	Dimensions
	1	Simu 2 TEP 4	27x35x23
	2	Simu 2 TEP 1	15x15x15
	3	Simu 2 TEP 2	17x18x15
	4	Simu 2 TEP 3	17x23x18
	5	Simu 1 PET2	32x26x34
	6	Simu 1 PET1	32x26x34
	7	phantom 1 PET2	22x17x17
	8	phantom 1 PET3	22x17x17
	9	phantom 1 PET1	22x17x17
	10	phantom 2 PET3	13x14x12
	11	phantom 2 PET1	13x14x12
	12	phantom 2 PET2	13x14x12
	13	phantom 3 PET3	21x21x13
	14	phantom 3 PET2	21x21x13
	15	phantom 3 PET1	21x21x13
	16	clinical 3 PET	28x23x28
	17	clinical 2 4mm PET	13x14x14
	18	clinical 2 2mm PET	25x25x23
	19	clinical 1 PET	19x18x15

Table 2 .

 2 Comparison of Dice coefficients of the different methods.

	Sequence	Chan-Vese 3D Poisson 3D	RW	LARW
	1	0.9040	0.8489	0.7783	0.8679
	2	0.9448	0.9441	0.9061	0.8764
	3	0.8674	0.8766	0.7813	0.8121
	4	0.8974	0.9764	0.8357	0.8017
	5	0.8789	0.8970	0.7491	0.8702
	6	0.9012	0.8924	0.7283	0.8747
	7	0.8221	0.7791	0.8070	0.7747
	8	0.8082	0.7824	0.8030	0.7445
	9	0.8025	0.7864	0.8057	0.7461
	10	0.7600	0.5112	0.8105	0.7972
	11	0.7614	0.4908	0.7122	0.8133
	12	0.7415	0.5229	0.7513	0.7778
	13	0.8235	0.5100	0.8192	0.8435
	14	0.7958	0.4946	0.8029	0.8584
	15	0.8014	0.5122	0.8346	0.8509
	16	0.8492	0.8792	0.8516	0.8469
	17	0.7630	0.7323	0.7301	0.6742
	18	0.7629	0.7688	0.6907	0.7171
	19	0.6444	0.5714	0.5570	0.6743
	Mean	0.8158	0.7251	0.7766	0.8014

Table 3 .

 3 Comparison of Sensitivity (TPR) and Specificity (TNR) parameters.

		Chan-Vese 3D	Poisson 3D	LARW
	Sequence	TPR	TNR	TPR	TNR	TPR	TNR
	1	0.8431	0.9927	0.9770	0.8929	0.7967	0.9870
	2	0.9032	0.9990	0.9912	0.9878	0.7801	1
	3	0.8607	0.9770	0.9694	0.9551	0.8036	0.9675
	4	0.8139	1	0.9581	0.9992	0.6696	0.9998
	5	0.8053	0.9965	0.9268	0.9823	0.8609	0.9851
	6	0.8750	0.9915	0.9670	0.9746	0.8901	0.9816
	7	0.7428	0.9930	0.9357	0.9494	0.6608	0.9951
	8	0.7283	0.9920	0.9164	0.9538	0.6254	0.9941
	9	0.7154	0.9927	0.9084	0.9564	0.6190	0.9956
	10	0.95	0.9791	1	0.9273	0.7125	0.9971
	11	0.9375	0.98	1	0.9211	0.7625	0.9957
	12	0.95	0.9767	1	0.9306	0.7	0.9962
	13	0.9739	0.9920	1	0.9607	0.8435	0.9968
	14	0.9826	0.99	1	0.9582	0.8435	0.9975
	15	0.9826	0.9904	1	0.9610	0.8435	0.972
	16	0.7551	0.9949	0.9625	0.9490	0.8071	0.9778
	17	0.7523	0.9799	0.8692	0.9537	0.5607	0.9906
	18	0.7178	0.9765	0.8442	0.9495	0.6596	0.9741
	19	0.7973	0.9466	0.8640	0.9085	0.7867	0.9569
	Mean	0.9460	0.9863	0.9521	0.9511	0.7487	0.9874