Stéphane Durand
email: stephane.durand@inria.fr

Bruno Gaujal
email: bruno.gaujal@inria.fr

Complexity and Optimality of the Best Response Algorithm in Random Potential Games

In this paper we compute the worst-case and average execution time of the Best Response Algorithm (BRA) to compute a pure Nash equilibrium in finite potential games. Our approach is based on a Markov chain model of BRA and a coupling technique that transform the average execution time of this discrete algorithm into the solution of an ordinary differential equation. In a potential game with N players and A strategies per player, we show that the worst case complexity of BRA (number of moves) is exactly N A N -1 , while its average complexity over random potential games is equal to e γ N + O(N ), where γ is the Euler constant. We also show that the effective number of states visited by BRA is equal to log N + c + O(1/N ) (with c e γ ), on average. Finally, we show that BRA computes a pure Nash Equilibrium faster (in the strong stochastic order sense) than any local search algorithm over random potential games.

Introduction

The question of computing Nash Equilibria (NE) in games is a central question in algorithmic game theory and has been investigated of many papers. The most classical result is in [START_REF] Daskalakis | The complexity of computing a nash equilibrium[END_REF], showing that the problem of computing NE in finite games is PPAD complete.

Potential games have been introduced in [START_REF] Rosenthal | A class of games possessing pure-strategy nash equilibria[END_REF] and have proven very useful, especially in the context of routing problems in networks, first mentioned in [START_REF] Beckman | Studies in the Economics of Transportation[END_REF] and exhaustively studied ever since, in the transportation as well as computer science literature, see for example [START_REF] Orda | Competitive routing in multeuser communication networks[END_REF][START_REF] Gallager | A minimum delay routing algorithm using distributed computation[END_REF][START_REF] Wardrop | Some theoretical aspects of road traffic research. part ii[END_REF]. They have also been heavily investigated in the context of distributed optimization (see for example [START_REF] Roughgarden | Selfish Routing and the Price of Anarchy[END_REF]). In [START_REF] Fabrikant | The complexity of pure nash equilibria[END_REF][START_REF] Ackermann | On the impact of combinatorial structure on congestion games[END_REF] the authors show that the computation of NE for general potential games is PLS complete (Polynomial Local Search complete). As for PPAD, this complexity class is believed to be different from P .

The best response dynamics is one of the most basic tool in game theory. The original proof of the existence of a Nash Equilibrium by Nash [START_REF] Nash | Equilibrium poits in n-person games[END_REF] can be seen as the proof of existence of a fixed point of the best response correspondence (best response is called countering in [START_REF] Nash | Equilibrium poits in n-person games[END_REF]). It has been well-known for a long time that the Best Response Algorithm converges in finite time to a pure NE in potential games [START_REF] Monderer | Potential games[END_REF]. So BRA is a natural candidate for computing Nash equilibria.

In this paper, we analyze the performance of BRA over potential games with N players, each with A possible strategies. It is well known that the convergence time of BRA over potential games can be exponential in the number of players (see for example [START_REF] Fabrikant | The complexity of pure nash equilibria[END_REF]). Here, we confirm this by showing that the worst case complexity of BRA (number of plays) is exactly N A N -1 . Special cases, such as graphical potential games have been analyzed in [START_REF] Babichenko | Graphical potential games[END_REF] by showing an equivalence between the potential of such games and Markov fields. In other special cases such as scheduling congestion games with identical tasks, it has been show that BRA takes at most N steps before finding a NE [START_REF] Even-Dar | Convergence time to nash equilibria[END_REF]. Extensions with positive and negative externalities also have a linear complexity [START_REF] Feldman | Convergence of best-response dynamics in games with conflicting congestion effects[END_REF].

However the average complexity of BRA over all potential games has attracted surprisingly little attention. Random (non potential) games with two players have been studied in [START_REF] Imre | Nash equilibria in random games[END_REF]: With two IID utility matrices of size A × A, the computation of a NE is O(A 3 log log A) with high probability using a rather sophisticated algorithm.

Our main contribution is to show that for potential games with N players, E[M BRA ], the average number of strategy profiles visited by BRA before convergence, is

E[M BRA ] = log(N ) + C + O(1/N ) (
where C e γ , γ being the Euler constant). We also show that the average number of comparisons performed by the algorithm is equal to e γ (A -1)(N -1) + o(AN ). This could be intuitively explained by the fact that random potential games have a lot of pure NE [START_REF] Rinott | On the number of pure strategy nash equilibria in random games[END_REF]. In our framework, potentials are IID random variables so that, on average, one action profile out of (A -1)N + 1 profiles is a NE while in the worst case, a potential game may have a single pure NE. This is only a partial explanation, however. This does not explain the fact that the complexity does not depend on the number of actions, nor the value of the constant factor, e γ ≈ 1.78.

We further show that the Best Response Algorithm computes a pure Nash Equilibrium faster than any algorithm based on player's local information, not only in average but also in the strong stochastic order sense.

Missing proofs and additional details (numerical simulations, analysis of alternative algorithms) are given in a long version of this paper, available in HAL Archive ( [START_REF] Durand | Complexity and Optimality of the Best Response Algorithm in Random Potential Games[END_REF]).

Coupling and Markovian Analysis

The main idea of our approach is to see the evolution of BRA in a random environment as a dynamical system, whose behavior can be computed using differential equations. This will allow us to compute the exact asymptotics of the average complexity in N and A, not only O(.) bounds. Second moments of T BRA and of M BRA can also be computed by the same approach (see [START_REF] Durand | Complexity and Optimality of the Best Response Algorithm in Random Potential Games[END_REF]).

The first step ( § 4.2) is to construct an approximation of the behavior of BRA over a potential game. This approximation is called IFA in the paper, for Intersection-Free Approximation because it discards strategies already explored by BRA. We show that the execution time of BRA is smaller than the execution time of its IFA approximation for the strong stochastic order. This is done by constructing a non-trivial coupling between both executions. This powerful technique is exploited to our great benefit here.

The second and most important step ( § 4.4) is to consider one run of the IFA approximation of BRA as a trajectory of a Markov chain over the continuous space of potentials. Doing so, the average complexity is transformed into the average hitting time of an absorbing state of the Markov chain. The theory of Markov chains implies that this average hitting time satisfies a Poisson differential equation. Thus, the average complexity of BRA is given by the solution of a system of ordinary differential equations. This system happens to have a solution in closed form whose asymptotics in N and A can be computed by taking integrals over initial states.

As for the proof of optimality of BRA among all local search algorithms ( § 5), our approach is based once again on a coupling argument. While using coupling techniques is more classical in this context (comparison of algorithms), this particular case retains some originality because the coupling used here is not built off-line but is being constructed on the fly while the algorithm runs.

Best Response Algorithm and Potential Games

We consider a game with a finite number N of players, each with A strategies.

Definition 1 (N -player game). A game is a tuple

G def = G(N , A, u) with
a finite set of players N = {1, . . . , N }; a finite set A k of pure strategies for each player k ∈ N .

The set of strategy profiles or states of the game is

A def = A 1 × A 2 × • • • A N .
-The players' payoff functions u k : A → R.

We define the best response correspondence br k (x) as the set of all strategies that maximizes the payoff for player k under profile x = (x 1 , . . . , x N ): br k (x)

def = argmax α∈A k u k (α; x -k ) .
A Nash equilibrium (NE) is a fixed point of this correspondence, i.e. a profile x * such that x * k ∈ br k (x * ) for every player k.

Definition 2 (Potential games and its generalizations).

A game is an (exact) potential game [START_REF] Monderer | Potential games[END_REF] if it admits a function (called the potential) Φ : A → R such that for any player k and any unilateral deviation of k from strategy

profile x to x : u k (α, x -k ) -u k (α , x -k ) = Φ(α, x -k ) -Φ(α , x -k ).
A game is a generalized ordinal potential game [START_REF] Monderer | Potential games[END_REF] (or G-potential game for short) if there exists a potential function Φ : A → R such that, for any player k

and any state x, u k (α, x -k ) > u k (α , x -k ) ⇒ Φ(α, x -k ) > Φ(α , x -k ).
A game is a best-response potential game [START_REF] Voorneveld | Best-response potential games[END_REF] (or BR-potential game for short) if there is Φ : A → R such that for any player k and strategy profile x,

br k (x) = argmax α∈A k Φ(α, x -k ) .
As shown in [START_REF] Voorneveld | Best-response potential games[END_REF], exact potential games are BR-potential games, but there exist G-potential games that are not BR-potential games. In the following, we will consider the most general case (i.e. all games that are either BR-potential or G-potential games) and call them potential games for simplicity.

