Stéphane Bessy

Fedor V Fomin

Serge Gaspers
email: sgaspers@dim.uchile.cl

Christophe Paul

Anthony Perez

Saket Saurabh

Stéphan Thomassé

Kernels for Feedback Arc Set In Tournaments

A tournament T = (V, A) is a directed graph in which there is exactly one arc between every pair of distinct vertices. Given a digraph on n vertices and an integer parameter k, the Feedback Arc Set problem asks whether the given digraph has a set of k arcs whose removal results in an acyclic digraph. The Feedback Arc Set problem restricted to tournaments is known as the k-Feedback Arc Set in Tournaments (k-FAST) problem. In this paper we obtain a linear vertex kernel for k-FAST. That is, we give a polynomial time algorithm which given an input instance T to k-FAST obtains an equivalent instance T ′ on O(k) vertices. In fact, given any fixed ǫ > 0, the kernelized instance has at most (2 + ǫ)k vertices. Our result improves the previous known bound of O(k 2 ) on the kernel size for k-FAST. Our kernelization algorithm solves the problem on a subclass of tournaments in polynomial time and uses a known polynomial time approximation scheme for k-FAST.

Introduction

Given a directed graph G = (V, A) on n vertices and an integer parameter k, the Feedback Arc Set problem asks whether the given digraph has a set of k arcs whose removal results in an acyclic directed graph. In this paper, we consider this problem in a special class of directed graphs, tournaments. A tournament T = (V, A) is a directed graph in which there is exactly one directed arc between every pair of vertices. More formally the problem we consider is defined as follows.

k-Feedback Arc Set in Tournaments (k-FAST): Given a tournament T = (V, A) and a positive integer k, does there exist a subset F ⊆ A of at most k arcs whose removal makes T acyclic.

In the weighted version of k-FAST, we are also given integer weights (each weight is at least one) on the arcs and the objective is to find a feedback arc set of weight at most k. This problem is called k-Weighted Feedback Arc Set in Tournaments (k-WFAST).

Feedback arc sets in tournaments are well studied from the combinatorial [START_REF] Erdös | On sets on consistent arcs in tournaments[END_REF][START_REF] Jung | On subgraphs without cycles in tournaments[END_REF][START_REF] Reid | Disproof of a conjecture of Erdös and Moser on tournaments[END_REF][START_REF] Seshu | Linear Graphs and Electrical Networks[END_REF][START_REF] Spencer | Optimal ranking of tournaments[END_REF][START_REF] Younger | Minimum feedback arc sets for a directed graph[END_REF], statistical [START_REF] Slater | Inconsistencies in a schedule of paired comparisons[END_REF] and algorithmic [START_REF] Ailon | Aggregating inconsistent information: ranking and clustering[END_REF][START_REF] Alon | Ranking tournaments[END_REF][START_REF] Coppersmith | Ordering by weighted number of wins gives a good ranking for weighted tournaments[END_REF][START_REF] Kenyon-Mathieu | How to rank with few errors[END_REF][START_REF] Van Zuylen | Deterministic approximation algorithms for ranking and clusterings[END_REF][START_REF] Van Zuylen | Deterministic pivoting algorithms for constrained ranking and clustering problems[END_REF] points of view. The problems k-FAST and k-WFAST have several applications. In rank aggregation we are given several rankings of a set of objects, and we wish to produce a single ranking that on average is as consistent as possible with the given ones, according to some chosen measure of consistency. This problem has been studied in the context of voting [START_REF] Borda | Mémoire sur les élections au scrutin[END_REF][START_REF] Condorcet | Essai sur l'application de l'analyse à la probabilité des décisions rendues à la pluralité des voix[END_REF], machine learning [START_REF] Cohen | Learning to order things[END_REF], and search engine ranking [START_REF] Dwork | Rank aggregation revisited[END_REF][START_REF] Dwork | Rank aggregation methods for the web[END_REF]. A natural consistency measure for rank aggregation is the number of pairs that occur in a different order in the two rankings. This leads to Kemeny rank aggregation [START_REF] Kemeny | Mathematics without numbers[END_REF][START_REF] Kemeny | Mathematical models in the social sciences[END_REF], a special case of k-WFAST.

The k-FAST problem is known to be NP-complete by recent results of Alon [START_REF] Alon | Ranking tournaments[END_REF] and Charbit et al. [START_REF] Charbit | The minimum feedback arc set problem is NP-hard for tournaments[END_REF] while k-WFAST is known to be NP-complete by Bartholdi III et al. [START_REF] Bartholdi | Voting schemes for which it can be difficult to tell who won the election[END_REF]. From an approximation perspective, k-WFAST is APX-hard [START_REF] Speckenmeyer | On feedback problems in digraphs[END_REF] but admits a polynomial time approximation scheme when the edge weights are bounded by a constant [START_REF] Kenyon-Mathieu | How to rank with few errors[END_REF]. The problem is also well studied in parameterized complexity. In this area, a problem with input size n and a parameter k is said to be fixed parameter tractable (FPT) if there exists an algorithm to solve this problem in time f (k)•n O (1) , where f is an arbitrary function of k. Raman and Saurabh [START_REF] Raman | Parameterized algorithms for feedback set problems and their duals in tournaments[END_REF] showed that k-FAST and k-WFAST are FPT by obtaining an algorithm running in time O(2.415 k • k 4.752 + n O (1) ). Recently, Alon et al. [START_REF] Alon | Fast FAST[END_REF] have improved this result by giving an algorithm for k-WFAST running in time 1) ). This algorithm runs in sub-exponential time, a trait uncommon to parameterized algorithms. In this paper we investigate k-FAST from the view point of kernelization, currently one of the most active subfields of parameterized algorithms.

O(2 O( √ k log 2 k) + n O(
A parameterized problem is said to admit a polynomial kernel if there is a polynomial (in n) time algorithm, called a kernelization algorithm, that reduces the input instance to an instance whose size is bounded by a polynomial p(k) in k, while preserving the answer. This reduced instance is called a p(k) kernel for the problem. When p(k) is a linear function of k then the corresponding kernel is a linear kernel. Kernelization has been at the forefront of research in parameterized complexity in the last couple of years, leading to various new polynomial kernels as well as tools to show that several problems do not have a polynomial kernel under some complexity-theoretic assumptions [START_REF] Bodlaender | On problems without polynomial kernels (extended abstract)[END_REF][START_REF] Bodlaender | kernelization[END_REF][START_REF] Bousquet | A polynomial kernel for multicut in trees[END_REF][START_REF] Dom | Incompressibility through colors and IDs[END_REF][START_REF] Thomassé | A quadratic kernel for feedback vertex set[END_REF]. In this paper we continue the current theme of research on kernelization and obtain a linear vertex kernel for k-FAST. That is, we give a polynomial time algorithm which given an input instance T to k-FAST obtains an equivalent instance T ′ on O(k) vertices. More precisely, given any fixed ǫ > 0, we find a kernel with a most (2 + ǫ)k vertices in polynomial time. The reason we call it a linear vertex kernel is that, even though the number of vertices in the reduced instance is at most O(k), the number of arcs is still O(k 2 ). Our result improves the previous known bound of O(k 2 ) on the vertex kernel size for k-FAST [START_REF] Alon | Fast FAST[END_REF][START_REF] Dom | Fixed-parameter tractability results for feedback set problems in tournaments[END_REF]. For our kernelization algorithm we find a subclass of tournaments where one can find a minimum sized feedback arc set in polynomial time (see Lemma 3.8) and use the known polynomial time approximation scheme for k-FAST by Kenyon-Mathieu and Schudy [START_REF] Kenyon-Mathieu | How to rank with few errors[END_REF]. The polynomial time algorithm for a subclass of tournaments could be of independent interest.

The paper is organized as follows. In Section 2, we give some definition and preliminary results regarding feedback arc sets. In Section 3 we give a linear vertex kernel for k-FAST. Finally we conclude with some remarks in Section 4.

Preliminaries

Let T = (V, A) be a tournament on n vertices. We use T σ = (V σ , A) to denote a tournament whose vertices are ordered under a fixed ordering σ = v 1 , . . . , v n (we also use D σ for an ordered directed graph). We say that an arc v i v j of T σ is a backward arc if i > j, otherwise we call it a forward arc. Moreover, given any partition P := {V 1 , . . . , V l } of V σ , where every V i is an interval according to the ordering of T σ , we use A B to denote all arcs between the intervals (having their endpoints in different intervals), and A I for all arcs within the intervals. If T σ contains no backward arc, then we say that it is transitive.

For a vertex v ∈ V we denote its in-neighborhood by

N -(v) := {u ∈ V | uv ∈ A} and its out-neighborhood by N + (v) := {u ∈ V | vu ∈ A}. A set of vertices M ⊆ V is a module if and only if N + (u) \ M = N + (v) \ M for every u, v ∈ M .
For a subset of arcs A ′ ⊆ A, we define T [A ′ ] to be the digraph (V ′ , A ′ ) where V ′ is the union of endpoints of the arcs in A ′ . Given an ordered digraph D σ and an arc e = v i v j , S(e) = {v i , . . . , v j } denotes the span of e. The number of vertices in S(e) is called the length of e and is denoted by l(e). Thus, for every arc e = v i v j , l(e) = |i -j| + 1. Finally, for every vertex v in the span of e, we say that e is above v.

In this paper, we will use the well-known fact that every acyclic tournament admits a transitive ordering. In particular, we will consider maximal transitive modules. We also need the following result for our kernelization algorithm.

Lemma 2.1. ( [START_REF] Raman | Parameterized algorithms for feedback set problems and their duals in tournaments[END_REF]) Let D = (V, A) be a directed graph and F be a minimal feedback arc set of D. Let D ′ be the graph obtained from D by reversing the arcs of F in D, then D ′ is acyclic.

In this paper whenever we say circuit, we mean a directed cycle. Next we introduce a definition which is useful for a lemma we prove later. Definition 2.2. Let D σ = (V σ , A) be an ordered directed graph and let f = vu be a backward arc of D σ . We call certificate of f , and denote it by c(f ), any directed path from u to v using only forward arcs in the span of f in D σ .

Observe that such a directed path together with the backward arc f forms a directed cycle in D σ whose only backward arc is f . Definition 2.3. Let D σ = (V σ , A) be an ordered directed graph, and let F ⊆ A be a set of backward arcs of D σ . We say that we can certify F whenever it is possible to find a set F = {c(f ) : f ∈ F } of arc-disjoint certificates for the arcs in F .

Let D σ = (V σ , A) be an ordered directed graph, and let F ⊆ A be a subset of backward arcs of D σ . We say that we can certify the set F using only arcs from A ′ ⊆ A if F can be certified by a collection F such that the union of the arcs of the certificates in F is contained in A ′ . In the following, f as(D) denotes the size of a minimum feedback arc set, that is, the cardinality of a minimum sized set F of arcs whose removal makes D acyclic. Lemma 2.4. Let D σ be an ordered directed graph, and let P = {V 1 , . . . , V l } be a partition of D σ into intervals. Assume that the set F of all backward arcs of D σ [A B ] can be certified using only arcs from

A B . Then f as(D σ ) = f as(D σ [A I ]) + f as(D σ [A B ]).
Moreover, there exists a minimum sized feedback arc set of D σ containing F .

Proof. For any bipartition of the arc set

A into A 1 and A 2 , f as(D σ ) ≥ f as(D σ [A 1 ]) + f as(D σ [A 2 ]
). Hence, in particular for a partition of the arc set A into A I and A B we have that f as

(D σ ) ≥ f as(D σ [A I ]) + f as(D σ [A B ]). Next, we show that f as(D σ ) ≤ f as(D σ [A I ]) + f as(D σ [A B ]
). This follows from the fact that once we reverse all the arcs in F , each remaining circuit lies in D σ [V i ] for some i ∈ {1, . . . , l}. In other words once we reverse all the arcs in F , every circuit is completely contained in D σ [A I ]. This concludes the proof of the first part of the lemma. In fact, what we have shown is that there exists a minimum sized feedback arc set of D σ containing F . This concludes the proof of the lemma.

Kernels for k-FAST

In this section we first give a subquadratic vertex kernel of size O(k √ k) for k-FAST and then improve on it to get our final vertex kernel of size O(k). We start by giving a few reduction rules that will be needed to bound the size of the kernels. Rule 3.1. If a vertex v is not contained in any triangle, delete v from T . Rule 3.2. If there exists an arc uv that belongs to more than k distinct triangles, then reverse uv and decrease k by 1.

We say that a reduction rule is sound, if whenever the rule is applied to an instance (T, k) to obtain an instance (T ′ , k ′ ), T has a feedback arc set of size at most k if and only if T ′ has a feedback arc set of size at most k ′ . Lemma 3.1. ( [START_REF] Alon | Fast FAST[END_REF][START_REF] Dom | Fixed-parameter tractability results for feedback set problems in tournaments[END_REF]) Rules 3.1 and 3.2 are sound and can be applied in polynomial time.

The Rules 3.1 and 3.2 together led to a quadratic kernel for k-WFAST [START_REF] Alon | Fast FAST[END_REF]. Earlier, these rules were used by Dom et al. [START_REF] Dom | Fixed-parameter tractability results for feedback set problems in tournaments[END_REF] to obtain a quadratic kernel for k-FAST. We now add a new reduction rule that will allow us to obtain the claimed bound on the kernel sizes for k-FAST. Given an ordered tournament T σ = (V σ , A), we say that P = {V 1 , . . . , V l } is a safe partition of V σ into intervals whenever it is possible to certify the backward arcs of T σ [A B ] using only arcs from A B . Rule 3.3. Let T σ be an ordered tournament, P = {V 1 , . . . , V l } be a safe partition of V σ into intervals and F be the set of backward arcs of T σ [A B ]. Then reverse all the arcs of F and decrease k by |F |. Lemma 3.2. Rule 3.3 is sound.

Proof. Let P be a safe partition of T σ . Observe that it is possible to certify all the backward arcs, that is F , using only arcs in A B . Hence using Lemma 2.4 we have that f as(T σ ) = f as(T σ [A I ]) + f as(T σ [A B ]). Furthermore, by Lemma 2.4 we also know that there exists a minimum sized feedback arc set of D σ containing F . Thus, T σ has a feedback arc set of size at most k if and only if the tournament T ′ σ obtained from T σ by reversing all the arcs of F has a feedback arc set of size at most k -|F |.

A subquadratic kernel for k-FAST

In this section, we show how to obtain an O(k √ k) sized vertex kernel for k-FAST. To do so, we introduce the following reduction rule. Rule 3.4. Let V m be a maximal transitive module of size p, and I and O be the set of in-neighbors and out-neighbors of the vertices of V m in T , respectively. Let Z be the set of arcs uv such that u ∈ O and v ∈ I. If q = |Z| < p then reverse all the arcs in Z and decrease k by q.

I O Proof. We first prove that the partition P = {I, V m , O} forms a safe partition of the input tournament. Let V ′ m = {w 1 , . . . , w q } ⊆ V m be an arbitrary subset of size q of V m and let

T m |T m | = p q < p
Z = {u i v i | 1 ≤ i ≤ q}. Consider the collection F = {v i w i u i | u i v i ∈ Z, w i ∈ V ′
m } and notice that it certifies all the arcs in Z. In fact we have managed to certify all the backwards arcs of the partition using only arcs from A B and hence P forms a safe partition. Thus, by Rule 3.3, it is safe to reverse all the arcs from O to I. The time complexity follows from the fact that computing a modular decomposition tree can be done in O(n + m) time on directed graphs [START_REF] Mcconnell | Linear-time modular decomposition of directed graphs[END_REF].

We show that any Yes-instance to which none of the Rules 3.1, 3. Proof. Let S be a feedback arc set of size at most k of T and let T ′ be the tournament obtained from T by reversing all the arcs in S. Let σ be the transitive ordering of T ′ and T σ = (V σ , A) be the ordered tournament corresponding to the ordering σ. We say that a vertex is affected if it is incident to some arc in S. Thus, the number of affected vertices is at most 2|S| ≤ 2k. The reduction Rule 3.1 ensures that the first and last vertex of T σ are affected. To see this note that if the first vertex in V σ is not affected then it is a source vertex (vertex with in-degree 0) and hence it is not part of any triangle and thus Rule 3.1 would have applied. We can similarly argue for the last vertex. Next we argue that there is no backward arc e of length greater than 2k + 2 in T σ . Assume to the contrary that e = uv is a backward arc with S(e) = {v, x 1 , x 2 , . . . , x 2k+1 , . . . , u} and hence l(e) > 2k +2. Consider the collection T = {vx i u | 1 ≤ i ≤ 2k} and observe that at most k of these triples can contain an arc from S \{e} and hence there exist at least k +1 triplets in T which corresponds to distinct triangles all containing e. But then e would have been reversed by an application of Rule 3.2. Hence, we have shown that there is no backward arc e of length greater than 2k + 2 in T σ . Thus e∈S l(e) ≤ 2k 2 + 2k. We also know that between two consecutive affected vertices there is exactly one maximal transitive module. Let us denote by t i the number of vertices in these modules, where i ∈ {1, . . . , 2k -1}. The objective here is to bound the number of vertices in V σ or V using 2k-1 i=1 t i . To do so, observe that since T is reduced under the Rule 3.4, there are at least t i backward arcs above every module with t i vertices, each of length at least t i . This implies that 2k-1 i=1 t 2 i ≤ e∈S l(e) ≤ 2k 2 + 2k. Now, using the Cauchy-Schwarz inequality we can show the following.

2k-1 i=1 t i = 2k-1 i=1 t i • 1 ≤ 2k-1 i=1 t 2 i • 2k-1 i=1 1 ≤ (2k 2 + 2k) • (2k -1) = 4k 3 + 2k 2 -k.
Thus every reduced Yes-instance has at most

√ 4k 3 + 2k 2 -k + 2k = O(k √ k) vertices.

A linear kernel for k-FAST

We begin this subsection by showing some general properties about tournaments which will be useful in obtaining a linear kernel for k-FAST.

Backward Weighted Tournaments

Let T σ be an ordered tournament with weights on its backward arcs. We call such a tournament a backward weighted tournament and denote it by T ω , and use ω(e) to denote the weight of a backward arc e. For every interval I := [v i , . . . , v j ] we use ω(I) to denote the total weight of all backward arcs having both their endpoints in I, that is, ω(I) = e=uv w(e) where u, v ∈ I and e is a backward arc.

Definition 3.5. (Contraction) Let T ω = (V σ , A) be an ordered tournament with weights on its backward arcs and I = [v i , . . . , v j ] be an interval. The contracted tournament is defined as

T ω ′ = (V σ ′ = V σ \ {I} ∪ {c I }, A ′ ).
The arc set A ′ is defined as follows.

• It contains all the arcs

A 1 = {uv | uv ∈ A, u / ∈ I, v / ∈ I} • Add A 2 = {uc I | uv ∈ A, u / ∈ I, v ∈ I} and A 3 = {c I v | uv ∈ A, u ∈ I, v / ∈ I}.
• Finally, we remove every forward arc involved in a 2-cycle after the addition of arcs in the previous step.

The order σ ′ for T ω ′ is provided by σ ′ = v 1 , . . . , v i-1 , c I , v j+1 , . . . , v n . We define the weight of a backward arc e = xy of A ′ as follows.

w ′ (xy) =    w(xy) if xy ∈ A 1 {xz∈A | z∈I} w(xz) if xy ∈ A 2 {zy∈A | z∈I} w(zy) if xy ∈ A 3
We refer to Figure 2 for an illustration. Next we generalize the notions of certificate and certification (Definitions 2.2 and 2.3) to backward weighted tournaments. Definition 3.6. Let T ω = (V σ , A) be a backward weighted tournament, and let f = vu ∈ A be a backward arc of T ω . We call ω-certificate of f , and denote it by C(f ), a collection of ω(f ) arc-disjoint directed paths going from u to v and using only forward arcs in the span of f in T ω . Definition 3.7. Let T ω = (V σ , A) be a backward weighted tournament, and let F ⊆ A be a subset of backward arcs of T ω . We say that we can ω-certify F whenever it is possible to find a set F = {C(f ) : f ∈ F } of arc-disjoint ω-certificates for the arcs in F . Lemma 3.8. Let T ω = (V σ , A) be a backward weighted tournament such that for every interval I := [v i , . . . , v j ] the following holds:

v p v p c I 2 v j v i
2 • ω(I) ≤ |I| -1 (1)
Then it is possible to ω-certify the backward arcs of T ω .

Proof.

Let V σ = v 1 , . . . , v n .
The proof is by induction on n, the number of vertices. Note that by applying (1) to the interval I = [v 1 , . . . , v n ], we have that there exists a vertex v i in T ω that is not incident to any backward arc. Let T ′ ω = (V ′ σ , A ′ ) denote the tournament T ω \ {v i }. We say that an interval I is critical whenever |I| ≥ 2 and 2 • ω(I) = |I| -1. We now consider several cases, based on different types of critical intervals.

(i) Suppose that there are no critical intervals. Thus, in T ′ ω , every interval satisfies (1), and hence by induction on n the result holds.

(ii) Suppose now that the only critical interval is I = [v 1 , . . . , v n ], and let e = vu be a backward arc above v i with the maximum length. Note that since v i does not belong to any backward arc, we can use it to form a directed path c(e) = uv i v, which is a certificate for e. We now consider T ′ ω where the weight of e has been decreased by 1. In this process if ω(e) becomes 0 then we reverse the arc e. We now show that every interval of T ′ ω respects [START_REF] Ailon | Aggregating inconsistent information: ranking and clustering[END_REF]. If an interval I ′ ∈ T ′ ω does not contain v i in the corresponding interval in T ω , then by our assumption we have that 2 • ω(I ′ ) ≤ |I ′ | -1. Now we assume that the interval corresponding to

I ′ in T ω contains v i but either u / ∈ I ′ ∪ {v i } or v / ∈ I ′ ∪ {v i }. Then we have 2 • ω(I ′ ) = 2 • ω(I) < |I| -1 = |I ′ | and hence we get that 2 • ω(I ′ ) ≤ |I ′ | -1.
Finally, we assume that the interval corresponding to

I ′ in T ω contains v i and u, v ∈ I ′ ∪ {v i }. In this case, 2 • ω(I ′ ) = 2 • (ω(I) -1) ≤ |I| -1 -2 < |I ′ | -1.
Thus, by the induction hypothesis, we obtain a family of arc-disjoint ω-certificates F ′ which ω-certify the backward arcs of T ′ ω . Observe that the maximality of l(e) ensures that if e is reversed then it will not be used in any ω-certificate of F ′ , thus implying that F ′ ∪ c(e) is a family ω-certifying the backward arcs of T ω .

(iii) Finally, suppose that there exists a critical interval I V σ . Roughly speaking, we will show that I and V σ \ I can be certified separately. To do so, we first show the following.

Claim. Let I ⊂ V σ be a critical interval. Then the tournament T ω ′ = (V σ ′ , A ′ ) obtained from T ω by contracting I satisfies the conditions of the lemma.

Proof. Let H ′ be any interval of T ω ′ . As before if H ′ does not contain c I then the result holds by hypothesis. Otherwise, let H be the interval corresponding to H ′ in T ω . We will show that 2ω(H ′ ) ≤ |H ′ | -1. By hypothesis, we know that 2ω(H) ≤ |H| -1 and that 2ω(I) = |I| -1. Thus we have the following.

2ω(H ′ ) = 2 • (ω(H) -ω(I)) ≤ |H| -1 -|I| + 1 = (|H| + 1 -|I|) -1 = |H ′ | -1
Thus, we have shown that the tournament T ω ′ satisfies the conditions of the lemma.

We now consider a minimal critical interval I. By induction, and using the claim, we know that we can obtain a family of arc-disjoint ω-certificates F ′ which ωcertifies the backward arcs of T ω ′ without using any arc within I. Now, by minimality of I, we can use (ii) to obtain a family of arc-disjoint ω-certificates F ′′ which ω-certifies the backward arcs of I using only arcs within I. Thus, F ′ ∪ F ′′ is a family ω-certifying all backward arcs of T ω .

This concludes the proof of the lemma.

In the following, any interval that does not respect condition ( 1) is said to be a dense interval. Lemma 3.9. Let T ω = (V σ , A) be a backward weighted tournament with |V σ | ≥ 2p + 1 and ω(V σ ) ≤ p. Then there exists a safe partition of V σ with at least one backward arc between the intervals and it can be computed in polynomial time.

Proof. The proof is by induction on n = |V σ |. Observe that the statement is true for n = 3, which is our base case.

For the inductive step, we assume first that there is no dense interval in T ω . In this case Lemma 3.8 ensures that the partition of V σ into singletons of vertices is a safe partition. So from now on we assume that there exists at least one dense interval.

Let I be a dense interval. By definition of I, we have that ω(I) ≥ 1 2 • |I|. We now contract I and obtain the backward weighted tournament T ω ′ = (V σ ′ , A ′ ). In the contracted tournament T ω ′ , we have:

|V σ ′ | ≥ 2p + 1 -(|I| -1) = 2p -|I| + 2; ω ′ (V σ ′ ) ≤ p -1 2 • |I|. Thus, if we set r := p -1 2 • |I|, we get that |V σ ′ | ≥ 2r + 1 and ω ′ (V σ ′ ) ≤ r. Since |V σ ′ | < |V σ |,
by the induction hypothesis we can find a safe partition P of T ω ′ , and thus obtain a family F ω ′ that ω-certifies the backward arcs of T ω ′ [A B ] using only arcs in A B .

We claim that P ′ obtained from P by substituting c I by its corresponding interval I is a safe partition in T ω . To see this, first observe that if c I has not been used to ω-certify the backward arcs in T ω ′ [A B ], that is, c I is not an end point of any arc in the ω-certificates, then we are done. So from now on we assume that c I has been part of a ω-certificate for some backward arc. Let e be a backward arc in T ω ′ [A B ], and let c ω ′ (e) ∈ F ω ′ be a ω-certificate of e. First we assume that c I is not the first vertex of the certificate c ω ′ (e) (with respect to ordering σ ′ ), and let c 1 and c 2 be the left (in-) and right (out-) neighbors of c I in c ω ′ (e). By definition of the contraction step together with the fact that there is a forward arc between c 1 and c I and between c I and c 2 in T ω ′ , we have that there were no backward arcs between any vertex in the interval corresponding to c I and c 1 and c 2 in the original tournament T ω . So we can always find a vertex in I to replace c I in c ω ′ (e), thus obtaining a certificate c(e) for e in T ω [A B ] (observe that e remains a backward arc even in T ω ). Now we assume that c I is either a first or last vertex in the certificate c ω ′ (e). Let e ′ be an arc corresponding to e in T ω ′ with one of its endpoints being e I ∈ I. To certify e ′ in T ω [A B ], we need to show that we can construct a certificate c(e ′ ) using only arcs of T ω [A B ]. We have two cases to deal with. (ii) If c I is the last vertex of c ω ′ (e) then let c q be its left neighbor in c ω ′ (e). Once again, we have that there are only forward arcs between c q and vertices in I, and thus between c q and e I . So using this we can construct a ω-certificate for e ′ in T ω .

Notice that the fact that all ω-certificates are pairwise arc-disjoint in T ω ′ [A B ] implies that the corresponding ω-certificates are arc-disjoint in T ω [A B ], and so P ′ is indeed a safe partition of V σ .

We are now ready to give the linear size kernel for k-FAST. To do so, we make use of the fact that there exists a polynomial time approximation scheme for this problem [START_REF] Kenyon-Mathieu | How to rank with few errors[END_REF]. Theorem 3.10. For every fixed ǫ > 0, there exists a vertex kernel for k-FAST with at most (2 + ǫ)k vertices that can be computed in polynomial time.

Proof. Let (T = (V, A), k) be an instance of k-FAST. For a fixed ǫ > 0, we start by computing a feedback arc set S of size at most (1 + ǫ 2 )k. To find such a set S, we use the known polynomial time approximation scheme for k-FAST [START_REF] Kenyon-Mathieu | How to rank with few errors[END_REF]. Then, we order T with the transitive ordering of the tournament obtained by reversing every arc of S in T . Let T σ denote the resulting ordered tournament. By the upper bound on the size of S, we know that T σ has at most (1 + ǫ 2 )k backward arcs. Thus, if T σ has more than (2 + ǫ)k vertices then Lemma 3.9 ensures that we can find a safe partition with at least one backward arc between the intervals in polynomial time. Hence we can reduce the tournament by applying Rule 3.3. We then apply Rule 3.1, and repeat the previous steps until we do not find a safe partition or k = 0. In the former case, we know by Lemma 3.9 that T can have at most (2 + ǫ)k vertices, thus implying the result. In all other cases we return No.

Conclusion

In this paper we obtained linear vertex kernel for k-FAST, in fact, a vertex kernel of size (2 + ǫ)k for any fixed ǫ > 0. The new bound on the kernel size improves the previous known bound of O(k 2 ) on the vertex kernel size for k-FAST given in [START_REF] Alon | Fast FAST[END_REF][START_REF] Dom | Fixed-parameter tractability results for feedback set problems in tournaments[END_REF]. It would be interesting to see if one can obtain kernels for other problems using either polynomial time approximation schemes or a constant factor approximation algorithm for the corresponding problem. An interesting problem which remains unanswered is, whether there exists a linear or even a o(k 2 ) vertex kernel for the k-Feedback Vertex Set in Tournaments (k-FVST) problem. In the k-FVST problem we are given a tournament T and a positive integer k and the aim is to find a set of at most k vertices whose deletion makes the input tournament acyclic. The smallest known kernel for k-FVST has size O(k 2 ).

Figure 1 :Lemma 3 . 3 .

 133 Figure 1: A transitive module on which Rule 3.4 applies.

Theorem 3 . 4 .

 34 2 and 3.4 could be applied has at most O(k √ k) vertices. Let (T = (V, A), k) be a Yes-instance to k-FAST which has been reduced according to Rules 3.1, 3.2 and 3.4. Then T has at most O(k √ k) vertices.

Figure 2 :

 2 Figure 2: Illustration of the contraction step for the interval I := [v i , . . . , v j ].

  (i) If c I is the first vertex of c ω ′ (e) then let c 1 be its right neighbor in c ω ′ (e). Using the same argument as before, there are only forward arcs between any vertex in I and c 1 . In particular, there is a forward arc e I c 1 in T ω , meaning that we can construct a ω-certificate for e ′ in T ω by setting c(e ′ ) := (c ω ′ (e) \ {c I }) ∪ {e I }.

Figure 3 :

 3 Figure 3: On the left, the ω-certificate c ω ′ (e) ∈ F ω ′ . On the right, the corresponding ω-certificate obtained in T ω by replacing c I by the interval I.


Oscar E Ramos
email: oramos@laas.fr

Nicolas Mansard
email: nmansard@laas.fr

Olivier Stasse
email: ostasse@laas.fr

Philippe Souères
email: soueres@laas.fr

Walking on Non-planar Surfaces using an Inverse Dynamic Stack of Tasks

This paper presents a method to handle walking on non-planar surfaces. The method obtains the trajectory of the center of mass and the next position of the foot from a pattern generator. Then, an inverse dynamics control scheme with a quadratic programming optimization solver is used to let the foot go from its initial to its final position, controlling also the center of mass and the waist. This solver is able to handle an arbitrary number of contact points. When the swinging foot is going down, collision points are detected and they are added as contact points to the model as soon as they appear. If there are three or more contact points, the foot can safely step, but if there are one or two contact points, the foot rotates properly to generate the largest support polygon. Using this heuristic, the foot can stand on non-planar surfaces. The results show the simulation of HRP-2 walking on a surface with obstacles.

I. INTRODUCTION

Biped walking is a sequencial task that alternates the motion of both legs and that strongly relies on contact forces between the feet and the ground in order to realize stable motions. Static balance keeps the center of mass (CoM) inside the support polygon, whereas dynamic balance keeps the zero moment point (ZMP) inside the support polygon allowing the CoM to be outside of it for a certain period of time. Dynamic balance, which represents the way that humans walk, is the preferred walking strategy in robotics research but classical schemes allow mainly walking only on planar surfaces.

However, there are several works that have succeded to walk on uneven terrain. In [START_REF] Yamaguchi | Development of a biped walking robot adapting to a horizontally uneven surface[END_REF], a control method for adapting to a horizontally composed plane with unknown steps height is proposed. The method allows walking on non-flat floor by using the sole sensors that measure a distance between the sole itself and the ground before landing. In [START_REF] Nishiwaki | High frequency walking pattern generation based on preview control of zmp[END_REF], a pattern generation system that updates the CoM trajectory in a short cycle using preview control and a feedack controller that tracks the ZMP reference was proposed. This method also enabled a robot to change the walking pattern according to the change of applied force. A method that uses preview control from the current inclination of the upper body to a stable future state in order to update the upper body trajectory was proposed in [START_REF] Wieber | Trajectory free linear model predictive control for stable walking in the presence of strong perturbations[END_REF]. Another method [START_REF] Hashimoto | Landing pattern modification method with predictive attitude and compliance control to deal with uneven terrain[END_REF] is based on a predictive attitude compensation control and a nonlinear compliance control using only force sensors and a particular biped foot system that adapts to a rough terrain.A technique to stabilize the upper body postion under varying ground reaction forces is proposed in [START_REF] Takenaka | Real time motion generation and control for biped robot-4th report: Integrated balance control[END_REF]. This technique generates moments to handle Fig. 1: The robot HRP-2 stepping over a non-planar obstacle with its right foot. disturbances that cannot be handled by ground reaction forces by using horizontal and rotational acceleration of the upper body and step duration change: as a result, a landing position is locally modified.

More recently, a strategy for adjusting the ZMP reference to generate online repetitive walking has been proposed [START_REF] Nishiwaki | Strategies for adjusting the zmp reference trajectory for maintaining balance in humanoid walking[END_REF]. The strategy consists on changing the reference ZMP inside the sole area, changing the position of the next step and changing the duration of the current step. In [START_REF] Morisawa | Reactive biped walking control for a collision of a swinging foot on uneven terrain[END_REF], proper impedance gains of the feet and updates of the desired landing position of the CoM pattern after a contact is detected are used as a method for a collision of a swinging foot on uneven terrain.

In this paper, we propose to use a generic hierarchical optimization solver to track the position of the CoM and to autonomously control the swinging foot based on the next foot position on the ground given by a pattern generator. The considered solver consists of the inverse-dynamics control cascade proposed in [START_REF] Saab | Generation of Dynamic Motion for Anthropomorphic Systems under Prioritized Equality and Inequality Constraints[END_REF], [START_REF] Saab | Generic Dynamic Motion Generation with Multiple Unilateral Constraints[END_REF], and modified in [START_REF] Mansard | A dedicated solver for fast operational-space inverse dynamics[END_REF] to allow faster computations. As the foot goes down, possible collisions are detected and handled by the solver as contact points, and, if necessary, the foot is properly rotated to move towards the ground without losing dynamic stability. Figure 1 shows an example of an obstacle on which the robot can stand keeping its stability.

This paper is organized as follows. Section II presents the inverse dynamics control including the contact constraints. Section III gives a brief summary of the pattern generator used in this work. The usage of the inverse dynamics scheme with the pattern generator is presented in Section IV. Finally, the results of the simulation are presented in Section V.

II. INVERSE DYNAMICS CONTROL

A. Generic Dynamic Model

The humanoid robot configuration will be represented by the generalized coordinates q = [x b q a ] T , where x b represents the position and orientation of the robot (more specifically, the base, free-floating or free-flyer) with respect to the world, and q a represents the n actuated joints. It will be assumed that there are n f points in contact with the environment. Let

x p = [x 1 x 2 • • • x n f ] T ∈ R 3n f represent the contact points in the world frame, and f = [f 1 f 2 • • • f n f ] T ∈ R 3n f
represent the punctual contact forces, with f i acting at point x i . The dynamic model of the robot is

S T τ = Aq + b + J T f f (1)
where S = [0 I] is a matrix that selects the actuated joints, τ is the torque vector on the joints, A is the whole-body inertia matrix, b is the drift including Coriolis, centrifugal and gravity forces, and J f = ∂xp ∂q is the Jacobian of the contact points. Additionally, for rigid contacts, the contact points must not move and the restriction ẍp = 0 must be added. This restriction is equivalent to

J f q + Jf q = 0 (2)
provided that ẋp = J f q. Let S n be a matrix that selects the normal component (third component) of each force in f and of each point in x p , with its elements given by {s n ij } = δ(3i -j), where δ represents the Kronecker delta function. When considering the model of "point contact" [START_REF] Featherstone | Rigid Body Dynamics Algorithms[END_REF], the entire contact constraint is described simultaneously by S n x p = 0 and S n f ≥ 0, if the contact persists. The first part of the constraint is already satisfied by [START_REF] Nishiwaki | High frequency walking pattern generation based on preview control of zmp[END_REF], and the second part relating the contact forces

S n f ≥ 0 (3)
must be explicitly considered in the model. The unidirectionallity of (3) guarantees that the ZMP will exist inside the support polygon, ensuring the dynamic stability condition [START_REF] Saab | Generic Dynamic Motion Generation with Multiple Unilateral Constraints[END_REF].

B. Decoupled Dynamics

The explicit separation of the dynamic equations into motion and actuation has been proposed in [START_REF] Mansard | A dedicated solver for fast operational-space inverse dynamics[END_REF]. The decoupling allows a faster computation since the decoupled variables are of lower dimension than the coupled original variables. This section will give a brief review of the decoupled model. Using the Choleski decomposition, the inverse of the inertia matrix A can be expressed as A -1 = BB T , and thus A = B -T B -1 , with B an invertible triangular matrix. Multiplying the dynamic model (1) by B T on the left we obtain:

B T S T τ = B -1 q + B T b + G T f (4)
where G = J f B. The result of multiplying (4) by G on the left, using (2) to guarantee no kinematic motion at the rigid contact points, is:

GB T S T τ = GB T b + GG T f -Jf q (5)
and solving for B T S T τ gives

B T S T τ = B T b + G T f + δ c + V u (6)
where δ c = -G + Jf q is the contact drift, V is a basis of the kernel of G such that V V T = I -G + G, and u is a vector in this kernel. The matrix S = [I 0] is defined so that it cancels out S T ( SS T = 0), and thus the torque. By multiplying both terms of ( 6) by SB -T on the left, a torquefree relation between u and the contact forces f is obtained:

0 = Sb + SJ T f f + SB -T δ c + SB -T V u (7)
This expression is composed of both ( 1) and ( 2), but explicitly using u and f as variables. The other contact condition related to the unidirectionallity of the forces (3) can be directly included in the decoupled model as S n f ≥ 0.

It is also possible to recover the torques from this model multiplying ( 6) by SB -T on the left to obtain:

τ = Sb + SJ T f f + SB -T (δ c + V u) (8)
where the fact that SS T = I has been used.

C. Representation of the tasks

For any observable s i (q) of the robot configuration, a desired position and orientation s * i can be specified without loss of generality by the i th task e i = s i -s * i . A usual choice for the reference acceleration of the task is the proportional derivative (PD) control law ë * i = -K p e i -K v ėi with K p > 0 and K v = 2 K p , to make the error decrease exponentially to zero. At the acceleration level, the classical relation between the task space and the joint space becomes:

ëi = Ji q + J i q (9)
where J i = ∂ei ∂q is the i th task Jacobian. The equivalence between ( 4) and ( 6) leads to the direct relation between u and the acceleration q as:

q = Bδ c + BV u (10)
so that u can be interpreted as a "motion" variable. Then, using [START_REF] Mansard | A dedicated solver for fast operational-space inverse dynamics[END_REF], the generic task representation (9) can be equivalently expressed in terms of the decoupled variable u as:

ëi = Ji q + J i Bδ c + J i BV u (11)

D. Spatial-force reduction

Let x c ∈ R 6 represent the position and orientation of a central point rigidly attached to the foot (for instance, its center of mass, or the ankle). The 6D Jacobian associated with this frame is J c = ∂xc ∂q , and the relation between J c and the Jacobian J f in (1) is J f = XJ c , with [START_REF] Mansard | A dedicated solver for fast operational-space inverse dynamics[END_REF]:

X =    I -p1 . . . . . . I -p n f    (12)
where p i is the position of the contact point i with respect to the frame x c , and pi represents the skew-symmetric matrix obtained from p i equivalent to a cross product pre-multiplication. When there are three or more contact points on the same body, X is full-column rank and the null space of J f is equal to the null space of J c . In this case, J c can be used to compute the V basis as V = ker(J c B) instead of V = ker(J f B), and the contact drift as δ c = -(J c B) + Jc q. However, the null spaces are not equal when there are only 1 or 2 contact points, which is a case that will be discussed in section IV-B

E. QP for the dynamic inversion

The robot motion implies the fulfillment of one or more tasks satisfying at the same time the dynamic model and the contact constraints. Under this premise, the optimization problem for a task i is expressed using ( 11), [START_REF] Morisawa | Reactive biped walking control for a collision of a swinging foot on uneven terrain[END_REF], and (3) as

min u,f J i BV u + Ji q + J i Bδ c -ë * i s.t. ( SB -T V )u + ( SJ T f )f = -S(b + B -T δ c ) S n f ≥ 0
In [START_REF] Saab | Generation of Dynamic Motion for Anthropomorphic Systems under Prioritized Equality and Inequality Constraints[END_REF], it has been proposed to use a hierarchical QP solver (HQP) [START_REF] Escande | Fast resolution of hierarchized inverse kinematics with inequality constraints[END_REF] to solve the operational-space inverse dynamics. The hierarchy is denoted by ≺ where a ≺ b means that a has more priority than b. With this notation, the dynamic stack of m tasks is (7) ≺ (3) ≺ (11.1) ≺ • • • ≺ (11.m), and the optimization variables determined by the HQP solver are the decoupled u and f . Using these optimization variables, the torque τ obtained from (8) can be used as input to a torquecontrolled robot or simulator, and the acceleration q from (10) can be integrated to be used as input to a position-controlled robot.

III. PATTERN GENERATOR

The pattern generator assumes that the robot CoM is kept at a certain constant distance z cm above the ground, so that the position of the CoM can be represented only by two components (x cm , y cm ). For simplicity in notation,

x cm k = x cm (kT ) = x cm (t k ),
where T is the sampling period and k is the k-th sample, and the state variables are

xk = [x cm k ẋcm k ẍcm k ] T .
The ZMP position on the ground plane is represented as (z x , z y ). In the following treatment only the x components will be explicited, but the y components are obtained in a similar way. The discrete dynamic system relating the ZMP and the CoM is [START_REF] Kajita | Biped walking pattern generation by using preview control of zero-moment point[END_REF]:

  x cm k+1 ẋcm k+1 ẍcm k+1   =   1 T T 2 /2 0 1 T 0 0 1     x cm k ẋcm k ẍcm k   +   T 3 /6 T 2 /2 T   ... x cm k (13) z x = 1 0 -z cm /g   x cm k ẋcm k ẍcm k   (14)
Considering a prediction horizon of N samples, and using the previous dynamics recursively, it can be shown [START_REF] Wieber | Trajectory free linear model predictive control for stable walking in the presence of strong perturbations[END_REF] that the velocity of the CoM from time t k+1 to t k+N is expressed as:

Ẋk+1 = ẋcm k+1 • • • ẋcm k+N T = P vs xk + P vu ... Xk (15)
where ...

Xk = [ ... x cm k • • • ... x cm k+N -1 ]
T is the jerk of the CoM from time t k to t k+N -1 , and

P vs =    0 1 T . . . . . . . . . 0 1 N T    , P vu =    T 2 /2 0 0 . . . . . . 0 (1 + 2N )T 2 /2 • • • T 2 /2   
Also, the ZMP from time t k+1 to time t k+N is expressed as:

Z x k+1 = [z x k+1 • • • z x k+N ] T = P zs xk + P zu ... Xk (16)
where

P zs =    1 T T 2 /2 -z cm /g . . . . . . . . . 1 N T N 2 T 2 /2 -z cm /g    P zu =    T 3 /6 -T z cm /g 0 0 . . . . . . . . . (1 + 3M + 3M 2 ) T 3 6 -T z cm g • • • T 3 6 -T z cm g    with M = N -1.
The pattern generator used in this work is the one proposed in [START_REF] Herdt | Online walking motion generation with automatic footstep placement[END_REF], which basically regulates the speed of the CoM and obtains the foot placement as output of the optimization process. The considered optimization variable is

u k = [ ... Xk ... Y k X f k Y f k ] T , where (X f k , Y f k )
corresponds to the positions on the ground of the following m foot steps. The optimization problem in [START_REF] Herdt | Online walking motion generation with automatic footstep placement[END_REF] is stated as:

min u k α 2 Ẋk+1 -Ẋref k+1 2 + Ẏk+1 -Ẏ ref k+1 2 + β 2 Z x k+1 -Z x ref k+1 2 + Z y k+1 -Z y ref k+1 2 + γ 2 ... Xk 2 + ... Y k 2 ( 17
)
where Ẋref k+1 and Ẏ ref k+1 are the desired mean value for the speed of the CoM, and Z x ref k+1 , Z y ref k+1 are the references for the ZMP. These ZMP references are not fixed in advanced but are permanently recomputed from the feet position decided by the algorithm so that the ZMP lies in the middle of the foot. Let (X f c k , Y f c k ) be the current position of the foot on the ground. Then, the ZMP references are given by:

Z x ref k+1 = U c k+1 X f c k + U k+1 X f k (18) Z y ref k+1 = U c k+1 Y f c k + U k+1 Y f k ( 19
)
with

U c k+1 =                  1 . . . 1 0 . . . 0 0 . . . 0                  U k+1 =                  0 0 0 • • • 0 . . . . . . . . . . . . . . . 0 0 0 • • • 0 1 0 0 • • • 0 . . . . . . . . . . . . . . . 1 0 0 • • • 0 0 1 0 • • • 0 . . . . . . . . . . . . . . . 0 1 0 • • • 0                 
Additionally, during the single support phase, the constraint on the ZMP to ensure that it lies inside the support polygon is expressed as:

d x (θ) d y (θ) z x -x f z y -y f ≤ b(θ) (20)
where (x f , y f ) is the position of the foot, θ is its orientation, d x (θ), d y (θ) are column vectors containing the x, y coordinates of the normal vectors to the edges of the feet, and b(θ) is the column vector containing their position with a security margin.

IV. USING THE PATTERN GENERATOR WITH THE INVERSE DYNAMIC CONTROL SCHEME

The pattern generator gives the desired trajectory of the CoM and of the swinging foot. Classical schemes track these elements using only kinematics and no notion of contact with the ground. The method proposed here uses the inverse dynamics stack of tasks control scheme described in section II with the online output of the pattern generator. While the CoM and the wais are tracked as tasks, the swinging foot is not tracked. Since the solver can handle the notion of contact with a surface (the ground or a local irregularity), only the next position of the swinging foot is considered. Then, an interpolation task to take the foot from its initial position to its final position is used. This task is interrupted as soon as a contact appears leading to a possible rotation of the foot to reach stability.

A. Interpolation task

The interpolation task is used to take the swinging foot from its initial position to its final position. The objective is to achieve a desired final position and velocity in a specific period of time, assuming that the current position and velocity are known. Let the evolution of a generic feature be described by x(t), and let the initial time be t 0 with corresponding initial position and velocity given by x 0 = x(t 0 ) and ẋ0 = ẋ(t 0 ), respectively. After a desired time duration T , the desired position and velocity will be x f = x(t 0 + T ) and ẋf = ẋ(t 0 + T ). The control is performed via the acceleration ẍ(t) which is set to a linear function of the form [START_REF] De Lasa | Feature-Based Locomotion Controllers[END_REF]: whose solution leads to quadratic velocity and cubic position trajectories given by:

ẍ(t) = ẍ0 + ẍf -ẍ0 T (t -t 0 ) (21)
ẋ(t) = ẍf -ẍ0 2T (t -t 0 ) 2 + ẍ0 (t -t 0 ) + ẋ0 (22) x(t) = ẍf -ẍ0 6T (t -t 0 ) 3 + ẍ0 2 (t -t 0 ) 2 + ẋ0 (t -t 0 ) + x 0 (23)
The unknown parameters ẍ0 and ẍf can be obtained evaluating ( 22) and ( 23) at time t f . However, as time passes by, the initial position and velocity as well as the time remaining for t f are updated, and a new acceleration trajectory ( 21) is computed. Then, since ( 21) is permanently updated, only the initial acceleration ẍ0 is needed (accelerations at the following instants of time are of no interest), and the dynamic task in [START_REF] Featherstone | Rigid Body Dynamics Algorithms[END_REF] is simply defined as ë * = ẍ0 , which can be obtained by solving the system ( 22) and ( 23) at time t f resulting:

ë * = ẍ0 = 6 T 2 (x f -x 0 ) - 2 T ( ẋf + 2 ẋ0 ) (24)
It should be noted that T becomes smaller after each iteration and is eventually close to zero as the initial time reaches the final time. To handle this situation, when T < t h , where t h is a threshold, the previous value of ẍ0 is kept, which in practice gives good results.

B. Handling 1 and 2 contact points

One of the main objectives of this control scheme is to explicitly handle the case of a rigid body with only 2 or 1 contact point(s) with the environment. In this situation, X in ( 12) is not full-column rank, and the reductions where J c could be used instead of J f do not necessarily hold since both matrices do not have the same null space anymore.

In this case, the equivalence J f = XJ c has to be used and the computation of V will be explicitly ker(XJ c B).

The relation between the derivative of both Jacobians is

Jf = X Jc + ẊJ c , but the second term is null provided that Ẋ is null (the contact points do not move) leading to Jf = X Jc . Then, for the contact drift, δ c = -(XJ c B) + X Jc q.

C. Walking generation scheme

In order to make the robot walk, the output of the pattern generator is used as input to the inverse dynamics solver as Fig. 3: A single contact point on the sole of the robot.

shown in Figure 2. The dynamic solver (section II-E) considers the following three tasks:

• Task to track the trajectory of the CoM. The pattern generator gives the desired trajectory for the center of mass assuming a constant height. This task controls only the x and y components of the CoM.

• Task to partially track the waist trajectory. This task aims to control the height (position in z) of the waist at a certain constant value (since this is the basic premise of the pattern generator), as well as the orientation in x, y. It is partial because not all the six degrees of freedom are controlled.

• Interpolation tasks on the swinging foot. The pattern generator gives the next final position of the foot on the ground. This task takes the swinging foot from its initial position to its desired final position. There are two interpolation steps: the first task takes the foot from the initial position to an intermediate position that lies halfway between the initial and final positions with a predetermined height. The second interpolation takes the foot from the intermediate to the final position, if possible. If a contact is detected before arriving to the final position, the collision is handled as described in the following section. The interpolation task (section IV-A) guarantees that the time constraint is fulfilled according to the pattern generator.

D. Non-horizontal foot contacts handling

The second interpolation for the foot, described in the previous section, would ideally end without interruption if the ground is flat. However, if there is a non-planar surface, generated by the own irregular ground or by a small obstacle (whose size is assumed to be smaller than the step height), the interpolation will not end but will be interrupted as soon as some point(s) on the foot sole touch(es) a rigid object. In this case, there are three possibilities for the contact situation:

• There is only one contact point • There are two contact points • There are three or more contact points If there are more than three contact points, the foot can step on those points, which will generate the support polygon. When there are less than 3 contact points, the foot still needs to move to find at least one more contact. In this case, if the foot is left uncontrolled after the first contact(s), the robot dynamics might take the foot to an unstable position. In order to avoid these instabilities, the foot extremes have to be controlled so that the maximum support polygon is obtained.

If there is only one contact point, the foot is in the situation shown by figure 3, where the single contact point is p c 1 and the foot extremes are denoted by p e 1 , p e 2 , p e 3 and p e 4 . Four triangles are formed by joining the contact point to two consecutive foot extremes. The area of the triangle formed by the consecutive foot extremes p e i and p e i+1 is given by

A i = 0.5 (p c 1 -p e i ) × (p c 1 -p e i+1 ) (25)
The triangle with the greatest area will contain the extremes of the foot that are farther from the contact point, and thus, it is desirable to take those extremes to the ground so that the largest support polygon is obtained. To that end, a task is assigned to each of these extremes controlling only the vertical z position so that they go to the ground. While these points are moving, if a contact point is detected, it is added as a contact to the solver and its position is checked against the points that were going to the ground. The task for the closest extreme is removed. The foot continues its rotation until another contact point is detected, in which case, the remaining extreme task is removed.

If there are only two contact points, there are two possibilities shown in figure 4. To obtain the case, a line L is passed through both contact points and the extremes of the foot are determined to lie in one side or on the other side of the line. Let the contact points be p c 1 and p c 2 , the vector joining these points be v = (v x , v y ) = p c 1 -p c 2 , and the vector joining one point with one extreme of the foot be v i = (v ix , v iy ) = p e i -p c 2 . The idea is to rotate the points so that the line L is aligned with the vertical line. Then, the sign of the arctangent can be used to determine the side of the line in which a point lies. The angle that line L must rotate to be aligned with the vertical is θ = atan2(v y , v x ). Then, each v i is rotated by the angle θ as: After this rotation, the angle to the line is determined as φ = atan2(v f iy , v f ix ), the sign of φ indicating whether the point is on one side or the other of line L. If only one point is on one side and three are on the other side (Figure 4a), the three points are taken to the ground. If two points are on each side (Figure 4b), then, the area of the quadrilateral formed by the contact points and the extremes on each side is computed and the largest area indicates that those extremes are farther and must be taken to the ground. As in the case of a single contact, as soon as a new contact appears, it is added as a contact to the solver, and the task for the closest extreme is removed.

v f ix = v ix cos(-θ) -v iy sin(-θ) (26) v f iy = v ix sin(-θ) + v iy cos(-θ) (27)

V. RESULTS

To validate the proposed scheme, a simulation of the robot HRP-2 considering its full dynamics was performed on a ground where random objects were presented. The height of the obstacles was always smaller than the step height. As described in the paper, only the CoM and the waist are tracked by the solver. The foot position are specified using two interpolation tasks. The robot first walks in a straight line and then slightly turns left. Some snapshots of the robot are shown in Figure 5, and a video of the robot walking can be found in homepages.laas.fr/oramos/humanoids2012/walking.avi. The snapshots show the right and left foot stepping over different obstacles in different conditions. Since there is no control on the arms, they move freely to help keep the CoM position.

Figure 6 shows the tracking of the center of mass using the inverse dynamics stack of tasks. It can be observed that the CoM is properly tracked by the system. Figure 7 shows the trajectory of the right foot. The dashed line shows the trajectory that is obtained from the pattern generator, which assumes a constant ground height. As stated before, the right foot does not follow this trajectory but only uses the final position of the foot. In this experience, the two first steps of the right foot find an obstacle on the ground. It can be observed in the curves with solid line that even though the foot is not tracked, the interpolation task generates a trajectory similar to that of the pattern generator. The difference in the height z for the two first steps is due to the obstacle found by the foot, which is properly handled by the system without losing the dynamic balance. Figure 8 shows the same information for the left foot. In this case, the second, third and fourth steps found an obstacle on the ground and that is the reason why there is an observable difference in the z axis between the output of the pattern generator and the trajectory generated by the dynamic solver.

VI. CONCLUSION

The control scheme presented in this paper makes the humanoid robot able to walk on non-planar surfaces by detecting collision points and moving the foot properly to reach a larger support polygon. This approach is grounded on a mathematical reasoning and has been tested in a simulation environment where collision detection is available. The method is limited to horizontal surfaces and would fail if the ground has a large slope, in which case, a modification in the pattern generator would be needed. The implementation of the controller on the real robot requires a strategy to efficiently detect contact points. This constitutes our current work.

Fig. 2 :

 2 Fig. 2: Scheme of the usage of the pattern generator with the inverse dynamics stack of tasks.

Fig. 4 :

 4 Fig. 4: Two contact points on the sole of the robot.

Fig. 6 :

 6 Fig. 6: Evolution of the center of mass in the XY (horizontal) plane. The dashed line shows the reference of the pattern generator and the solid line shows the tracking result of the controller.

Fig. 7 :

 7 Fig. 7: Trajectory of the right foot in xyz. The dashed line shows the trajectory of the pattern generator and the solid line shows the trajectory generated by the inverse dynamics controller.

Fig. 8 :

 8 Fig. 8: Trajectory of the left foot in xyz. The dashed line shows the trajectory of the pattern generator and the solid line shows the trajectory generated by the inverse dynamics controller.

Fig. 5 :

 5 Fig. 5: Robot walking on a non-planar surface.


V Bocchetti
email: virgile.bocchetti@u-cergy.fr

H T Diep
email: diep@u-cergy.fr

Melting of Rare-Gas Crystals: Monte Carlo Simulation versus Experiments

Keywords: numbers: 34.20.Cf, 64.70.Dv, 65.20.+w, 65.40.-b:

We study the melting transition in crystals of rare gas Ar, Xe, and Kr by the use of extensive Monte Carlo simulations with the Lennard-Jones potential. The parameters of this potential have been deduced by Bernardes in 1958 from experiments of rare gas in the gaseous phase. It is amazing that the parameters of such a popular potential were not fully tested so far. In order to carry out precise tests, we have written a high-performance Monte Carlo program which allows in particular to take into account correctly the periodic boundary conditions to reduce surface effects and to reduce CPU time. Using the Bernardes parameters, we find that the melting temperature of several rare gas is from 13 to 20% higher than that obtained from experiments. We have throughout studied the case of Ar by examining both finite-size and cutoff-distance effects. In order to get a good agreement with the experimental melting temperature, we propose a modification of these parameters to describe better the melting of rare-gas crystals.

I. INTRODUCTION

Melting of crystals has always been a fascinating subject for more than a century since the discovery of the empirical Lindemann's criterion [1]. The Lindemann's criterion says that if the average of vibration amplitude u, namely u 2 , exceeds a certain value, usually 10% of the distance between nearest-neighbors, then the melting occurs. Times and over again, many authors have tried to find out microscopic precursor mechanisms that lead a crystal to melt. Until 30 years ago, one of the favorite pictures of melting is the softening of a phonon mode due to the temperature (T ). The atoms have no longer restoring forces which keep them staying close to their equilibrium positions: they move around and the system goes to a liquid state. The soft-mode picture has encountered some scepticism because in real crystals as well as in simulations one observes that well below the melting temperature (T m ), many defects, dislocations, interstitial atoms ... are excited. Therefore, it is hard to believe that the system stays in a periodic structure with propagating phonon modes up to T m . Evidence of defects is found in many works [2][3][4][5]. Another question that is unsolved in a clear manner is the form of the potential that binds the atoms together in a given lattice structure. In a microscopic point of view, the potential should come mainly from the symmetry of atomic orbitals. But ab-initio calculations are still far away from being able to use realistic hypotheses [START_REF] Pettifor | Bonding and Structure of Molecules and Solids[END_REF]. Empirical potentials have been used instead to study melting. One can mention the popular 6-12 power Lennard-Jones (LJ) potential [START_REF] Aschrof | Solid State Physics[END_REF][START_REF] Ch | Introduction to Solid State Physics, 8-th edition[END_REF], various similar power potentials, the many-body Gupta's potential [START_REF] Gupta | [END_REF], the Stillinger-Weber (SW) potential [10], and the Tersoff potential [11,12]. Two-body potentials such as the LJ one crystallize atoms in the FCC at low temperatures and nothing else; this comes from the fact that LJ potential is isotropic so the atoms are crystallized in the most dense isotropic structure, namely the FCC lattice. In order to stabilize other structures, several phenomenological potentials have been introduced, often without a microscopic justification. For example, the SW potential or the Tersoff potential stabilize the diamond structure at low temperatures. These potentials have been used with success to calculate properties of Si clusters [13] and amorphous Si crystals [14].

In this paper, we use the LJ potential to study the melting of rare-gas crystals which have the FCC lattice structure at low temperatures. It is amazing that such a simple question was not studied with precision so far in spite of an abundance of experimental data on rare gas such as Ar, Xe, Ne and Kr. Most of the melting studies concerning rare gas were done in particular cases: small clusters [14], adlayers on a substrate, etc. The main reason to avoid to study the bulk melting may be due to some technical difficulties such as periodic boundary conditions, volume expansion with temperature, etc. Previous Monte Carlo (MC) studies of bulk melting have been carried out with LJ potential but emphasize was put on the melting mechanism rather than on the precise melting temperature in real materials [2,3].

The purpose of this paper is therefore to test whether or not the experimental T m can be reproduced by MC simulation using the values of the LJ parameters deduced for rare gas in the gaseous state long time ago [15]. We will show here that by appropriate choices of technical procedures, we are able to obtain melting temperature for various rare gas directly from our simulations, unlike previous simulations [16][17][18][19] which have had recourse to various means and some thermodynamic functions to deduce it. We find in this work the melting temperatures for several rare gas higher than experimental values. A revision of the values of LJ parameters widely used in the literature for more than 50 years should be made in order to better describe the solid state of rare gas. Note that in a recent work [20], a hypothetical thermodynamic integration path is used to find the relative free energies of the solid and liquid phases, for various system sizes, at constant cutoff radius, in an attempt to explain the overestimate of the melting temperature with the LJ potential. However, due to various approximations, several results were not physically clear, in particular why the melting temperature oscillates with increasing cutoff distance.

Section II is devoted to a description of the model and our MC technique. The results are shown and discussed in Sect. III. Concluding remarks are given in Sect. IV.

II. MODEL AND METHOD

The LJ potential is given by

U = 1 2 i =j V ij (1)
with

V ij = 4ǫ a σ r ij 12 -b σ r ij 6 (2)
where ǫ , a, b and σ are constants, r ij denotes the distance between two atoms at r i and r j . Here a = 1 and b = 1. We list for convenience the values of the above constants for various rare gas and their melting temperatures, using the data from Ref. [START_REF] Ch | Introduction to Solid State Physics, 8-th edition[END_REF]. Note that the constants ǫ and σ are those deduced with some approximations using experimental data of the gaseous phase of rare gas which have been described in the pioneer paper of Bernardes [15]. The listed values of the constants should be therefore viewed as approximate values. To our knowledge, there were not papers using these constants to verify if the melting temperature is correctly obtained, either by theoretical calculations or by simulations. The absence of works dealing with this point has motivated our present work. In one of our previous works [3], we have used the LJ potential to investigate the mechanism which initiates the melting. We have counted the number of defects created with increasing temperature (T ) and found that the melting occurs when these defects interact with each other forming a kind of chains of defects which break the solid periodic state at a temperature well below the melting point. In that work, we did not use the constants of the LJ potential for any specific material. We used instead the reduced temperature and dimensionless parameters. It was not our purpose to clarify the value of melting temperature for a given kind of crystal. It is now the time to verify those constants to see whether or not they yield the correct value of T m for the melting of the solid phase of rare gas. I. Let us describe our MC technique in the study of crystal melting. In our simulations we take a system of FCC lattice of size N = 4 × L × L × L where L = 4, 5, 6, 8 and 10 namely a system of N = 256 , 500, 1024, 2048 and 4000 atoms. Unlike spin systems on lattice where the system volume does not change with T , a system of vibrating atoms expands its volume with increasing T . We use the method described in Ref. [START_REF] Allen | Computer simulation of liquids[END_REF] but with a modification that the displacement of an atom i from its current position r i is taken in a sphere of very small radius r with two arbitrary angles (θ, φ). To respect the uniform spatial distribution around r i , we use the Jacobian, namely we take an arbitrary cos(θ) between -1 and +1, instead of an arbitrary θ. φ is taken at random between 0 and 2π. Note that we first move all atoms to new positions, and change the system volume for a small arbitrary amount, then calculate the new system energy. We tune the magnitude of displacement (r) of the atoms and the volume variation in order to have an acceptation rate of displacement between 30% and 50% . This criterion is empiric but it is frequently used in MC simulation. A MC step consists of moving all atoms, each with an arbitrary displacement, and changing the system volume once. The transition probability to the new state is exp(-W kB T ) where

W = P (V new -V old ) + 1 2 (U new -U old ) + N k B T ln V old V new ( 3
)
where P is the pressure which is set to zero here (constant pressure), V old and V new are old and new system volumes, U old and U new old and new system energies. One of the most difficult tasks in melting simulations is the application of the periodic boundary conditions. This is due to random positions of atoms at the crystal boundaries. We have taken the following actions to shorten, without loosing accuracy, our simulation CPU time:

(i) each atom has a list of neighbors up to a distance r d longer than the potential cutoff distance r c . To establish the list for the first time, we have to calculate all distances and we select neighbors at r ≤ d. It takes time. The fact to choose r d > r c is to ensure that for small displacements, some neighbors with r c < r < r d can enter inside the cutoff sphere without the need to reestablish the list of neighbors (ii) when an atom moves outside a sphere of radius r 1 around its equilibrium position, its list of neighbors is updated by recalculating all pairs, neighbors up to r d fill the new list. r 1 is chosen to be equal to 20% of r 0 the nearest-neighbor distance at equilibrium (T = 0) (iii) periodic boundary conditions are applied by translating the system in all directions in a manner that avoids mismatches at the boundaries.

The last step is the most important technically. We show in Fig. 2 a snapshot of the system and its translated atoms by periodic boundary conditions. We see that if they are of the same color, one cannot distinguish the system boundaries.

In the following, we will show results with a cutoff distance r c equal to twice ℓ, the lattice constant of the FCC cell, namely r c = 2ℓ = 2r 0 √ 2 where r 0 is the nearestneighbor distance. As seen in section III D this value of r c is large enough to ensure a correct value of the melting temperature.

III. MELTING OF RARE-GAS CRYSTALS A. The case of Ar

Let us show first the results for Ar obtained by using the values listed in Table I. We will discuss next the modification necessary for obtaining the result in agreement with experiment. In order to take a correct average of physical quantities, we record spontaneous values of all physical quantities during each MC run. We have to go to several millions of MC steps before observing statistical fluctuations around equilibrium. We show an example of the energy per atom E versus MC time in Fig. 3 for N = 256 atoms at two temperatures T = 92 K and 94 K. At T = 92 K, the system is still in the solid phase. Its energy is stabilized after about one million MC steps per atom. However, at T = 94 K, the system melts after three millions of MC steps per atom: E is stabilized in the liquid phase only after such a long MC time. It is very important to emphasize that the convergence time to equilibrium depends on various MC parameters such as value of displacement magnitude r, volume variation δ etc. So, nothing can replace an observation of the timedependence of physical parameters during the simulation, although this consumes a huge computer memory. We show in Fig. 4 the energy per atom E versus T in the case of Ar with N = 256 atoms, using the parameters given in Table I. We observe here that the melting occurs at T m ≃ 93 K with a large latent heat. This value of T m is rather far from the experimental data given in Table I. We return to this point later. We show now the snapshots of the system for different temperatures in Fig. 5 below and above the melting. As seen, the system just starts to be spatially disordered at 92 K, and is well in the disordered phase (liquid) at 100 K. The radial distribution for Ar is shown in Fig. 6 at several temperatures. One sees clearly the peaks at first, second and third neighbor distances for T < T m indicating the crystalline phase, while there is no such clear distinction for T > T m where a liquid phase is set in.

We display in Fig. 7 the curve which shows the dilatation of the simulation box against the temperature.

Let us discuss now the difference between our value of T m (N = 256) = 93 K with the experimental value T m (exp)=84 K. Clearly, the parameters given in Table I which have been deduced with experimental data for gaseous Ar do not describe well the melting of solid Ar. Note however that the above value of T m is for N = 256 atoms. If we consider larger samples, T m should increase further.

In order to modify correctly the LJ parameters while respecting the known properties of Ar, we use the listed values of LJ parameters to estimate the size effect, namely to find the value of T m (N = ∞) for an infinite crystal. Once this task is done, we then we look for the parameters (ǫ, σ) which give correctly the experimental T m (exp). The effect of the cutoff distance will be examined.

B. Size effect

The size effect is an important fact when we work with numerical simulations because of the finite size of the simulation cell. In second-order phase transitions, finitesize scaling laws allow us to calculate the critical exponents by varying the system size [START_REF] Laudau | Monte Carlo Simulation in Statistical Physics[END_REF][START_REF] Hohenberg | [END_REF][24][25]. The correlation length is infinite at the critical point in the thermodynamic limit. However, in first-order transitions, the correlation length is finite at the transition point, the two phases coexist [START_REF] Barber | Phase Transitions and Critical Phenomena[END_REF][START_REF] Binder | [END_REF] and the energy is discontinuous. The correlation is finite at the transition temperature means that when the system size exceeds the correlation length the transition temperature does not depend anymore on the system size. We need not to go to the infinite limit to find saturated transition temperature. Note that if the size of the system is too small, a first-order phase transition can appear as a second-order phase transition if the correlation length exceeds the system size. As we can see in Fig. 8, fortunately when the size of the simulation cell is larger than 5 lattice constants, melting temperature T m is bounded by 102 K. It means that the correlation length at the melting temperature is smaller than 5 lattice spacings. The saturated melting temperature obtained by our MC simulations is thus 18 K (or 20%) higher than the experimental one. This disagreement come from the values of the Bernardes parameters used for the LJ potential. The reason is that these values have not been calculated by fitting with the melting temperature but they were calculated with a low-density gas using the second virial coefficient [15]. The Bernardes parameters are therefore questionable as we can see in Refs. [28][29][30]. Taken from those articles, the values of the nearest-neighbor distance, the cohesive energy and the bulk modulus calculated using the Bernardes parameters are listed in Table II. We can see the differences between experimental values and theoretical ones for all cases Ne, Ar, Kr and Xe. The Bernardes parameters also yield the high melting temperature found in our simulations with respect to the experimental one.

At this stage, it is interesting to note that Molecular Dynamics (MD) simulation of melting of a perfect crystal with periodic boundary conditions produce superheating. There is an empirical rule which states that the melting temperature of a crystal without any surface and any defects, is 20% higher than the true thermodynamic melting temperature T m . However in MC simulations, defects and dislocations are naturally created in the crystal by means of random numbers used in every MC step for atom displacements. Thus, the superheating should not exist. Agrawal et al. have shown in Ref. [31] that for Ar, with MC simulation, T m is about 15% higher than the experimental value. With our results for Ar, we find an increase of about 20%. For Kr, we find an increase of 13% for N = 256 atoms after 23 × 10 6 MC steps per atom. This increase is more important if we consider larger sizes as seen in the case of Ar. As said, the high values of T m in MC simulations are not due to the superheating as in MD simulations. Rather, we believe that these high values are due to the inaccuracy of the listed Bernardes parameters. We will propose a modification in the following

C. Modification of Lennard-Jones parameters

In order to reduce the T m value, we propose now to modify the value of ǫ, the prefactor of the LJ potential, and the coefficient σ.

Note that in papers dealing with melting in other materials by means of MD calculations or MC simulations, there have been several propositions to modify constants appearing in potentials in order to obtain a correct agreement with experimental value of T m . This is because these constants are often deduced from experimental data which are not valid for the whole temperature range. Among these papers, we can mention the case of melting of Si crystal studied by MD [32,33] and MC simulations [31,34], using the Tersoff potential [11]. Our proposition to modify some constants of the LJ potential when applied to a rare-gas crystal is certainly a necessity in order to reproduce the experimental T m .

We have done simulations with different pairs of (ǫ, σ). It turned out that T m depends essentially on ǫ. There is however an optimal value of σ which is σ = 3.44 Å corresponding to the experimental nearest-neighbor distance r 0 = 3.75 Å of Ar (cf. Table II). We show in Fig. 9 the curves obtained for two selected values of ǫ =0.008767853 and 0.008951794 which give respectively T m = 83 K and 86 K. These values of T m are in agreement with the experimental value 84 K within statistical errors. Note that the modified ǫ is about 15% smaller than the original Bernardes value ǫ = 0.01042332126.

D. Effect of cutoff distance

At this stage, a natural question we ask ourselves is "what is the effect of the cutoff distance r c ?". We know that for a long-range interaction, the longer the interaction range is the lower the energy becomes. As a consequence, the melting transition is higher. However, as r c increases, the contribution of neglected neighbors becomes smaller. From a certain value of r c , T m does not vary significantly. This is observed in Fig. 10 where T m is saturated for r c ≥ 2ℓ, i. e. r c 10.6 Å. All the results shown above for r c are valid in the discussion of the size effect and the modification of ǫ and σ.

E. The case of other rare gases

In order to show that our algorithm works well with other rare gas, we have plotted the curve of energy versus temperature, obtained for Krypton and Xenon in Fig. 11. Again here, we see that T m , even for a small size, is already higher than the experimental value for each crystal. We think that the Bernardes parameters for Kr and Xe should be modified to get an agreement with experiments as what proposed above for Ar.

IV. CONCLUSION

We have shown in this paper results of the melting temperature for rare gas, by performing extensive MC simulations with the LJ potential. We have obtained directly from our simulations physical quantities such as internal energy, lattice constant and radial distribution as functions of temperature. We have shown that melting occurs with a large latent heat and a jump in the lattice constant. Effects of system size and cutoff distance have been investigated. Let us emphasize that the theoretical LJ parameters widely used in the literature [15] yield a melting temperature higher than the experimental one as seen above, from ∼ 15% for Kr to ∼ 20% for Ar. This is not a surprise because those LJ parameters already yield theoretical NN distance, cohesive energy and bulk modulus different from corresponding experimental ones (see Table II). We have demonstrated that, in order to reduce the melting temperature to fit with experiments, it is necessary to modify the original Bernardes LJ parameter ǫ in such a way to reduce the energy at T = 0.

The effect of σ, within possible values of NN distance, is very small on T m . A good agreement on T m between experiments and simulations for Ar is obtained with the modified values given in Fig. 9. For other rare gas such as Kr and Xe, we have to proceed to a modification of their Barnardes parameters in a manner similar to that done above for Ar if we want to get a good agreement with experiments.

Figure 1 .

 1 Figure 1. (Color online) The Lennard-Jones potential for Ar with the parameters listed in TableI.

Figure 2 .

 2 Figure 2. (Color online) Upper: Snapshot of the system (magenta) with its translated atoms (gray) for periodic boundary conditions, Lower: The same snapshot of the system with its translated atoms but with the same color (online): boundaries are undistinguishable. Argon crystal at T = 95 K with N = 500 atoms.

Figure 3 .

 3 Figure 3. (Color online) Energy per atom E versus number of MC steps per atom t for an Argon crystal at T = 92 K (red) and T = 94 K (green) with N = 256 atoms.

Figure 4 .

 4 Figure 4. (Color online) Energy per atom E versus temperature T for an Argon crystal of 256 atoms, obtained with the parameters listed in TableI. One observes that the melting temperature is Tm ≃ 93 K while the experimental value is 84 K. Run of 50 millions MC steps/atom. See text for comments.

Figure 5 .

 5 Figure 5. (Color online)Instantaneous pictures of the supercell of Ar for different temperatures, respectively from left to right : 75 K (far below Tm), 92 K (close to Tm) and 100 K (above Tm).

Figure 6 .

 6 Figure 6. (Color online) Radial distribution of Ar for different temperatures: below melting T = 74 K (discontinued red) and 92 K (discontinued green), above melting T = 94 K (solid blue) and 100 K (solid magenta), with N = 256 atoms.

Figure 7 .

 7 Figure 7. (Color online)Evolution of the lattice constant against temperature T . As we can see, there is a jump of this quantity which occurs at the same temperature of the jump of the energy in Fig. 4. Argon crystal, N = 256.

Figure 8 .

 8 Figure 8. (Color online)Energy per atom versus T with various system sizes from N = 256 (4 FCC lattice cells) to N = 4000 (10 FCC lattice cells). The arrows indicate Tm for the smallest and largest sizes. See text for comments.

Figure 9 .

 9 Figure 9. (Color online) Energy per atom versus T for two selected values of ǫ, with σ = 3.44 Å, for an Ar crystal with N = 500. The arrows indicate Tm for the two indicated values of ǫ. See text for comments.

Figure 10 .

 10 Figure 10. (Color online) Energy per atom versus T with various values of rc for Ar crystal with N = 500. The left arrow indicates Tm for rc = 1.5ℓ and the right arrow indicates Tm for rc = 2ℓ ∼ 10.6 Åand 3ℓ ∼ 15.9 Å. Note that ℓ is the FCC lattice constant which is equal to r0 √ 2 where r0 = 3.75 Å is the NN distance. See text for comments.

Figure 11 .

 11 Figure 11. (Color online) Upper: Energy versus temperature for a Krypton perfect crystal with N = 256 atoms. One observes that Tm = 132 K while the experimental value is 117 K. Lower: Energy versus temperature for a Xenon perfect crystal with N = 500 atoms. One observes that Tm = 191 K while the experimental value is 161 K. These curves have been obtained with the Bernardes values of parameters.

Table I .

 I Lenard-Jones parameters a

	Element	ǫ (eV)	σ( Å)	Tm (K)
	Ar	0.01042332126	3.40	84
	Kr	0.01404339691	3.65	117
	Ne	0.00312075487	4.74	24
	Xe	0.01997283116	3.98	161
	He	0.00087381136	2.56	-

a Reference

[START_REF] Ch | Introduction to Solid State Physics, 8-th edition[END_REF]

Table II .

 II Experimental (Exp.) and theoretical (Th.) values of NN distance r0, cohesive energy u0 and bulk modulus B0 a for comparison. See text for comments.

			Ne	Ar	Kr	Xe
	r0	Exp.	3.13	3.75	3.99	4.33
		Th.	2.99	3.71	3.98	4.34
	u0	Exp.	-0.02	-0.08	-0.11	-0.17
		Th.	-0.027	-0.0889	-0.12	-0.172
	B0	Exp.	1.1	2.7	3.5	3.6
		Th.	1.81	3.18	3.46	3.81

a Reference

[28][29][30]