We consider a general version of the Best Response Algorithm (BRA) with uniform choice over all possible best responses when ties occur and where the next player is selected according to a revision function R(.), that may depend of the whole past of the algorithm. We assume that this function is weakly fair: each player appears infinitely often in the sequence of plays induced by R, almost surely. This revision function can be deterministic (for example, round-robin: R(t) = t mod N ) or random (for example, Bernoulli where the next player is chosen according to an probability distribution ρ (the revision law): ∀k ∈ N , P(R(t) = k) = ρ k ). In that case, weak fairness implies that the probability of choosing any player k is strictly positive (∀k ∈ N , ρ k > 0). It is well known (see [START_REF] Monderer | Potential games[END_REF]) that for any potential game G, Algorithm 1 converges in finite time, almost surely, to a Nash Equilibrium of G.

Worst Case Complexity

In this section, we analyze the time complexity of BRA. More precisely, we consider three measures (related to each other). The first one is T BRA , the number of iterations (or the number of times that the function br is called) before BRA reaches a Nash equilibrium. A related measure is the total number of comparisons used by BRA (denoted C BRA ). One should expect that C BRA ≈ (A -1)T BRA . Finally, another interesting quantity is the number of different states visited by BRA (denoted M BRA ). This is called the number of moves done by BRA before convergence to a Nash equilibrium (NE). Of course, M BRA T BRA .

These quantities depend on the game over which BRA is run, on the initial state x 0 and on the infinite sequence of revision players R. It should also be clear that they are functions of the game only through the potential Φ, so we denote by T BRA (Φ, x 0 , R) the number of steps before convergence of Algorithm BRA for a game with potential Φ, starting in state x 0 , under the condition that the sequence of players is given by R.

In the worst case, for some weakly fair revision functions R, T BRA (Φ, x 0 , R) can be unbounded because the revision sequence induced by R can be arbitrarily bad: one player might appear too few times to guarantee convergence in any bounded time. When R is the round-robin function, the time for convergence is finite but can still be very large, as shown in the following theorem.

Theorem 1. In the worst case,

T BRA (Φ, x 0 , round-robin) = N A N -1 .
It is well known that the worst case complexity of BRA is exponential in the number of players (see for example [START_REF] Fabrikant | The complexity of pure nash equilibria[END_REF]). The version of this result given here (for round robin revision and generalized potential games) is only given for the record (the proof is given in [START_REF] Durand | Complexity and Optimality of the Best Response Algorithm in Random Potential Games[END_REF]).

Average Complexity of BRA

Randomization

In the following we will randomize over the potential games over which BRA is used. Since the behavior of BRA only depends on the potential function, we randomize directly over the potential Φ.

We consider a randomization over all games, uniformly over all possible orders for the potentials. On one hand, this is the classical average complexity approach when no additional information is known about the games (the same approach is used in [START_REF] Imre | Nash equilibria in random games[END_REF] for 2 player games for example). This yields IID potential for all profiles, as explained below. On the other hand, some may argue that uniformly random games are not generic in some sense and a good performance of BRA on average does not necessarily translate in good performances for "real word" games. In any case, this is a first step that must be taken in absence of additional information about specific games that one may want to study.

There are several equivalent ways to do this randomization. The first one is based on the fact that the complexity of the algorithm does not depend on the actual values of the potential of the states but only on the comparisons between them. When two potentials are equal, a strict order between them is chosen uniformly. Therefore, the natural randomization is to consider the linear extensions (total orders) of all possible partial orders over the set A and pick one uniformly. The number of total orders on A is the number of permutations on A, namely (A N )!.

The second (equivalent) randomization is the following: The potentials of all states x are chosen independent, identically distributed according to an arbitrary distribution F admitting a density w.r.t. the Lebesgue measure.

Both randomizations are equivalent. Indeed, take any k states x 1 , . . . , x k in A. In both cases, P(Φ(x

1 ) > Φ(x 2 ) > • • • > Φ(x k )) = 1/k!. Now, since F is increasing, F -1
is well-defined and we get P(Φ(x) > Φ(x )) = P(F -1 (Φ(x)) > F -1 (Φ(x ))). Note that F -1 (Φ(x)) is uniformly distributed on [0, 1]. Therefore, with no loss of generality, one can assume that the potential of all the states are i.i.d., uniformly distributed on [0, 1]. This randomization is used in the following.

Intersection-Free Approximation

The direct analysis of the behavior of BRA over a random potential is difficult because, over time, more and more states have been visited by the algorithm. Thus, its behavior is non-homogeneous in time. To avoid this difficulty, we consider a new model, called the Intersection-Free Approximation (IFA) in the following. Under the Intersection-Free Approximation, every time a new player (say k) has to compute its best response in a state (say x), it compares Φ(x) with the potential of its A -1 other possible strategies, as for the real BRA. Here however, we assume that those A -1 states have not yet been visited during the previous steps of the algorithm. Note that under the real behavior of BRA, it could happen that some of these states have already been compared at a previous step of the algorithm, by another player (this will be called an intersection in the following). Under the Intersection-Free Approximation, the states visited by the algorithm are always "new" states, never compared before with any other states.

More formally, the algorithm BRA under IFA can be written as follows.

= t + 1; if k / ∈ L then Generate IID potentials Φ(α, x -k ), α ∈ A k \ {x k } unif. on [0, 1]; Compute best response: α k := argmax β∈A k Φ(β; x -k ); if α k = x k then L := L ∪ {k} else L := {k}; x k := α k ; until size(L) = N ;
Let us recall that C BRA (resp. T BRA , M BRA ) is the number of comparisons (resp. number of steps, number of moves) taken by BRA before convergence and let us define C IF A (resp. T IF A , M IF A ) to be the number of comparisons (steps, moves) of BRA under the intersection-free approximation. By definition, the worst case complexity of IFA under a round-robin revision sequence is infinite. However, its average complexity is the same as for BRA, as shown by the following lemma.

Lemma 1 (BRA and IFA are asymptotically equivalent). Under the foregoing notations and using a round-robin revision function, the following comparisons hold, where st is the strong stochastic order:

1. C BRA st C IF A (equivalently, [19] ∀t ∈ R, P(C BRA > t) P(C IF A > t)).
2. If I is the total number of intersections in BRA, then

T BRA st T IF A + I A-1 . 3. E[T BRA ] = E[T IF A ] + o(1) and E[C BRA ] = E[C IF A ] + o(1), 4. E[M BRA ] = E[M IF A ] + o(1).
The proof of the lemma is available in [START_REF] Durand | Complexity and Optimality of the Best Response Algorithm in Random Potential Games[END_REF]. It is based on the construction of a coupling between the executions of BRA with and without IFA. The assumption that the revision function is round-robin for BRA and for IFA does not play a big role in the proof, and it could be removed. However, the following section, asserting the optimality of round-robin implies that extending the proof to more general revision functions has a limited interest.

Round-Robin and Other Revision Sequences

As for the worst case analysis, the revision sequence influences the average time complexity of the algorithm. We show that on average round-robin is asymptotically the best one.

Lemma 2 (Asymptotic optimality of round-robin). For any revision func

- tion R, E T BRA (Φ, x 0 , round-robin) E T BRA (Φ, x 0 , R) + (N )
, where the expectation is taken over all potentials Φ and all initial states x 0 and (N ) goes to zero when N goes to infinity.

The proof is again available in [START_REF] Durand | Complexity and Optimality of the Best Response Algorithm in Random Potential Games[END_REF]. It uses the comparison with IFA. In the rest, we focus on round-robin revision functions and omit it in the notations, unless specified otherwise.

Complexity Analysis

We will be analyzing the intersection-free approximation of the behavior of BRA, under a round-robin revision sequence, with no further reference to this.

Let us consider the intersection-free approximation and let y be the potential of the current state x: (y def = Φ(x)). Let k be the number of players that have already played best response without changing the profile. This number of "satisfied" players can replace the explicit set L used in Algorithm 2 when the revision sequence is round-robin. The evolution at the next step of BRA under IFA is as follows. The kth player computes its best response. The player has A -1 new strategies whose potential must be compared with the current potential (y). As mentioned before, we can assume that the potentials of those

a def = A -1 strategies are IID, uniformly distributed in [0, 1].
With probability y a none of the new strategies beat the current choice. The state remains at y, one more player is satisfied and it is the turn of the k + 1-st player to try its best response.

With probability 1 -y a , one of the new strategies is the best response. The current state moves to a new state where the number of satisfied players is set back to 1 and the potential increases to a value larger than u > y with probability 1 -u a .

Let Y t be the potential at step t (Y t ∈ [0, 1]) and K t be the current number of consecutive players whose best response did not change the current potential (K t ∈ {1, 2, . . . , N }) (number of satisfied players). The previous discussion says that the couple (Y t , K t ) is a discrete-time, continuous-space Markov chain whose kernel is:

P (Y t+1 , K t+1 ) = (y, k + 1) (Y t , K t ) = (y, k) = y a ,
and, for any u > y,

P (Y t+1 , K t+1 ) ∈ ([u, 1], 1) (Y t , K t ) = (y, k) = 1 -u a .
All the other transitions have a null probability.

Let m(y, k) be the number of moves of IFA before convergence when the current state of the Markov chain is equal to (y, k).

With probability y a , the next player does not change its choice so that m(y, k) = m(y, k + 1). With probability density au a-1 the next player finds a new best response with potential u so that one move is taken and m(y, k By definition, the boundary conditions are: ∀y, M (y, N ) = 0 (the current state is NE when all players agree on this) and ∀k, M (1, k) = 0 (the potentials are all bounded by 1, so a state with potential 1 is guaranteed to be a NE).

) = 1 + m(u, 1). Let M (y, k) = E[m(y, k)].
By setting B(y) def = 1 y au a-1 (M (u, 1) + 1)du, we get the following system of integral equations

               M (y, 1) = y a M (y, 2) + B(y), M (y, 2) = y a M (y, 3) + B(y), . . . = . . . M (y, N -2) = y a M (y, N -1) + B(y), M (y, N -1) = B(y). (1)
Successive substitution of M (y, 2), . . . , M (y, N -1) in the first equality yields M (y, 1) = B(y)H(y) where H(y) -2) . Differentiating w.r.t. y, one gets an ordinary differential equation in M (y, 1):

def = 1 + y a + • • • + y a(N
dM (y, 1) dy + (ay a-1 H - 1 H dH dy )M (y, 1) = -ay a-1 H.
The equation is of the form ḟ + gf = h. Using the boundary condition M (1, 1) = 0, its generic solution is

M (y, 1) = e -Q(y) 1 y au a-1 H(u)e Q(u) du. (2)
where

Q(y) def = y 0 au a-1 H(u) - 1 H(u) dH(u) du du = -log(H(y))+ y 0 au a-1 H(u)du.
(3) The average number of profile changes in the execution of the algorithm starting from an arbitrary profile is E[M IF A ] = 1 0 M (y, 1)dy. Since M (y, 1) is decreasing in y, E[M IF A ] is upper-bounded by M (0, 1). Using Q(0) = 0, H(0) = 1 and replacing Q and H by their values,

M (0, 1) = 1 0 exp N -2 i=0 u a(i+1) i + 1 au a-1 du = 1 0 exp N -2 i=0 v i+1 i + 1 dv (with v = u a ) = 1-1 N 0 exp N -1 i=1 v i i dv + 1 1-1 N exp N -1 i=1 v i i dv 1-1 N 0 exp ∞ i=1 v i i dv + 1 N exp N -1 i=1 1 i (4)
Furthermore, this bound is tight, up to an additive constant (see [START_REF] Durand | Complexity and Optimality of the Best Response Algorithm in Random Potential Games[END_REF]).

Let us now consider the average number of comparisons made by BRA under the intersection-free assumption. Let C(y, k) be the average number of comparisons starting in a state with potential y and k players have played without changing their strategy. The Poisson equation for C(y, k) is : The solution of this differential system can be obtained in closed form, using a similar approach as for M (y, 1).

C(y, k) = y a (C(y, k + 1) + a) +
C(y, 1) = a N -2 i=0 y ai exp - N -1 i=1 y ai -1 i .

The average number of comparisons is E[C

IF A ] = 1 0 C(y, 1)dy.
For all y < 1,

C(y, 1) = a ∞ i=0 y ai exp N -1 i=1 1/i exp - ∞ i=1 y ia /i + o(aN ) (6) = a 1 1 -y a (N -1)e γ (1 -y a ) + o(aN ) + O(1) (7)
= a(N -1)e γ + o(aN ),

where γ is the Euler constant (γ ≈ 0.5772...). Therefore, the same equality holds for the integral, equal to

E[C IF A ].
The results of this section, together with Lemma 1, lead to the following theorem, the main result of the section. The average complexity E[T BRA ] can be split into two parts: The number of plays before reaching a NE and the number of plays needed to check if a state is indeed a NE. This last part takes exactly N -1 steps in Algorithm 1: The players have to play one by one to fill up set L. This means that a NE equilibrium is reached on average as soon as e γ -1 ≈ 78% of the players have played once. The second moments of the number of steps and the number of moves under IFA can be computed similarly (see [START_REF] Durand | Complexity and Optimality of the Best Response Algorithm in Random Potential Games[END_REF]). In both cases, the standard deviations are of the same order as the means.

Optimality of BRA

In this section, we prove that BRA finds a Nash equilibrium faster than any local search algorithm (defined in Section 5), in the strong stochastic order sense.

By definition a Local Search Algorithm can only access the payoff matrix, one player at a time. This access is often called a query in the literature. Once the payoff of a strategy profile has been obtained, it is stored in memory and can re-used later by the algorithm without an additional query.

In addition to queries, a local search algorithm can use any arithmetic operation, draw random variables and choose a strategy for all players.

Any local search algorithm can be written in the following form, based on the history of the execution, H t , that corresponds to the amount of information gathered by the algorithm up to step t. The functions J and R used in the inner loop are arbitrary functions that choose the next state as well as the next player to play, according to the whole history of the process. These functions can be deterministic or random. Testing if x(t + 1) is NE is not detailed. Notice, however, that it can only be done when all the payoff vectors for all the players in state x(t + 1) have been stored in memory.

The complexity of a local search algorithm A is defined as the total number of its payoff vector queries (denoted T A ).

Theorem 3 (Optimality of BRA). Let A be any local search algorithm that computes a Nash Equilibria in potential games. Under the foregoing randomization, and choosing the starting point x 0 uniformly among all states, ∀t 0, P(T BRA t|R = R A ) P(T A t), where R A is the revision sequence constructed in A.

The proof is reported in [START_REF] Durand | Complexity and Optimality of the Best Response Algorithm in Random Potential Games[END_REF]. Combining this theorem with Lemma 2 establishes the optimality of BRA with round-robin.

Conclusion and Perspectives

The best response algorithm is one on the most basic object in game theory. In this paper, we prove it has a linear complexity on average over uniformly randomized potential games. Furthermore, BRA is optimal in the class of local search algorithms when one has no information about the structure of the potential game. Does all this make BRA the perfect algorithm to compute NE in general? We believe that the answer is no because BRA suffers from several drawbacks. First, it does not tolerate simultaneous plays. Second, it requires to know the entire payoff vector of a player before choosing its strategy. Other drawbacks include high sensitivity on the order of play and on noisy perturbations on the payoffs. Designing algorithms that do not suffer from these drawbacks is the object of our future investigations.

Algorithm 1 :

 1 Best Response Algorithm (BRA) Input: Game utilities (u k (•)); Initial state (x := x(0)); Weakly fair revision function R; List of satisfied customers, initially empty: L := ∅; repeat Pick next player k := R(t); t := t + 1; if x k ∈ br k (x) then Update strategy for player k to any x k ∈ br k (x); L := ∅; L := L ∪ {k}; until size(L) = N ;

Algorithm 2 :

 2 BRA algorithm under IFA Input: Initial state (x(0)); Revision function R; Set of satisfied players, initially empty L := ∅. repeat Pick next player k := R(t); t :

1 y

 1 The previous one step analysis of m(y, k) makes M (y, k) satisfy a forward Poisson equation: M (y, k) = y a M (y, k + 1) + au a-1 (M (u, 1) + 1)du.

= 1 - 1 N 0 dv 1

 1101 -v + e γ + O(1/N ) = log(N ) + e γ + O(1/N ).

1 y

 1 au a-1 (C(u, 1) + a)du, with the boundary conditions C(1, 1) = a(N -1) and C(y, N ) = 0.

Theorem 2 (

 2 Average complexity of BRA). Under the round-robin revision sequence, the average complexity of BRA over a potential game satisfies: (i) Average number of moves: E[M BRA ] = log(N ) + c + O(1/N )., where c e γ (ii) Average number of comparisons : E[C BRA ] = e γ AN + o(AN ). (iii) Average number of steps: E[T BRA ] = e γ N + o(N ).

  Algorithm 3: A general local search algorithm Initial storage reduced to the initial profile: H 0 := {(x(0)}. repeat Select next player: k := R(H t ); Query payoff vector of k under current state: u k (•, x -k (t)); Store the new visited states and their payoffs in memory: H t+1 := H t ∪ {((α, x -k (t)), u k (α, x -k (t))) α∈A k }; jump to next state x(t + 1) := J(H t+1 ); Set stop := 1 if the current state is a NE; t := t + 1; until stop;

Acknowledgement

This work was partially supported by LabEx Persyval-Lab


Stéphane Lanteri

Claire Scheid

Jonathan Viquerat

High order DGTD

teaching and research institutions in France or abroad, or from public or private research centers.

Nanophotonics is the field of science and technology which aimed at establishing and using the peculiar properties of light and light-matter interaction in various nanostructures. Nanophotonics includes all the phenomena that are used in optical sciences for the development of optical devices. Because of its numerous scientific and technological applications (e.g. in relation to telecommunication, energy production and biomedicine), nanophotonics represents an active field of research increasingly relying on numerical modeling beside experimental studies. Plasmonics is a related field to nanophotonics involving metalic nanostructures whose optical scattering is dominated by the response of the conduction electrons [START_REF] Maier | Plasmonics-Fundamentals and applications[END_REF]. The numerical modeling of light interaction with metallic nanostructures requires to solve the system of time-domain Maxwell equations coupled to appropriate models of physical dispersion in the metal such as the Drude and Drude-Lorentz models. In this talk, we will report on our recent efforts aiming at the development of a family of high order finite element type time-domain solvers for the numerical treatment of nanophotonics applications in the linear regime. The proposed method is an extension of the so-called DGTD (Discontinuous Galerkin Time-Domain) method that was initially proposed for the simulation of electromagnetic wave propagation in non-dispersive heterogeneous media at microwave frequencies [START_REF] Hesthaven | Nodal High-Order Methods on Unstructured Grids. I. Time-Domain Solution of Maxwell's Equations[END_REF]- [START_REF] Fezoui | Convergence and Stability of a Discontinuous Galerkin Time-Domain Method for the 3D Heterogeneous Maxwell Equations on Unstructured Meshes[END_REF]. For the numerical treatment of dispersion models in metals, we have adopted an Auxiliary Differential Equation (ADE) technique leading to solve the time-domain Maxwell equations coupled to a ODEs.

  HIGH ORDER DGTD SOLVER FOR THE NUMERICAL MODELING OF NANOSCALE LIGHT/MATTER INTERACTION ST ÉPHANE LANTERI * , CLAIRE SCHEID † AND JONATHAN VIQUERAT * * Inria Sophia Antipolis-Méditerranée Research Center, France. email: {Stephane.Lanteri,Jonathan.Viquerat}@inria.fr † University of Nice-Sophia Antipolis, France. email: Claire.Scheid@unice.fr


Stéphane Éphane Lanteri
email: stephane.lanteri@inria.fr

Raphaël Léger
email: raphael.leger@inria.fr

Diego Paredes
email: diego.paredes@pucv.cl

Claire Scheid
email: claire.scheid@unice.fr

Frédéric Valentin
email: valentin@lncc.br

Rapha Ël

L Éger

F R Éd

Éric Valentin

A multiscale hybrid-mixed method for the Maxwell equations in time-domain

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

This work proposes a Multiscale Hybrid-Mixed (MHM) method for the Maxwell equation in time domain. The MHM method is a consequence of a hybridization procedure, and emerges as a method that naturally incorporates multiple scales while provides solutions with high-order precision. The computation of local problems is embedded in the upscaling procedure, which are completely independent and thus may be naturally obtained using parallel computation facilities. In this talk, we present the new MHM method for the two-dimensional Maxwell equations in time domain (Transverse Magnetic mode). We address some theoretical aspects of the method and propose an extensive numerical validation. We conclude that the MHM method is naturally shaped to be used in parallel computing environments and appears to be a highly competitive option to handle realistic multiscale hyperbolic boundary value problems with precision on coarse meshes.


Daniele De Seta

Yann Nguyen

Damian Bonnard

Evelyne Ferrary

Benoit Godey

David Bakhos

Michel Mondain

Olivier Deguine

Olivier Sterkers

Daniele Bernardeschi

Isabelle Mosnier
email: isabelle.mosnier@aphp.fr

The Role of Electrode Placement in Bilateral Simultaneously Cochlear Implanted Adult Patients

Keywords: bilateral implantation, speech perception, electrode position, cochlear implant, angular depth of insertion, cochlear size

Objective: To evaluate the influence on hearing performance of the electrode placement in adult patients simultaneously and bilaterally cochlear implanted.

Study Design: Case series with planned data collection

Setting: Tertiary referral university centers

Subjects and Methods:

The postoperative CT scan of nineteen patients simultaneously and bilaterally implanted with a long straight electrode array was studied. The size of the cochlea was measured considering the major cochlear diameter and the cochlear height. The electrode-to-modiolus distance for the electrodes positioned at 180-and 360-degrees, and the angular depth of insertion of the array were also measured. Speech perception was assessed at 1-year and at 5-years postimplantation using disyllabic words lists in quiet and in noise, with the speech coming from the front, and a cocktail-party background noise coming from 5 loudspeakers.

Results: At 1-year postimplantation, the electrode-to-modiolus distance at 180-degrees was correlated with the speech perception scores in both quiet and noise. In patients with a full electrode insertion, no correlation was found between the angular depth of insertion and hearing performance. The speech perception scores in noise gradually declined as a function of the number of inserted and active electrodes. No relationship between electrode position and speech scores was found at 5-years postimplantation.

Conclusion:

In adult patients simultaneously and bilaterally implanted, the use of a long straight array, the full electrode array insertion, and the proximity to the modiolus might be determining factors to obtain the best speech performance at 1-year, without influence on the speech scores after long-term use.

INTRODUCTION

The preservation of the inner ear structures during the insertion of cochlear implant, together with the identification of the ideal site of stimulation in the cochlea, should allow the best hearing performance. As a consequence, the quality of insertion of the cochlear implants has been extensively studied during the last decades [START_REF] Shepherd | Electrical stimulation of the auditory nerve: the effect of electrode position on neural excitation[END_REF][START_REF] Skinner | Determination of the position of nucleus cochlear implant electrodes in the inner ear[END_REF][START_REF] Skinner | CT-derived estimation of cochlear morphology and electrode array position in relation to word recognition in Nucleus-22 recipients[END_REF][START_REF] Aschendorff | Quality control after insertion of the nucleus contour and contour advance electrode in adults[END_REF][START_REF] Holden | Factors affecting open-set word recognition in adults with cochlear implants[END_REF] . In this context, three parameters have been more accurately investigated: the translocation of the array with the subsequent basilar membrane rupture, the depth of insertion of the electrode array, and the proximity of the electrodes to the spiral ganglion cells. To date, it is not clear how the position of the electrode in the cochlea can impact the hearing performance results, since many variables may influence this outcome. All the currently available electrode arrays have their own specific length, diameter, shape, and physical properties that influence the trajectory during the insertion and determine the final position in the cochlear lumen. Furthermore, variations in human cochlear anatomy, as well as the intersubject variability, have been described in several studies [START_REF] Erixon | Variational anatomy of the human cochlea: implications for cochlear implantation[END_REF][START_REF] Rask-Andersen | Human cochlea: anatomical characteristics and their relevance for cochlear implantation[END_REF][START_REF] Martinez-Monedero | Cochlear coiling pattern and orientation differences in cochlear implant candidates[END_REF] , whereas little is known about the intrasubject difference, i.e. the differences between the two ears.

Considering the hearing performance after cochlear implantation, intraindividuals variability in speech perception scores has been demonstrated among bilaterally cochlear implanted recipients [START_REF] Litovsky | Simultaneous bilateral cochlear implantation in adults: a multicenter clinical study[END_REF][START_REF] Mosnier | Speech performance and sound localization in a complex noisy environment in bilaterally implanted adult patients[END_REF] . In fact, in a prospective multicenter study, poor performance of one or both ears was reported at 1-year postimplantation in about 40% of simultaneously implanted patients with similar hearing loss history between the two ears (hearing deprivation, duration of deafness, etiology) [START_REF] Mosnier | Speech performance and sound localization in a complex noisy environment in bilaterally implanted adult patients[END_REF] . An explanation for poor hearing performance and/or asymmetry between the two ears could be differences in the electrode position within the cochlea [START_REF] Holden | Factors affecting open-set word recognition in adults with cochlear implants[END_REF] . The aim of the present study is to explore the correlation between speech performance and electrode placement parameters in patients simultaneous and bilaterally implanted, and to investigate whether cochlear anatomy differences could explain inter-and intraindividual differences in hearing performance.

MATERIALS AND METHODS

Selection criteria and subjects

Study participants were 19 adult patients presenting a post-lingual bilateral profound or total hearing loss. Specific subject demographics are summarized in Table 1. The duration of deafness, of hearing deprivation, of hearing aid use, and the etiologies were similar for both ears. Enrolling criteria, speech perception evaluation setting, and results at 1-and at 5-years have been previously reported [START_REF] Mosnier | Speech performance and sound localization in a complex noisy environment in bilaterally implanted adult patients[END_REF][START_REF] Seta | Does the Anatomy of the Cochlea influences the Hearing Outcomes in Simultaneous Bilateral Cochlear Implanted Adults? 38 th ARO Meeting Abstract[END_REF] . To be implanted, patients were required to have a maximum of 10% open set disyllabic word recognition score in quiet at 60 dB in the bestaided condition, a difference of profound hearing loss duration between the two ears of less than 5 years, and no malformations of the cochlea. Speech perception tests in quiet and in noise (SNR of +15 dB, +10 dB and +5 dB) were performed before implantation, at 1-year, and 5-years after activation. Responses were scored as the percentage of words correctly identified. All patients underwent bilateral implantation by expert otologists (more than 100 CI procedures) in a simultaneous surgical procedure with the same device (MED-EL Combi 40+, Standard Electrode Array, 31 mm length; Innsbruck, Austria).

A multi-slice helical CT scan (500 µm slice thickness), was realized in the immediate postoperative period.

All participants gave their informed written consent, and the study was approved by the local ethical committee (Saint-Louis, Paris, No. 61D0/22/A).

Radiological analysis

The DICOM (Digital Imaging and Communications in Medicine) data were analyzed by Osirix program (Osirix v 4.0 64-bit; Pixmeo Sarl, Bernex, Switzerland). This program allowed multiplanar reconstructions of cochlear anatomy and position of the arrays in the cochlea. All the images, acquired by different CT scans in the different centers, were reconstructed with 0.1 mm increments in order to standardize the measurement technique and reduce the error of measurement. To examine the cochlear sizes and their relationship with the insertion depth, a three-dimensional coordinate system was used, in accordance with the consensus of cochlear coordinates [START_REF] Verbist | Consensus panel on a cochlear coordinate system applicable in histologic, physiologic, and radiologic studies of the human cochlea[END_REF] , with the exception of the cochlear height that was measured in a reformatted coronal view. The largest cochlear diameter (distance A) going from the center of the round window membrane to the opposite lateral wall [START_REF] Escudé | The size of the cochlea and predictions of insertion depth angles for cochlear implant electrodes[END_REF] , was calculated on a plane perpendicular to the modiolus axis and coplanar to the basal turn, named 'cochlear view' by Xu et al. [START_REF] Xu | Cochlear view: postoperative radiography for cochlear implantation[END_REF] (Fig. 1A). The cochlear height was measured from the mid-point of the basal turn to the mid-point of the apical turn on a coronal section [START_REF] Purcell | Establishment of normative cochlear and vestibular measurements to aid in the diagnosis of inner ear malformations[END_REF][START_REF] Mori | CT analysis demonstrates that cochlear height does not change with age[END_REF] (Fig. 1B). The electrode-to-modiolus distances (EMD) for electrodes positioned at 180-and 360-degrees were measured on the mid-modiolar plane, crossing the mid of the round window (Fig. 1C).

The angular depth of insertion of the array was measured in the 'cochlear view' (slice thick of 5 mm), considering the mid-point of the round window as the 0-degrees reference (Fig. 1D).

To minimize the error, all the measurements were performed blindly by an otologist, each measurement was repeated three times in nonconsecutive days, and the mean value was then considered.

Statistical analysis

Values are expressed as means  standard error of the mean (SEM). For correlations between cochlear anatomy and cochlear array localization, and its relation with speech perception scores, Pearson's correlation coefficient (r) was calculated, and the ANOVA was used to test the slope of the linear regression line. One-way ANOVA was used to analyze the influence of the number of activated electrodes on speech performance. Student's t-test was used for comparisons between groups (male/female, right/left cochleae, full/partial insertions). For all comparisons, p<0.05 was considered as significant. All statistical analyses were performed using IBM SPSS for Windows (v 22.0, SPSS Inc., Chicago, Illinois, USA).

RESULTS

The mean speech performance in quiet and noise have previously been reported [START_REF] Mosnier | Speech performance and sound localization in a complex noisy environment in bilaterally implanted adult patients[END_REF] . At 1-year post-implantation, 7 patients were poor performers (speech perception scores in quiet < 60% in bilateral condition). Among the good performers, 9 patients obtained asymmetrical performance (difference of speech scores in quiet between the two ears ≥20%).

Cochlear anatomy and electrode position

The cochlear anatomical data are reported in Table 2. The distance A was positively correlated with the cochlear height (r=0.52, p=0.0007, data not shown). Surprisingly, the distance A and the cochlear height were different between the two ears (difference of mean distance A: 0.220.05 mm, p=0.04; difference of mean cochlear height: 0.30.06 mm, p=0.001, Student's t tests); no right or left ear predominance was observed. The distance A and the cochlear height were different as well between male and female ears, the males having a diameter and a cochlear height greater than females (p=0.0001, Student's t tests).

A full insertion of the electrode array was achieved in 26 ears, and a partial insertion in 12 ears (3 patients with a bilateral partial insertion, and 6 patients with a unilateral partial insertion). In ears with an incomplete insertion, the number of extra-cochlear electrodes ranged from 1 to 4. The size of the cochlea (i.e. distance A and cochlear height) was similar between the ears with a full insertion and ears with a partial insertion (Table 3).

In the 26 ears with a full electrode insertion, the angular depth of insertion in the cochlea varied widely [510-880-degrees] (Fig 2 ), and was negatively correlated with the distance A (r=-0.55, p=0.003) (Fig. 3A), on the other hand no correlation was found with the cochlear height (Fig 3B). The EMD was positively correlated with the distance A at both 180-(r=0.47, p=0.0004) and 360-degrees (r=0.66, p=0.0002, Fig. 3C), and with the cochlear height at 360degrees (r=0.6, p=0.001, Fig. 3D). The EMD distance at 180-and at 360-degrees was not correlated with the angular depth of insertion. These results indicate that in large cochleae (distance A), the electrode array was less deeply inserted and more distant from the modiolus at the basal turn (EMD at 180-degrees and 360-degrees). In the present study, the distance A was sufficient to define the cochlear size and reliable for the prediction of the position of the implant within the cochlea.

Correlation between electrode position and speech perception

At 1-year after cochlear implantation (38 implanted ears), speech perception scores were negatively correlated with EMD at 180-degreess both in quiet (r=-0.34, p=0.02) and in noise (SNR +15 dB: r=-0.44, p=0.006; SNR +10 dB: r=-0.63, p=0.0005; SNR+5 dB: r = -0.52, p=0.01, Fig. 4). The greater the EMD was, the poorer was the performance. No correlation was observed at 360-degrees. The number of inserted electrodes was correlated with speech perception in noise at SNR +15 dB and SNR +10 dB (ANOVA, p=0.02); the speech perception scores in noise gradually decreased as a function of the number of inserted electrodes (post hoc Dunnett's t test p=0.02) (Table 3). Considering the obvious interdependence between the number of intracochlear electrodes and the depth of insertion, we analyzed the influence of electrode position on hearing outcomes among the 26 ears with a full insertion of the electrode array. No correlation was found between the speech perception scores and the angular depth of insertion, both in quiet and in noise, whereas the speech perception scores were negatively correlated with EMD at 180-degreess both in quiet (r=-0.38, p=0.048) and in noise (SNR +15 dB: r=-0.4, p=0.049; SNR +10 dB: r = -0.62, p=0.006; SNR+5 dB: r=-0.51, p=0.032, data not shown).

A multifactorial ANOVA was performed and failed to demonstrate that the anatomic cochlear variations (distance A, cochlear height), and the different electrode position (EMD at 180-and 360-degrees) between the two ears, were the reason of the asymmetric speech score (difference ≥20% between better and poorer ear) at 1-year in 9 patients.

At 5-years post-implantation, most of the patients (85%) achieved good speech performance (speech perception score ≥60% in quiet in bilateral condition); the speech score of the poorer ear in noise continued to improve over time, and the majority of the patients with poor speech scores improved their performance both in quiet and in noise [START_REF] Seta | Does the Anatomy of the Cochlea influences the Hearing Outcomes in Simultaneous Bilateral Cochlear Implanted Adults? 38 th ARO Meeting Abstract[END_REF] . Studying the relationship between the electrode insertion parameters and the hearing outcomes, no correlation was found at 5-years postimplantation between speech perception scores and the angular depth of insertion, both in the entire sample and in the group with full insertion of the electrode array.

In contrast to what observed at 1-year postimplantation, the EMD was not correlated with speech perception scores, both at 180-degrees and 360-degrees (data not shown).

DISCUSSION

We have previously shown that in adult patients simultaneously and bilaterally implanted, poor or asymmetrical hearing performance at 1-year postimplantation are present in 40% of cases, and that the speech scores of the poorer ear continues to improve over time [START_REF] Mosnier | Speech performance and sound localization in a complex noisy environment in bilaterally implanted adult patients[END_REF][START_REF] Seta | Does the Anatomy of the Cochlea influences the Hearing Outcomes in Simultaneous Bilateral Cochlear Implanted Adults? 38 th ARO Meeting Abstract[END_REF] . In the present study, we demonstrate that both the distance between electrode array and modiolus at 180-degrees, and the number of inserted electrodes, are important variables that influence the early achievement of the best speech perception scores. The variability in cochlear anatomy could explain the differences in hearing outcomes between patients; nevertheless we failed to demonstrate an influence of cochlear geometry on intraindividual speech perception asymmetry, probably due to the small number of patients with asymmetric speech scores.

The variability in cochlear anatomy influences electrode array position

Several studies investigated the influence of cochlear anatomy on electrode array position within the cochlea [START_REF] Kawano | Computer-aided three-dimensional reconstruction in human cochlear maps: measurement of the lengths of organ of Corti, outer wall, inner wall, and Rosenthal's canal[END_REF][START_REF] Medina | Is electrode-modiolus distance a prognostic factor for hearing performances after cochlear implant surgery?[END_REF][START_REF] Verbist | Anatomic considerations of cochlear morphology and its implications for insertion trauma in cochlear implant surgery[END_REF][START_REF] Van Der Marel | Diversity in cochlear morphology and its influence on cochlear implant electrode position[END_REF][START_REF] Franke-Trieger | Insertion depth angles of cochlear implant arrays with varying length: a temporal bone study[END_REF] . Important variations in the first segment of the scala tympani, such as unusual narrowing or constriction, have been reported. The basal end of the cochlea is in fact of major interest in cochlear implant surgery; it bends in three dimensions, resembling to a ''fish hook'', and in some cases its anatomical variations lead to a difficulty for the surgeon to choose the ideal cochleostomy site in order to reach the scala tympani without damaging any inner ear structure [START_REF] Rask-Andersen | Human cochlea: anatomical characteristics and their relevance for cochlear implantation[END_REF] .

In this study the cochlear size was assessed using the major cochlear diameter of the basal turn, that is assumed to be a good predictor of the length of the two first turns of the cochlea [START_REF] Erixon | How to predict cochlear length before cochlear implantation surgery[END_REF][START_REF] Singla | Surgical anatomy of the basal turn of the human cochlea as pertaining to cochlear implantation[END_REF][START_REF] Alexiades | Method to Estimate the Complete and Two-Turn Cochlear Duct Length[END_REF] , and using the cochlear height; our results are in line with the data present in literature [START_REF] Erixon | Variational anatomy of the human cochlea: implications for cochlear implantation[END_REF][START_REF] Martinez-Monedero | Cochlear coiling pattern and orientation differences in cochlear implant candidates[END_REF][START_REF] Escudé | The size of the cochlea and predictions of insertion depth angles for cochlear implant electrodes[END_REF][START_REF] Purcell | Establishment of normative cochlear and vestibular measurements to aid in the diagnosis of inner ear malformations[END_REF][START_REF] Mori | CT analysis demonstrates that cochlear height does not change with age[END_REF] . These two measures are clearly correlated to each other, meaning that a greater basal turn diameter is associated to a higher cochlea. Both distance A and cochlear height vary with sex, males having bigger cochlea compared to females, as already described in the literature [START_REF] Escudé | The size of the cochlea and predictions of insertion depth angles for cochlear implant electrodes[END_REF][START_REF] Mori | CT analysis demonstrates that cochlear height does not change with age[END_REF][START_REF] Van Der Marel | Diversity in cochlear morphology and its influence on cochlear implant electrode position[END_REF] . Additionally, we observed an asymmetry between the two ears in distance A (0.22 mm), that was only described by Escude et al. [START_REF] Escudé | The size of the cochlea and predictions of insertion depth angles for cochlear implant electrodes[END_REF] , and in cochlear height (0.3 mm). No ear predominance was found, as previously reported [START_REF] Mori | CT analysis demonstrates that cochlear height does not change with age[END_REF][START_REF] Van Der Marel | Diversity in cochlear morphology and its influence on cochlear implant electrode position[END_REF][START_REF] Singla | Surgical anatomy of the basal turn of the human cochlea as pertaining to cochlear implantation[END_REF][START_REF] Alexiades | Method to Estimate the Complete and Two-Turn Cochlear Duct Length[END_REF][START_REF] Pelliccia | Cochlea size variability and implications in clinical practice[END_REF] .

In patients implanted with long (31 mm) and straight electrode arrays, we demonstrated that as expected, the smaller the diameter of the cochlea is, the closer is the electrode array to the modiolus at the basal turn, and the deeper is the array insertion. The depth of array insertion was strongly correlated (r=-0.63) with the major cochlear diameter measurement, with a shallower insertion in bigger cochlea and deeper insertion in smaller cochlea. Van der Marel et al. [START_REF] Van Der Marel | Diversity in cochlear morphology and its influence on cochlear implant electrode position[END_REF] found a weaker correlation (Pearson's r = -0.3) analyzing 362 cochleae implanted with Advanced Bionics implants. In other studies, a more significant correlation between depth of insertion and cochlear diameter was found using straight electrodes [START_REF] Franke-Trieger | Insertion depth angles of cochlear implant arrays with varying length: a temporal bone study[END_REF][START_REF] Franke-Trieger | Estimation of insertion depth angle based on cochlea diameter and linear insertion depth: a prediction tool for the CI422[END_REF] .

An incomplete insertion of the electrode array was observed in 12/38 ears (32%). This observation is in accordance with a histopathological study, which reported the 52% of incomplete insertion in absence of intrascalar bony or soft tissue that could explain a partial insertion [START_REF] Lee | Factors associated with incomplete insertion of electrodes in cochlear implant surgery: a histopathologic study[END_REF] . The anatomical study of Rask-Andersen et al. [START_REF] Rask-Andersen | Human cochlea: anatomical characteristics and their relevance for cochlear implantation[END_REF] describes a narrowing of the cochlear duct or a sharp bend of cochlear coiling between the first and the second turn as another possible cause for incomplete insertion. No significant difference in the size of the cochlea between ears with incomplete and complete insertions was found in our study, nevertheless it should be noticed that the three cochleae with 4 electrodes outside, had a smaller distance A than the other ears (see Table 3). On the base of the cochlear length equation based on distance A value (Alexiades et al. [START_REF] Alexiades | Method to Estimate the Complete and Two-Turn Cochlear Duct Length[END_REF] ), we can assume that a 31 mm length array was too long to be totally inserted in these three ears. At the present, different lengths of cochlear arrays are available, and it is crucial to measure the distance A before implantation in order to adapt the type (and length) of the electrode array to be implanted.

Is the electrode position related to speech perception?

If we consider the ears with full insertion of the electrode array, despite a large variation of the angular depth of insertion, no correlation was found between this variable and the hearing performance. This observation is consistent with a histological analysis over a series of 27 temporal bone specimens of subjects with cochlear implant [START_REF] Lee | Depth of electrode insertion and postoperative performance in humans with cochlear implants: a histopathologic study[END_REF] . Van der Marel et al. [START_REF] Van Der Marel | The Influence of Cochlear Implant Electrode Position on Performance[END_REF] analyzed six position-related variables including the angular and linear insertion depth of the array and did not find any correlation with speech outcomes at 2-years postoperative. In a prospective randomized study including 13 patients, Buchman et al. [START_REF] Buchman | Pillsbury HC Influence of cochlear implant insertion depth on performance: a prospective randomized trial[END_REF] didn't find a difference in speech scores between MedEl standard array (mean angular depth of insertion 657-degrees) and medium array (mean angular depth of insertion 423-degrees), although better performance was found in the standard array group when 6 more patients were included retrospectively. On the contrary, other studies reported poorer performance in case of deeper insertions [START_REF] Skinner | In vivo estimates of the position of advanced bionics electrode arrays in the human cochlea[END_REF] , explained by the increased number of electrodes in the scala vestibuli, reduced pitch discrimination, decreased basal stimulation [START_REF] Finley | Role of electrode placement as a contributor to variability in cochlear implant outcomes[END_REF] , and pitch confusion at apical contacts [START_REF] Gani | Implications of deep electrode insertion on cochlear implant fitting[END_REF] . The negative correlation between the electrode angular depth of insertion and hearing outcomes found by Yukawa et al. [START_REF] Yukawa | Effects of insertion depth of cochlear implant electrodes upon speech perception[END_REF] may be explained by the presence of confounding factors, such as the lower number of activated electrodes in case of partial insertion. Indeed, in the present study, in case of incomplete insertion, the speech perception scores in noise at 1-year decreased as a function of the number of inserted electrodes (see Table 3).

Considering the distance between the electrode array and the modiolus, it has been shown that a closer position to the spiral ganglion cells is associated with better speech perception [START_REF] Medina | Is electrode-modiolus distance a prognostic factor for hearing performances after cochlear implant surgery?[END_REF][START_REF] Finley | Role of electrode placement as a contributor to variability in cochlear implant outcomes[END_REF] . This effect may be related to the minimization of channel interaction, which leads to reduction of electrical thresholds and/or improvement of the spatial selectivity. Our findings are in accordance with Esquia-Medina et al. [START_REF] Medina | Is electrode-modiolus distance a prognostic factor for hearing performances after cochlear implant surgery?[END_REF] who reported a correlation between speech perception scores and average EMD of the 6 most basal electrodes of MED-EL devices (corresponding approximately to the region from 0-to 180-degrees) at 6 months, whereas no correlation was found at 12 months. In this study, as well as the present one, such relationship was not present for the electrode at 360-degrees, possibly due to the narrowing of the scala tympani from base to apex [START_REF] Biedron | The internal dimensions of the cochlear scalae with special reference to cochlear electrode insertion trauma[END_REF] that reduces the variability of the array position. This relationship between the EMD and the hearing performance could point out a preferential use of perimodiolar electrode array in order to obtain a rapid hearing rehabilitation. Nevertheless, Doshi et al. [START_REF] Doshi | Straight versus modiolar hugging electrodes: does one perform better than the other?[END_REF] reported no differences between speech perception outcomes at 3-and 9-months in patients implanted with either straight or perimodiolar electrodes array. A reason could be the more frequent dislocation from scala tympani to scala vestibuli in case of perimodiolar electrodes [START_REF] Boyer | Scalar Localization by Cone-Beam Computed Tomography of Cochlear Implant Carriers: A Comparative Study Between Straight and Periomodiolar Precurved Electrode Arrays[END_REF] .

Although such scalar dislocation is difficult to assess in standard CT scan, it might negatively influence the cochlear implant outcome [START_REF] Aschendorff | Quality control after insertion of the nucleus contour and contour advance electrode in adults[END_REF][START_REF] Holden | Factors affecting open-set word recognition in adults with cochlear implants[END_REF][START_REF] Gani | Implications of deep electrode insertion on cochlear implant fitting[END_REF][START_REF] Wanna | Impact of electrode design and surgical approach on scalar location and cochlear implant outcomes[END_REF] . An aspect that has not been explored in this study is the surgeon's gesture. A recent study described a high intra-and inter-individual variability of the insertion axis of the array into the cochlea; yet, this variability was reduced among expert surgeons [START_REF] Torres | Variability of the mental representation of the cochlear anatomy during cochlear implantation[END_REF] . Since all the participants to the present study were senior otologists, we estimate that this doesn't represent a great factor of bias of the study; furthermore, how the insertion axis influences the trajectory of insertion or the final position of the array has not yet been described or reported. An additional limitation of this study could be represented by the migration of the array that can occur between 1-and 5-years. Nevertheless, in all patients the most basal electrodes remained activated with stable impedance values over time and providing auditory responses, thus an extrusion of the electrodes should be unlikely [START_REF] Johnston | Computed Tomography Estimation of Cochlear Duct Length Can Predict Full Insertion in Cochlear Implantation[END_REF] .

In conclusion, whereas 1-year results suggest that the number of inserted electrodes and the distance electrode-to-modiolus are related to good performance, these parameters does not influence the speech scores after long term use. In order to obtain a rapid hearing rehabilitation and the best results at 1-year, the preoperative measurement of the cochlear diameter (distance A) may guide the choice of the correct array length allowing a complete insertion. In case of unilateral implantation the choice of the side to be implanted should be oriented, in presence of equal clinical and audiological conditions of the two ears, to the smaller cochlear diameter. The lines represent the significant linear regression.

FIGURES

Figure 1 :DFigure 2 :

 12 Figure 1:

Figure 3 :

 3 Figure 3: Correlation between the size of the cochlea (cochlear diameter, cochlear height) and

Figure 4 :Table 2 :Female (n = 28 )

 4228 Figure 4: Correlations between the electrode array position and the speech perception scores

Acknowledgments

DDS would like to thank the French Society of Otolaryngology (SFORL) for the 2013 Research Grant.

The authors report no conflicts of interest.


Frédéric Perales
email: frederic.perales@irsn.fr

Yann Monerie

Frédéric Dubois

Laurent Stainier

Application to fracture of hydrided Zircaloy

Keywords: Fissuration, non smooth dynamics, cohesive zone model, numerical simulation, hydrided Zircaloy-4

scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

INTRODUCTION

Zirconium alloys like Zircaloy-4 are metals used as the cladding material of nuclear fuel rods in Pressurized Water Reactor (PWR). During nuclear reactor operation, bulk microstructure of these metals evolves because of hydrogen which precipitates as zirconium hydrides. So, the microstructure of fuel cladding tubes appears as a bimaterial constituted of hydride embedded in a metal matrix. To analyse the effects of this microstructure heterogeneity on material behavior, we use a computational micromechanical approach.

The micromechanical approach considered in this work is based on the concept of Frictional Cohesive Zone Model (FCZM) which represents the mechanical processes in the fracture process zone in front of a crack tip. The FCZM allows simulation of crack initiation (without any ad hoc criteria), crack propagation (crack path is not known in advance) and non-smooth surfacic behavior such as unilateral contact and friction, in heterogeneous materials at several scales. The mechanical separation is described by a constitutive equation relating applied stress on the crack lips to the displacement jump between each element of a finite element modeling. Each mesh is, thus, considered as an independant body connected to another by non-linear and non-smooth relations.

Since velocity can become discontinuous during dynamic fracture, the Non-Smooth Contact Dynamic (NSCD) approach is used. This approach, initiated by [START_REF] Moreau | Non Smooth Mechanics and Applications[END_REF] and Jean (1999), consists in solving nonsmooth equations of frictional contact without any regularization nor penalization.

The numerical platform, developed by the French 'Institut de Radioprotection et de S ûreté Nucléaire' (IRSN) in the frame of its research program on nuclear fuel safety, is composed of three libraries : a Fortran90 library dedicated to surfacic behaviors related to FCZM in the NSCD approach (LMGC90), a C++ library dedicated to Finite Element modeling (PELICANS) and a C++ library dedicated to bulk constitutive models (MATLIB).

The ability of the numerical platform is illustred on the dynamic fracture of hydrided Zircaloy-4 under transient loading. A study is performed to determine the macroscopic behavior of hydrided Zircaloy-4 as function of hydrogen content.

FRICTIONAL COHESIVE ZONE MODEL (FCZM)

The FCZM considered here is a model coupling cohesion with unilateral contact and Coulomb friction inspired by the cohesive law of [START_REF] Raous | [END_REF]; [START_REF] Raous | 3d Contact Mechanics International Symposium, Collection Solid Mechanics and its Applications[END_REF]. A cohesive reaction, denoted R coh , is added to the unilateral contact problem (Eq. 1) and the Coulomb friction problem (Eq. 2).

0 ≤ (R N + R coh N )⊥u N ≥ 0 (1) R T + R coh T < µ|R N + R coh N | ⇒ uT = 0 R T + R coh T = µ|R N + R coh N | ⇒ ∃λ ≥ 0 such as uT = -λ(R T + R coh T ) (2)
Subscripts N and T respectively denote the normal and tangential component of the total force (R = R N n + R T ), the cohesive force (R coh = R coh N n + R T ) and the displacement jump ([u] = u N n + u T ), n is the unit normal vector of the cohesive zone. The friction coefficient is denoted µ.

In this study, in order to particularize the FCZM, we use the cohesive force described by Eq. 3, initially introduced by [START_REF] Fremond | [END_REF], based on a surfacic variable β (Eq. 4) (β = 1 : the cohesion is complete, 0 < β < 1 : the cohesion is partially broken and β = 0 : there is no cohesion).

R coh = K(β) • [u] if β > 0 or u N = 0 K(β) • [u] -p • n otherwise (3) β = min(g( [u] ), g( [u max ] )) (4)
where u max is the maximum value reached by u, K(β) is a matrix describing the progressive surfacic softening behavior during crack propagation, p denotes internal post-fracture pressure (MPa) and g is a decreasing function of [u] .

We use a variant of the surfacic damage law defined by [START_REF] Michel | Revue Europenne des Elements Finis[END_REF] : andC T are the initial normal and tangent stiffness to the interface if cohesion is complete (MPa/m), 0 ≤ β 0 ≤ 1 is a initial surface damage and w is a reference fracture energy (J/m 2 ). In a 2D case, Figure 1 shows the normal behavior (with u T = 0) and the tangential behavior (with R N constant).

K(β) = β C N n ⊗ n + C T u T ⊗ u T u T 2 g(x) =        β 0 if x ≤ δ 0 0 if x ≥ δ c β 0 (δ c -x) δ c -2δ 0 + x otherwise where δ 0 = w (9 -4ln4) 1 C N + 1 C T , δ c = 3δ 0 , C N
PSfrag replacements

R N u N /δ 0 p R T u T µR N (a) (b) 1 σ 0 0 PSfrag replacements R N u N /δ 0 p R T u T µR N (a) (b) 1 σ 0 0 Figure 1: 2D cohesive law : (a) normal behavior (u T = 0) (b) tangent behavior (R N constant and u N = 0).
When decohesion vanishes (β = 0 in this model), we obtain the classical Signorini problem for the normal behavior and the Coulomb friction problem for the tangential behavior.

NON-SMOOTH CONTACT DYNAMICS (NSCD)

The equation of dynamics

The dynamic problem is composed by the Eqs. 1-4 and the discrete dynamic equation :

M q = F (q, q, t) + r (5)
where q, q and q respectively are discrete displacement, velocity and acceleration, F (q, q, t) represents internal and external forces and r is representative of local reaction forces. Since shocks are expected, the velocity q can be a step function and the acceleration q, then, cannot be defined. In this way, Eq. 5 is solved using NSCD approach (Jean, 1999) within the functional framework proposed by [START_REF] Moreau | Non Smooth Mechanics and Applications[END_REF]. Particularly, the derivative in the dynamic equation are to be understood in the sense of distributions.

Time discretization

Denoting by a subscript i quantities at a time t i and by i + 1 quantities at time t i+1 , the numerical method (θ-method) can be written over the time interval ]t i , t i+1 ] as :

M i+1 ( qi+1 -qi ) = h[(1 -θ)F i + θF i+1 ] + hr i+1 q i+1 = q i + h[(1 -θ) qi + θ qi+1 ] (6)
where hr i+1 is the mean impulse value. The unknowns of the problem are the approximations of the velocity qi+1 and the impulse hr i+1 . A Newton-Raphson algorithm is used to solve the non-linear set of discretized equations ( 6) :

qk+1 i+1 = qk i+1 + ∆ qfree + hw k r k+1 i+1 ∆ qfree = w k i+1 {-M k i+1 ( qk i+1 -qi ) + h[(1 -θ)F i + θF k i+1 ] + hr k i+1 } (7)
where k stands for iterations, w k i+1 is the inverse of the iteration matrix.

Non-Smooth Contact Dynamics and Frictional Cohesive Zone Model

Using linear mappings between global values (velocity q and impulse hr) and their relative values ( u and hR), and a change of variable introduced by Monerie and Acary (2001), the three dimensional system to solve for each contact α writes (Jean et al., 2001) :

     U α -U α locf ree + W αα h(R -R coh ) α = 0 R α N -proj R -(R α N + ρU α N ) = 0 R α T -proj D(µ|R α N |) (R α T + ρU α T ) = 0 (8)
where ρ > 0 and D(µ|R α N |) is the section of Coulomb's cone for contact α i.e. the disc with center 0 and radius µ|R α N | and U α locf ree is the velocity at the contact α without any reaction coming from the FCZM. The system (8) may be written as a mapping φ :

φ(U α , R α ) = 0 (9)
and solved using the generalized Newton method of Alart and Curnier (1988).

NUMERICAL PLATFORM

The numerical platform is composed of three libraries with Object Oriented Programming (OOP) concepts. This design increases software flexibility (reusability makes the existing code easier to modify and extend its functionality and capabilities). The three libraries are coupled with an encapsulated design and each module has a clear meaning from a mechanical point of view :

• LMGC90 is a Fortran90 library dedicated to surfacic behaviors related to FCZM in the NSCD approach. It is developed at the Laboratory of Mechanics and Civil Engineering by [START_REF] Dubois | Actes du sixième colloque national en calcul de structures -CSMA-AFM-LMS[END_REF].

• PELICANS is a C++ Finite Element library developed by the French 'Institut de Radioprotection et de S ûreté Nucléaire' (IRSN) [START_REF] Piar | Trend in Physical and Numerical of Multiphase Industrial Flows[END_REF].

• MATLIB is a C++ library which allows to consider constitutive models written in the standard generalized thermodynamics framework as well as user-defined constitutive models. It is developed at the Department of Aerospace, Mechanics and Materials by [START_REF] Stainier | 6ème Colloque National en Calcul des Structures[END_REF].

The main idea is to include cohesive surfaces along all finite element boundaries [START_REF] Xu | [END_REF]. So, a multibody strategy is chosen, each finite element being a body (see Figure 2).

The software allows three dimensional finite deformation simulations in heterogeneous materials from crack initiation (without initiation criteria) to post-fracture behavior.

APPLICATIONS

Dynamic crack propagation in a representative volume element (RVE) of hydrided Zircaloy-4

Zircaloy cladding tubes at high burn-up are considered as a bimaterial which is formed by hydride inclusions (hydride content varies gradually in the tube thickness) surrounded by a metallic matrix. Two cases are considered : a RVE of Zircaloy-4 without hydrides and a RVE of Zircaloy-4 with hydrides distributed randomly (with a volume fraction of 30 % ∼ 5473 ppm) and oriented horizontally. The RVE is a square (length of a side = 10 µm). Matrix is assumed to be elastoplastic (J 2 plasticity) and hydrides neo-Hookean. The bonding strength value between the two phases is considered as "strong" (w i /w h 1, w i and w h denote the surface energy of interface matrix/hydride and matrix/matrix respectively). The material properties and the cohesive coefficients are given in Table 1 and2 Table 1: Material properties used in the finite element calculations. The material properties correspond to irradiated Zircaloy-4 matrix and δ-hydride.

C N = C T (P a/m) w (J/m 2 ) Matrix/Matrix 10 19 1 Matrix/Hydride 2x10 19 1 Hydride/Hydride 2x10 19 1
Table 2: Constitutive parameters for cohesive surfaces.

On the "strength vs imposed displacement" curve (in the left in Figure 3), the volume behavior (elastoplastiticy) associated with the FCZM (surfacic behavior described by a softening function) lead to a damage elastoplasticity behavior for the Zircaloy-4.

Figure 3 shows the weakening effect of the hydrides. The introduction of 30 % volume fraction of hydrides reduces the total elongation (TE) by a factor of 5 (TE for Zircaloy-4 24 %, TE for hydrided Zircaloy-4 4.8 %). In this case, the behavior of the hydrided Zircaloy-4 is brittle. These results are consistent with experimental results.

Hydrides have a weakening effect : the energy release for hydrided Zircaloy-4 is approximately 80 % lower than for the non hydrided material.

Macroscopic behavior of RVE of hydrided Zircaloy-4 depending on hydrogen content

This study deals with the macroscopic description of damage in hydrided Zircaloy-4 with various hydrogen contents. For that, we consider RVEs with different microstructures. Each microstructure is composed of a metallic matrix (Zircaloy-4) with rectangular elastic inclusions (hydrides) at a given volume fraction. RVEs are sufficiently large to ascertain statistical uniformity of distributed hydrides (square : length of a side = µm).

As the example 5.1, the matrix is assumed to be elastoplastic (J 2 plasticity), hydrides neo-Hookean and the bonding strength value between the two phases "strong". The material properties and cohesive coefficients respectively are given in Table 1 and2 In Figure 4, the weakening effect is accuratly represented by Eq. 10.

W W zr = 1 + W i W zr f 1 + W i W zr -1 f (10)
where W , W Zr and W i respectively are the macroscopical energy release rate in hydrided Zircaloy-4, Zircaloy-4 and hydrides (J), f = 91 1.66x10 6 [H] is the hydride volume fraction (%) and [H] is the hydrogen content (ppm).

CONCLUSION

This paper has presented a numerical framework for the simulation of the dynamic fracture of heterogeneous materials like hydrided Zircaloy-4. This approach is based on the coupling between the Cohesive Zone Model concept and the Non-Smooth Contact Dynamics method. The associated software is developed by coupling libraries with the Object Oriented method. The ability of this software has been illustrated on the dynamic crack propagation of hydrided Zircaloy-4. A damage elastoplastic behavior has been simulated for the Zircaloy-4 and the weakening effect of the hydrides has been emphasised. A study for the macroscopic behavior of hydrided Zircaloy-4 versus hydrogen content has been leaded.

Figure 2 :

 2 Figure 2: Multibody strategy : each finite element is a body

  . We assume there is low friction µ = 0.05, no post-fracture pressure p = 0 and same compliance for the normal and tangential behaviors C N = C T . We impose velocity boundary conditions along horizontal faces : V = 2 m/s (right) and V = -2 m/s (left).Young Modulus Poisson's ratio Yield in tension Hardening Modulus (

Figure 3 :

 3 Figure 3: Strength vs imposed displacement for Zircaloy-4 and hydrided Zircaloy-4 (5473 ppm) and rupture features of a RVE of hydrided Zircaloy-4.

Figure 4 :

 4 Figure 4: Influence of hydrogen content on energy release rate.

Figure 4

 4 Figure4shows the energy release rate versus hydrogen content. The energy release rate decreases rapidly with hydrogen content (180 ppm of hydrogen reduces the energy release rate by 20 %, 920ppm by 40 %, 7296 ppm by 80 %). These results show the importance of the hydride and, in particular, their deleterious influence on the material ductility.In Figure4, the weakening effect is accuratly represented by Eq. 10.