Raphaël Cerf

ASYMPTOTIC CONVERGENCE OF GENETIC ALGORITHMS

Keywords: 1991 Mathematics Subject Classification. 60F10, 60J10, 92D15 Freidlin-Wentzell theory, genetic algorithms, stochastic optimization

We study a markovian evolutionary process which encompasses the classical simple genetic algorithm. This process is obtained by perturbing randomly a very simple selection scheme. Using the Freidlin-Wentzell theory, we carry out a precise study of the asymptotic dynamics of the process as the perturbations disappear. We show how a delicate interaction between the perturbations and the selection pressure may force the convergence toward the global maxima of the fitness function. We put forward the existence of a critical population size, above which this kind of convergence can be achieved. We compute upper bounds of this critical population size for several examples. We derive several conditions to ensure convergence in the homogeneous case: these provide the first mathematically wellfounded convergence results for genetic algorithms.

I. Introduction

Stochastic optimization has been a field of intense research during the last decade: computers become more and more powerful and there is an increasing need for robust and efficient search techniques. In the jungle of optimization procedures, there exists a particular class of algorithms called genetic algorithms. Introduced by J.H. Holland in the seventies [START_REF] Holland | Adaptation in natural and artificial systems[END_REF], these algorithms are based on the genetic mechanisms which guide natural evolution: mutation, crossover and selection. Today, they are widely used to handle a large spectrum of optimization problems ranging from the classical traveling salesman problem to the design of network architecture [START_REF] Goldberg | Genetic algorithms in search, optimization and machine learning[END_REF].

However the behavior of genetic algorithms relies heavily on several control parameters (for instance the population size, the probabilities of mutation and crossover). Until now, no firm result concerning the choice of these parameters and the dynamics of the genetic algorithm has been available: practical know-how and experimental simulations were the only guides for handling concrete problems.

In this paper, we study a slightly more general evolutionary scheme. We consider first a very simple selection algorithm which is extremely rapidly trapped in bad points of the search space and we perturb randomly the whole mechanism. We focus on the asymptotic dynamics of the process as the random perturbations vanish. We show how a delicate interaction between mutations, crossovers and selection pressure ensures the convergence of the algorithm toward the set of the global maxima of the fitness function. We put forward the existence of a critical population size, above which this kind of convergence can be achieved. Our results seem to be the first well-founded convergence theorems of this kind concerning genetic algorithms.

To fulfill this program, we use the powerful tools developed by Freidlin and Wentzell in a much more general framework for the study of random perturbations of dynamical systems [START_REF] Freidlin-A | Random perturbations of dynamical systems[END_REF].

Classical genetic algorithms are not usually described as the random perturbation of a simple process. For a fixed level of intensity of the random perturbations, our model evolves exactly as a classical genetic algorithm. We are only able to analyze rigorously the asymptotic regime when the perturbations vanish. On one hand, the asymptotic picture sheds some light on the true dynamics of classical genetic algorithms. On the other hand, this approach suggests new ways of implementing genetic algorithms, using the theory of generalized annealing processes [START_REF] Deuschel | L 2 convergence of time nonhomogeneous Markov processes: I. Spectral Estimates[END_REF][START_REF] Trouvé | Asymptotical behaviour of several interacting annealing processes[END_REF].

To avoid lengthy and technical computations we focus here mainly on the study of the invariant measure of the process and we just quote without proof a result of convergence for the inhomogeneous case. In another work, we apply to our model the very technical tools developed by Catoni and Trouvé to control the speed of convergence of generalized annealing processes to equilibrium (see [START_REF] Cerf | The dynamics of mutation-selection algorithms with large population sizes[END_REF] and the references therein). A major difficulty with such dynamics is the lack of reversibility [START_REF] Deuschel | L 2 convergence of time nonhomogeneous Markov processes: I. Spectral Estimates[END_REF]. This paper has the following structure. We first describe our model. Then, we study the convergence of the algorithm in the homogeneous case and derive sufficient conditions to ensure that the population settles in a global maximum of the fitness function. We show that these conditions are fulfilled when the population size is large enough. We apply our results to the classical simple genetic algorithm. In the inhomogeneous case, we give several conditions on the rate of decrease of the perturbations to ensure the convergence of the population to the set of global maxima in finite time. We finally examine the value of the critical population size for a specific problem coming from statistical mechanics, namely the search of the ground state of the Ising model in a finite box.

II. General notations and conventions

The cardinality of a set X will be noted indifferently |X| or card X and its characteristic function 1 X . We adopt usual conventions concerning empty sets:

min ∅ = +∞, max ∅ = -∞, ∅ = 1, ∅ = 0.
If s is a real number, s denotes the unique integer such that s ≤ s < s + 1. We say that a real number s is positive if s > 0 and that it is non-negative if s ≥ 0. We consider a finite space of states E and a real-valued positive non-constant function f (which will be called the fitness function) defined on E. The letters i, j will denote elements of E. The set of global maxima of f is

f * = { i ∈ E : f (i) = max j∈E f (j) }.
By f (f * ) we mean the maximum of f over E i.e. max j∈E f (j). The Kronecker symbol δ(i, j) will be used to denote the identity matrix indexed by E: ∀i, j ∈ E δ(i, j) = 0 if i = j, δ(i, j) = 1 if i = j.

Let m be a positive integer. The state space of most of Markov chains under consideration will be E m , the set of m-uples of elements of E: these m-uples are called populations and their components individuals. They will be denoted by the letters x, y, z.

If x = (x 1 , • • • , x m ) belongs to E m , [x] is the set of individuals contained in x i.e. [x] = { x k : 1 ≤ k ≤ m }
and for i in E, x(i) is the number of occurrences of i in the population x:

x(i) = card { k : 1 ≤ k ≤ m, x k = i }.
The group of permutations S m of the set {1, • • • , m} operates on E m in the following way:

∀σ ∈ S m σ • x = σ • (x 1 , • • • , x m ) = (x σ(1) , • • • , x σ(m) ).
The set E m /S m is the set of equivalence classes associated with this group operation. We associate with f two functions defined on E m . The first one is still noted f and takes its values in (R * + ) m :

f (x) = f (x 1 , • • • , x m ) = (f (x 1 ), • • • , f (x m )).
The second one is noted f and is real-valued:

f (x) = f (x 1 , • • • , x m ) = max 1≤k≤m f (x k ).
For x in E m , x denotes the set of those elements of [x] which realize the value f (x):

x = { x k : 1 ≤ k ≤ m, f (x k ) = f (x) }.
For a point i of E, (i) is the m-uple whose m components are equal to i and A is the set of all such m-uples (which are called the uniform populations). By S we denote the set of equi-fitness populations, that is, populations whose individuals have the same fitness:

S = { x ∈ E m : f (x 1 ) = • • • = f (x m ) }.
We sometimes identify f * with { (i) : i ∈ f * } so that f * may be seen as a subset of A.

III. Description of the model

1. The unperturbed Markov chain (X ∞ n ) In the absence of perturbations, the process under study is a Markov chain (X ∞ n ) n≥0 with state space E m . The superscript ∞ reflects the fact that this process describes the limit behavior of our model, when all perturbations have disappeared. The transition probabilities of this chain are

P X ∞ n+1 = z/X ∞ n = y = 1 (card y) m m k=1 1 b y (z k )y(z k ) = 1 (card y) m i∈[z] 1 b y (i)y(i) z(i)
that is, the individuals of the population X ∞ n+1 are chosen randomly (under the uniform distribution) and independently among the elements of X ∞ n which are the best individuals of X ∞ n according to the fitness function f . Suppose the chain starts up with the initial population X 0 = x 0 . Then ∀n ≥ 1 [X ∞ n ] ⊂ x 0 and with probability one, after a finite number of steps N , the chain is absorbed in a state (i) where i belongs to x 0 . In particular, if x 0 is reduced to one point i, the chain is instantaneously absorbed in (i).

The perturbed Markov chain (X l

n ) The previous Markov chain (X ∞ n ) is randomly perturbed by three distinct mechanisms. The first two act directly upon the population and mimic the phenomena of mutation and crossover. The third one consists in loosening the selection of the individuals. The intensity of the perturbations is governed by an integer parameter l. As l grows toward infinity, the perturbations progressively disappear. The perturbed Markov chain (X l n ) is obtained through overlapping of several other chains (U l n ), (V l n ) which represent the successive populations obtained by applying the perturbing operators. More precisely, we decompose the transition from X l n to X l n+1 in three stages:

X l n mutation ------→ U l n crossover ------→ V l n selection ------→ X l n+1 .
We now proceed to a more detailed description of these three operations.

X l n -→ U l n : mutation. The mutation operator is modeled by random independent perturbations of the individuals of the population X l n . Such a perturbation is described by a markovian kernel p l on the space E, that is a function defined on

E × E with values in [0, 1] verifying ∀i ∈ E j∈E p l (i, j) = 1.
The transition probabilities from X l n to U l n are then given by ( 1)

P U l n = u/X l n = x = p l (x 1 , u 1 ) • • • p l (x m , u m ).
This perturbation is small whenever the matrix (p l (i, j)) (i,j)∈E×E is close to the identity matrix. To ensure the vanishing of mutations when l grows toward infinity, we will impose

(2) ∀i, j ∈ E lim l→∞ p l (i, j) = δ(i, j).
U l n -→ V l n : crossover. The crossover operator is modeled by random independent perturbations of the couples formed by consecutive individuals of the population (X l n ). As in the mutation case, such a perturbation is described by a markovian kernel q l on the space E × E, that is a function defined on (E × E) × (E × E) with values in [0, 1] verifying:

∀(i 1 , j 1 ) ∈ E × E (i 2 ,j 2 )∈E×E q l ((i 1 , j 1 ), (i 2 , j 2 )) = 1.
The transition probabilities from U l n to V l n are

(3)

P V l n = v/U l n = u = δ m (u m , v m ) 1≤k≤m/2 q l ((u 2k-1 , u 2k ), (v 2k-1 , v 2k ))
where δ m (i, j) = δ(i, j) if m is odd (the last individual of the population remains unchanged after crossover) and δ m (i, j) = 1 if m is even.

To ensure the vanishing of crossovers when l grows toward infinity, we will impose

(4) ∀(i 1 , j 1 ) ∈ E × E ∀(i 2 , j 2 ) ∈ E × E lim l→∞ q l ((i 1 , j 1 ), (i 2 , j 2 )) = δ(i 1 , i 2 )δ(j 1 , j 2 ).
V l n -→ X l n+1 : selection. In order to build our selection operator, we will use a selection function.

Definition. A selection function of order m is a function F defined on {1, • • • , m}×(R * + ) m with values in [0, 1] satisfying for each (f 1 , • • • , f m ) in (R * + ) m : a) m k=1 F (k, f 1 , • • • , f m ) = 1, b) ∀σ ∈ S m ∀k ∈ {1, • • • , m} F (σ(k), f σ(1) , • • • , f σ(m) ) = F (k, f 1 , • • • , f m ), c) f 1 ≥ f 2 ≥ • • • ≥ f m ⇒ F (1, f 1 , • • • , f m ) ≥ F (2, f 1 , • • • , f m ) ≥ • • • ≥ F (m, f 1 , • • • , f m ).
The value

F (k, f 1 , • • • , f m ) is the probability of choosing f k among the values f 1 , • • • , f m :
this probability does not depend on the indices and increases with the relative value of

f k among f 1 , • • • , f m . Now let F l be a selection function. The m individuals X l,1 n+1 , • • • , X l,m n+1
who make up the population X l n+1 are chosen randomly and independently in the population V l n from the law defined by F l :

∀ r ∈ {1, • • • , m} ∀i ∈ E P X l,r n+1 = i = h:V l,h n =i F l h, f (V l n ) . (recall that f (V l n ) is the m-uple f (V l,1 n ), • • • , f (V l,m n )
) so that the transition probabilities from V l n to X l n+1 are

(5)

P X l n+1 = x/V l n = v = i∈[x] k:v k =i F l (k, f (v)) x(i) = m r=1 k:v k =x r F l (k, f (v)).
The selection pressure is maximal if the individuals of X l n+1 are chosen randomly and uniformly from the fittest individuals of V l n . The unique selection function F ∞ which implements such a selection scheme is defined by

F ∞ (k, f (x)) = 1 b x (x k ) card x
i.e. we have then the uniform distribution over x.

To ensure the disappearance of the selection of individuals below peak fitness, we will impose the convergence of F l toward F ∞ on the set f (E) m :

(6) ∀x ∈ E m ∀k ∈ {1, • • • , m} lim l→∞ F l (k, f (x)) = F ∞ (k, f (x)).
The basic example of such a sequence of selection functions is given by ( 7)

F l (k, f 1 , • • • , f m ) = exp(f k ln w l )
m r=1 exp(f r ln w l ) where (w l ) l≥0 is an increasing sequence of positive real numbers which tends to infinity. We will only be concerned with such sequences of selection functions in the sequel. We can now evaluate the Transition probabilities of the chain (X l n ). The transition probability P X l n+1 = z/X l n = y is given by the sum

(u,v)∈(E m ) 2 P X l n+1 = z/V l n = v P V l n = v/U l n = u P U l n = u/X l n = y .
Conditions (2), ( 4) and ( 6) imply

∀(y, z) ∈ E m × E m lim l→∞ P X l n+1 = z/X l n = y = P X ∞ n+1 = z/X ∞ n = y
so that the transition probabilities of (X l n ) converge toward those of (X ∞ n ) as l tends to infinity. Thus the Markov chain (X l n ) appears as a perturbation of the Markov chain (X ∞ n ). In order to deal with the asymptotic dynamics of the chain (X l n ) we need more information about the way its transition probabilities converge to those of (X ∞ n ). We will make the following assumptions on the sequences of kernels p l and q l . Hypothesis H p . There exist a sequence (u l ) l≥0 of positive real numbers and a function α defined on E × E with values in R + satisfying a) lim

l→∞ u l = +∞ b) α is an irreducible kernel, that is ∀i, j ∈ E ∃i 0 , i 1 , • • • , i r i 0 = i, i r = j, 0≤k≤r-1 α(i k , i k+1 ) > 0 c) p l admits the development ∀i, j ∈ E ∀s > 0 p l (i, j) = α(i, j)(u l ) -1 + o ((u l ) -s ) if i = j 1 -α(i, j)(u l ) -1 + o ((u l ) -s ) if i = j
Hypothesis H q . There exist a sequence (v l ) l≥0 of positive real numbers and a function β defined on (E × E) × (E × E) with values in R + satisfying a) lim l→∞ v l = +∞ b) q l admits the development

∀(i 1 , j 1 ) ∈ E × E ∀(i 2 , j 2 ) ∈ E × E ∀s > 0 q l ((i 1 , j 1 ), (i 2 , j 2 )) = β ((i 1 , j 1 ), (i 2 , j 2 )) (v l ) -1 + o ((v l ) -s ) if (i 1 , j 1 ) = (i 2 , j 2 ) 1 -β ((i 1 , j 1 ), (i 2 , j 2 )) (v l ) -1 + o ((v l ) -s ) if (i 1 , j 1 ) = (i 2 , j 2 )
In fact, the asymptotic dynamics of the chain will be governed by the kernels α and β: the irreducibility condition on the kernel α (which implies the irreducibility of the kernel p l for l large enough) is essential in order to allow the population to visit all the space E even when the perturbations are small.

IV. Convergence of the homogeneous algorithm

This section is devoted to the study of the behavior of the chain (X l n ) when n first goes to infinity and then l does. We will be interested in the quantities lim

l→∞ lim n→∞ P X l n = z/X l 0 = y .
The kernel α is irreducible so that for l large enough, the chain (X l n ) is irreducible. Furthermore, for l sufficiently large, the diagonal coefficients P X l n+1 = x/X l n = x are positive so that the chain (X l n ) is aperiodic. In addition, the state space of (X l n ) is finite. The chain (X l n ) admits thus a unique invariant probability measure µ l [11, Theorems 4.1 and 4.2]. This measure charges all points of the space E m and we have

∀ y, z ∈ E m lim n→∞ P X l n = z/X l 0 = y = µ l (z).
We are thus interested in the convergence of the stationary measures (µ l ) l≥0 . We will note P x for the probability measure associated with the chain (X l n ) starting at X l 0 = x and E x will be the expectation with respect to P x .

1. The time of entrance in the set S of X l n Let τ l = min { n > 0 : X l n ∈ S } be the first entrance time of the chain (X l n ) into S and

s l = min x∈E m k:x k ∈b x F l (k, f (x)).
For each n in N and x in E m we have

P X l n+1 ∈ S/X l n = x ≥ min v∈E m P [X l n+1 ] ⊂ v/V l n = v ≥ (s l ) m .
From this we deduce that

(8) ∀q ∈ N ∀x ∈ E m P x τ l > q ≤ (1 -(s l ) m ) q .
Summing up this inequality from q = 0 to ∞ we get

∀x ∈ E m E x τ l ≤ (s l ) -m .
We proceed now to a more careful study of the sequence (s l ) l≥0 . With our choice for the sequence F l , it is clear that s l is strictly positive for every l. In addition, for a fixed population x, F l (k, f (x)) increases with l (recall that w l is an increasing sequence) whenever x k belongs to x. The sequence s l , being the minimum over E m of sums of such increasing sequences, is also increasing. Define

(9) δ = min{ |f (i) -f (j)| : i, j ∈ E, f (i) = f (j) }.
For x in E m , we have

1 - k:x k ∈b x F l (k, f (x)) = r:x r ∈b x exp((f (x r ) -f (x)) ln w l ) card x + r:x r ∈b x exp((f (x r ) -f (x)) ln w l ≤ (m -1)(w l ) -δ whence finally (10) 1 -s l ≤ (m -1)(w l ) -δ .

Concentration of µ l on S

We denote by (Z l n ) the Markov chain induced by (X l n ) on S and by ν l its invariant probability measure µ l S . More precisely, (Z l n ) is the Markov chain with state space S and with transition probabilities

∀x, y ∈ S P (Z l n+1 = y/Z l n = x) = P x (X τ l = y)
(we recall that τ l is the hitting time of S). We first use the standard formula for representing the invariant measure of (X l n ) with the help of the induced Markov chain (Z l n ) [7, Proposition 5.3]. We obtain [START_REF] Seneta | Non-negative matrices and Markov chains[END_REF] ∀x

∈ E m µ l (x) = µ l (S) y∈S ν l (y)E y   τ l -1 k=0 1 {X l k =x}   .
The study of the asymptotic behavior of µ l will entirely rely on this representation formula. If x does not belong to S then for y in S

E y   τ l -1 k=0 1 {X l k =x}   = E y   τ l -1 k=1 1 {X l k =x}   ≤ E y [τ l ] -1 ≤ (s l ) -m -1.
Reporting this inequality in the representation formula [START_REF] Seneta | Non-negative matrices and Markov chains[END_REF] yields

µ l (x) ≤ µ l (S) (s l ) -m -1 .
The convergence of the selection functions (F l ) l≥0 toward F ∞ implies that s l converges to one whence lim

l→∞ µ l (E m \ S) = 0. If x belongs to S then for y in S E y   τ l -1 k=0 1 {X l k =x}   = δ(x, y).
Reporting this equality in the representation formula (11) yields µ l (x) = µ l (S)ν l (x). We have just shown that µ l (S) converges to one as l goes to infinity. Therefore it remains to study the convergence of the sequence of probability measures ν l .

Asymptotics of the transition matrix of Z l n

To study the convergence of the measure ν l we will use the powerful machinery developed by Freidlin and Wentzell. The first step consists in evaluating the asymptotics of the transition probabilities of the induced Markov chain (Z l n ). Thus our next task will be to obtain estimates of the quantities P Z l n+1 = y/Z l n = x for x and y in S. We will successively study the asymptotics of the transition probabilities

P U l n = u/X l n = x , P V l n = v/U l n = u , P X l n+1 = x/V l n = v . Asymptotics of P U l n = u/X l n = x . Define for x and u in E m d(x, u) = card { k : 1 ≤ k ≤ m, x k = u k } and extend α to E m × E m by putting α(x, u) = k:x k =u k α(x k , u k ).
Define for x in E m the set U(x) as the set of populations reachable from x in one transition step through the kernel α that is U(x) = { u ∈ E m : α(x, u) > 0 }. Using the development of p l given by hypothesis H p in (1), we see that, as l → ∞, u) .

If u ∈ U(x) P U l n = u/X l n = x = o (u l ) -s ∀s > 0 If u ∈ U(x) P U l n = u/X l n = x ∼ α(x, u)(u l ) -d(x,
Asymptotics of P V l n = v/U l n = u . Define for u and v in E m d(u, v) = card { k : 1 ≤ k ≤ m/2, (u 2k-1 , u 2k ) = (v 2k-1 , v 2k ) } and extend β to E m × E m by putting β(u, v) = δ m (u m , v m ) k:(u 2k-1 ,u 2k ) =(v 2k-1 ,v 2k ) β((u 2k-1 , u 2k ), (v 2k-1 , v 2k )).
For u in E m , let V(u) be the set of populations reachable from u in one transition step through the kernel β that is

V(u) = { v ∈ E m : β(u, v) > 0 }. Using the development of q l
given by hypothesis H q in (3), we obtain, as l → ∞,

If v ∈ V(u) P V l n = v/U l n = u = o (v l ) -s ∀s > 0 If v ∈ V(u) P V l n = v/U l n = u ∼ β(u, v)(v l ) -d(u,v) .
Asymptotics of P X l n+1 = z/V l n = y . Formulas ( 5) and ( 7) yield

P X l n+1 = z/V l n = y = i∈[z] y(i) exp(f (i) ln w l ) m k=1 exp(f (y k ) ln w l ) z(i) whence P X l n+1 = z/V l n = y ∼ l→∞ 1 (card y) m i∈[z] y(i) z(i) exp   i∈[z] z(i)f (i) -m f (y) ln w l  
The important quantity in this formula is the exponent of w l which gives the rate of decreasing of the transition probabilities: it may be written

m 1 m m k=1 f (z k ) -f (y)
i.e. m multiplied by the difference between the mean of f over the population z and the maximum of f over the population y.

Choice of the sequences u l , v l and w l . We are interested in the dynamics of the chain (X l n ) when the perturbations become smaller and smaller that is when l grows to infinity. In this situation, with overwhelming probability, the chain (X l n ) behaves as the chain (X ∞ n ) would do. However the occurrence of rare events (which have very low probability) allows the chain (X l n ) to escape from the absorbing states of the chain (X ∞ n ) and to visit all the space. These rare events are caused by the three distinct mechanisms of perturbations described earlier. To ensure that these three mechanisms actually play their own role in the asymptotic dynamics of the chain, we must give the same order of intensity to the three kinds of perturbations. If one perturbation is negligible compared to another one, its asymptotic influence will be null. We will therefore focus on the situation where the three sequences u l , v l and w l which express the intensity of each perturbation have asymptotics logarithmically of the same order. In the sequel we make the choice u l = l a , v l = l b , w l = l c where a, b, c are three positive real numbers.

Let z 1 and z 2 be two arbitrary populations. As high as l might be, the probability for the chain (X l n ) to travel from z 1 to z 2 is non-zero. Let r be a positive integer. We will now study the

Asymptotics of P X l n+r = z 2 /X l n = z 1 .
This transition probability may be written as

P X l n+r = z 2 , V l n+r-1 = v r-1 , U l n+r-1 = u r-1 , X l n+r-1 = x r-1 , • • • , X l n+1 = x 1 , V l n = v 0 , U l n = u 0 , X l n = z 1
the sum being taken over all possible values of sequences of populations in

E m u 0 , • • • , u r-1 , v 0 , • • • , v r-1 and x 1 , • • • , x r-1 .
Anyway the sum is finite and to estimate it we need only to take into account the terms which prevail at infinity. Since the kernel α is irreducible, for r large enough, there is always a path leading from z 1 to z 2 which involves only transition probabilities of order l -s for some positive s. Therefore at least one term in the above sum is asymptotically of order l -t for some positive t.

If the path P:

z 1 = x 0 → u 0 → v 0 → x 1 → u 1 → • • • → x r-1 → u r-1 → v r-1 → x r = z 2
contains a transition of very low probability, the corresponding term in the sum will be negligible compared to any power of 1/l. More precisely:

If there exists k in {0, • • • , r -1} such that either u k ∈ U(x k ) or v k ∈ V(u k ), the product r-1 k=0 P X l n+1 = x k+1 /V l n = v k P V l n = v k /U l n = u k P U l n = u k /X l n = x k
is negligible with respect to l -s for any positive s.

If there exists an index k in {0, • • • , r -1} such that [x k+1 ] ⊂ [v k ],
the above term vanishes. Let D r (z 1 , z 2 ) be the set of paths

x 0 → u 0 → v 0 → x 1 → u 1 → • • • → x r-1 → u r-1 → v r-1 → x r satisfying x 0 = z 1 , x r = z 2 and ∀k ∈ {0, • • • , r -1} u k ∈ U(x k ), v k ∈ V(u k ), [x k+1 ] ⊂ [v k ]
(This set may be empty for small values of r).

It follows from the preceding discussion that

P X l n+r = z 2 /X l n = z 1 = D r (z 1 ,z 2 ) P X l n+r = z 2 /V l n+r-1 = v r-1 • • • P U l n = u 0 /X l n = z 1 + o l -s
for any positive s. In the above sum, the term associated with the path P of D r (z 1 , z 2 )

z 1 = x 0 → u 0 → v 0 → x 1 → u 1 → • • • → x r-1 → u r-1 → v r-1 → x r = z 2
is asymptotically of order l -t where t is equal to

V (P) = a r-1 k=0 d(x k , u k ) + b r-1 k=0 d(u k , v k ) + c r-1 k=0 m f (v k ) - m h=1 f (x k+1 h ) .
This quantity, which will be called the cost of the path, reflects the difficulty for the process to move along the path under consideration. We need only to consider the terms with the lowest exponents. Define the function

V r on E m × E m by V r (z 1 , z 2 ) = min P∈D r (z 1 ,z 2 ) V (P)
and let D * r (z 1 , z 2 ) be the set of those elements in D r (z 1 , z 2 ) which realize the minimum V r (z 1 , z 2 ). We define the length l(P) of a path P as the number of transitions of the chain (X l n ) it involves. With each path P we associate a constant C(P) defined by ( 12)

C(P) = 0≤k<l(P) α(x k , u k )β(u k , v k )(card v k ) -m i∈[x k+1 ] (v k (i)) x k+1 (i)
and we put

C r (z 1 , z 2 ) = P∈D * r (z 1 ,z 2 )
C(P).

In particular, D r (z 1 , z 2 ) is exactly the set of paths P of length r such that C(P) > 0.

With these notations we see that:

If D r (z 1 , z 2 ) is empty, ∀s > 0 P X l n+r = z 2 /X l n = z 1 = o(l -s ). If D r (z 1 , z 2 ) is not empty, P X l n+r = z 2 /X l n = z 1 = C r (z 1 , z 2 ) l -V r (z 1 ,z 2 ) + o l -V r (z 1 ,z 2 ) .
Remark that this formula is valid even if z 1 and z 2 differ only by a permutation (in this case V r (z 1 , z 2 ) may vanish).

Our next objective is to study the way the chain (X l n ) behaves when it escapes from the set S. We first examine how the chain may travel from z 1 to z 2 staying outside S.

Asymptotics of P X l n+r = z 2 , ∀k ∈ {1, • • • , r -1}, X l n+k ∈ S/X l n = z 1 . Let D r (z 1 , z 2 ) be the set of elements belonging to D r (z 1 , z 2 ) such that ∀k ∈ {1, • • • , r -1} x k ∈ S.
Analogously define the function

V r on E m × E m by V r (z 1 , z 2 ) = min e D r (z 1 ,z 2 ) a r-1 k=0 d(x k , u k ) + b r-1 k=0 d(u k , v k ) + c r-1 k=0 m f (v k ) - m h=1 f (x k+1 h )
and let D * r (z 1 , z 2 ) be the set of those elements in D r (z 1 , z 2 ) which realize the minimum.

Put C r (z 1 , z 2 ) = P∈ e D * r (z 1 ,z 2 )

C(P).

Using exactly the same technique as in the preceding case we get

P X l n+r = z 2 , ∀k ∈ {1, • • • , r -1}, X l n+k ∈ S/X l n = z 1 = (13) C r (z 1 , z 2 ) l -e V r (z 1 ,z 2 ) + o l -e V r (z 1 ,z 2 )
whenever D r (z 1 , z 2 ) is not empty.

The minimal costs V and V . The most probable path followed by the chain (X l n ) between two populations z 1 and z 2 will be a path of minimal cost. We consider two minimal costs, depending on whether the path is constrained outside S or not:

V (z 1 , z 2 ) = inf r∈N V r (z 1 , z 2 ), V (z 1 , z 2 ) = inf r∈N V r (z 1 , z 2 ) (We make the convention that V 0 (z 1 , z 2 ) = V 0 (z 1 , z 2 ) = ∞ if z 1 = z 2 and 0 if z 1 = z 2 ). Since D r (z 1 , z 2 ) ⊂ D r (z 1 , z 2 ), then V r (z 1 , z 2 ) ≤ V r (z 1 , z 2
) and taking the minimum over r we obtain

V (z 1 , z 2 ) ≤ V (z 1 , z 2 ).
We first state some elementary properties of the functions V r and V r .

Lemma 3.1. Let r be a positive integer.

For every integers k 1 and

k 2 such that 1 ≤ k 1 ≤ k 2 ≤ r -1, V r (z 1 , z 2 ) = min z k 1 ,z k 2 ∈ E m V k 1 -1 (z 1 , z k 1 ) + V k 2 -k 1 (z k 1 , z k 2 ) + V r-k 2 +1 (z k 2 , z 2 ) V r (z 1 , z 2 ) = min z k 1 ,z k 2 ∈E m \S V k 1 -1 (z 1 , z k 1 ) + V k 2 -k 1 (z k 1 , z k 2 ) + V r-k 2 +1 (z k 2 , z 2 ) Corollary 3.2. ∀r ∈ N * ∀z 1 , z 2 ∈ E m V r (z 1 , z 2 ) = min x 1 •••x r-1 ∈ E m V 1 (z 1 , x 1 ) + V 1 (x 1 , x 2 ) + • • • + V 1 (x r-1 , z 2 ) V r (z 1 , z 2 ) = min x 1 •••x r-1 ∈E m \S V 1 (z 1 , x 1 ) + V 1 (x 1 , x 2 ) + • • • + V 1 (x r-1 , z 2 )
Lemma 3.3. Let r be a positive integer. We have

∀z 1 ∈ E m ∀z 2 ∈ E m \ S V r (z 1 , z 2 ) ≥ V r (z 1 , z 2 ) > 0 ∀z 1 ∈ E m ∀z 2 ∈ E m V r (z 1 , z 2 ) = 0 =⇒ [z 2 ] ⊂ z 1 and z 2 ∈ S ∀z 1 ∈ A ∀z 2 ∈ E m z 1 = z 2 =⇒ V r (z 1 , z 2 ) > 0
The above results are immediate consequences of the definition of the functions V r and V r . The next lemma is of direct interest for V and V :

Lemma 3.4. Let r * = card (E m /S m ).
For every r strictly greater than r * we have

inf k≤r * V k (z 1 , z 2 ) ≤ V r (z 1 , z 2 ), inf k≤r * V k (z 1 , z 2 ) < V r (z 1 , z 2 )
Proof. Let r be an integer strictly greater than r * and consider an element of D r (z 1 , z 2 ):

u 0 , • • • , u r-1 , v 0 , • • • , v r-1 , x 0 , • • • , x r .
Since r > card (E m /S m ) then necessarily the sequence x 1 , • • • , x r contains two populations which are equivalent modulo S m :

∃ k 1 ∈ {1, • • • , r} ∃ k 2 ∈ {1, • • • , r} k 1 < k 2 ∃ σ ∈ S m σ • x k 1 = x k 2 .
In each of the three sequences u k , v k , x k we remove the elements whose index lies in

{k 1 , • • • , k 2 -1}. We obtain an element of D r-(k 2 -k 1 ) (z 1 , z 2 ) whose cost is less than the original one, whence V r-(k 2 -k 1 ) (z 1 , z 2 ) ≤ V r (z 1 , z 2 ).
Suppose now that the original sequence was in fact in the set D r (z 1 , z 2 ): necessarily x k 1 and therefore x k 2 were in E m \ S; by lemma 3.3, the cost of the path between x k 1 and x k 2 is strictly positive and the cost of the new path between z 1 and z 2 (i.e. the element of D r-(k 2 -k 1 ) (z 1 , z 2 ) obtained by removing the cycle

x k 1 → • • • → x k 2 ) is strictly smaller than the original one, whence V r-(k 2 -k 1 ) (z 1 , z 2 ) < V r (z 1 , z 2
). An immediate descending induction on r gives the desired inequalities.

As a consequence, we have the following

Corollary 3.5. Let r * = card (E m /S m ). Then V (z 1 , z 2 ) = inf k≤r * V k (z 1 , z 2 ), V (z 1 , z 2 ) = inf k≤r * V k (z 1 , z 2 ). Let I(z 1 , z 2 ) (respectively I(z 1 , z 2 )) be the set of integers k such that V (z 1 , z 2 ) = V k (z 1 , z 2 ) (respectively V (z 1 , z 2 ) = V k (z 1 , z 2 )). Let D(z 1 , z 2
) and D(z 1 , z 2 ) be the set of all paths realizing the minimum in the functions V and V :

D(z 1 , z 2 ) = k∈I(z 1 ,z 2 ) D * k (z 1 , z 2 ), D(z 1 , z 2 ) = k∈ e I(z 1 ,z 2 ) D * k (z 1 , z 2 ).
Lemma 3.4 and corollary 3.2 show that these sets are non-empty; furthermore the set D(z 1 , z 2 ) contains at least a path whose length is less than r * and the length of each path in D(z 1 , z 2 ) is less than r * (so that D(z 1 , z 2 ) is finite). Finally, we put

C(z 1 , z 2 ) = P∈D(z 1 ,z 2 ) C(P), C(z 1 , z 2 ) = P∈ e D(z 1 ,z 2 )
C(P).

We are now in position to estimate the transition probabilities of the chain (Z l n ). Let y be in E m and z in S. A straightforward application of the Markov property to the chain (Z l n ) shows that P Z l n+1 = z/Z l n = y = P y X l τ l = z . We decompose this quantity according to the possible values of τ l : (14)

P y X l τ l = z = q k=1 P y X l τ l = z, τ l = k + P y X l τ l = z, τ l > q .
However we have already shown in ( 8) and ( 10) that

P y τ l > q ≤ (1 -(s l ) m ) q ≤ m q (1 -s l ) q ≤ m q (m -1) q l -cδq .
We choose the integer q large enough to ensure that q > r * and cδq > V (y, z). Since

P y X l τ l = z, τ l = k = P y X l k = z, ∀h ∈ {1, • • • , k -1}, X l h ∈ S = C k (y, z) l -e V k (y,z) + o l -e V k (y,z)
then for each fixed integer q greater than r * ,

q k=1 P y X l τ l = z, τ l = k ∼ l→∞ C(y, z) l -e V (y,z) .
In addition, the second term in the right-hand side of ( 14) is dominated by l -cδq which is negligible compared to l -e V (y,z) whence finally (15)

P Z l n+1 = z/Z l n = y = P y X l τ l = z ∼ l→∞ C(y, z) l -e V (y,z) .
Now, to evaluate the stationary measure of the chain (Z l n ) we need

Some results from the Freidlin-Wentzell theory

In this section we restate word for word some key results from the Freidlin-Wentzell theory. This material is extracted from [START_REF] Freidlin-A | Random perturbations of dynamical systems[END_REF], chapter 6, Lemmas on Markov Chains.

Let H be a finite set and let a subset W be selected in H. A graph consisting of arrows i → j (i ∈ H \ W, j ∈ H, i = j) is called a W -graph if it satisfies the following conditions:

(1) every point i ∈ H \ W is the initial point of exactly one arrow;

(2) there are no closed cycles in the graph. We note that condition (2) can be replaced by the following condition:

(2 ) for any point i ∈ H \ W there exists a sequence of arrows leading from i to some point j ∈ W .

We denote by G(W ) the set of W -graphs; we shall use the letter g to denote graphs. If p ij (i, j ∈ H, j = i) are numbers, then (i→j)∈g p ij will be denoted by π(g). Lemma 4.1. Let us consider a Markov chain with set of states H and transition probabilities p ij and assume that every state can be reached from any other state in a finite number of steps. Then the stationary distribution of the chain is

{( i∈H Q i ) -1 Q i , i ∈ H}, where (16) Q i = g∈G{i} π(g).
Lemma 4.2. Let us be given a Markov chain on a phase space X divided into disjoint sets X i , where i runs over a finite set H. Suppose that there exist non-negative numbers p ij (j = i, i, j ∈ H) and a number a > 1 such that

(17) a -1 p ij ≤ P (x, X j ) ≤ ap ij (x ∈ X i , i = j)
for the transition probabilities of our chain. Furthermore, suppose that every set X j can be reached from any state x sooner or later (for this it is necessary and sufficient that for any j there exist a {j}-graph g such that π(g) > 0). Then

(18) a 2-2h i∈H Q i -1 Q i ≤ ν(X i ) ≤ a 2h-2 i∈H Q i -1 Q i
for any normalized invariant measure ν of our chain, where h is the number of elements in H and the Q i are defined by formula (16).

Convergence of the measure ν l

Put for y and z in S p l yz = C(y, z) l -e V (y,z) .

Since S is finite, from the estimation (15) we deduce that for any positive there exists an integer L such that

∀ l ≥ L ∀y, z ∈ S (1 + ) -1 p l yz ≤ P Z l n+1 = z/Z l n = y ≤ (1 + ) p l yz . Put for x in S Q l x = g∈G{x} π l (g)
where G{x} denotes the set of x-graphs over S and for g in G{x},

π l (g) = (y→z)∈g p l yz .
Application of lemma 4.2 to the chain (Z l n ) yields for l ≥ L

(19) (1 + ) 2-2|S| Q l x y∈S Q l y ≤ ν l (x) ≤ (1 + ) 2|S|-2 Q l x y∈S Q l y .
We study now the asymptotic behavior of Q l x . For g in G{x},

π l (g) = (y→z)∈g C(y, z) exp - (y→z)∈g
V (y, z) ln l so that the crucial quantity for computing the asymptotics of

Q l x is W (x) = min g∈G{x} (y→z)∈g V (y, z).
Let G * (x) be the set of x-graphs which realize the above minimum.

If g belongs to G * (x), π l (g) is asymptotically of order l -W (x) . If g does not belong to G * (x), π l (g) is negligible compared to l -W (x) . Finally

Q l x ∼ l→∞ g∈G * (x) (y→z)∈g C(y, z) l -W (x) . Put D(x) = g∈G * (x) (y→z)∈g C(y, z)
and let

W * = x ∈ S : W (x) = min y∈S W (y) .
By W (W * ) we mean the value min y∈S W (y). With these notations, we see that there exists an integer L such that for l ≥ L ,

(1 + ) -1 D(x) y∈W * D(y) l W (W * )-W (x) ≤ Q l x y∈S Q l y ≤ (1 + ) D(x) y∈W * D(y) l W (W * )-W (x)
whence for l ≥ max(L, L ), by ( 19),

(

+ ) 1-2|S| D(x) y∈W * D(y) l W (W * )-W (x) ≤ ν l (x) ≤ (1 + ) 2|S|-1 D(x) y∈W * D(y) l W (W * )-W (x) 1
We conclude that:

If x does not belong to W * then lim l→∞ ν l (x) = 0. If x belongs to W * then lim l→∞ ν l (x) = D(x) y∈W * D(y) = g∈G * (x) (y→z)∈g C(y, z) y∈W * g∈G * (y) (z 1 →z 2 )∈g C(z 1 , z 2 ) We denote by ν ∞ this limit measure. The complete formula for ν ∞ is ν ∞ (x) = g∈G * (x) (y→z)∈g P∈ e D(y,z) C(P) y∈W * g∈G * (y) (z 1 →z 2 )∈g P∈ e D(z 1 ,z 2 )
C(P) where C(P) is defined by [START_REF] Trouvé | Asymptotical behaviour of several interacting annealing processes[END_REF]. The measure ν ∞ is concentrated on the set W * and charges all points of W * . The crucial question now is: what is the set W * ?

The function W and the uniform populations

We derive some properties of the functional W . First, W may equivalently be defined through V instead of V : Lemma 6.1. The following equality holds for all x in E m :

(20) W (x) = min g∈G(x) (y→z)∈g V (y, z) = min g∈G(x) (y→z)∈g V (y, z).
Proof. The proof of this lemma can be found in [4, chapter 6, lemma 4.1].

The perturbed chain (X l n ) is attracted by the populations of the set S: these populations play a role analogous to the compacta K i 's (the ω-limit sets of the trajectories of the solutions of the deterministic system) in the work of Freidlin and Wentzell. They distinguish two kinds of compacta according to their stability properties. In our situation, the populations in A will play the role of the stable compacta whereas the populations in S \ A will behave as the unstable compacta. Lemma 6.2. If x belongs to S \ A then for each i in [x],

V (x, (i)) = V (x, (i)) = 0.
Proof. This result is straightforward and similar to lemma 4.2 of [START_REF] Freidlin-A | Random perturbations of dynamical systems[END_REF].

We restate now textually lemma 4.3 of [START_REF] Freidlin-A | Random perturbations of dynamical systems[END_REF]. Lemma 6.3. a) Among the x-graphs for which the minimum in W (x) is attained there is one in which from the index y, y = x of each population in S \ A an arrow y → z is issued with V (y, z) = 0 and z in A. b) For a population x in A, the value W (x) can be calculated according to (20), considering graphs on the set of populations in A.

c) If x is a population in S \ A then W (x) = min y∈A (W (y) + V (y, x)) .
Proof. The proof of this lemma is only a rewriting of the proof of lemma 4.3 of [START_REF] Freidlin-A | Random perturbations of dynamical systems[END_REF].

As a consequence, the set W * is included in the set of uniform populations A (lemmas 3.3 and 3.4 imply that for y in A and x = y, V (y, x) > 0).

Some more notations

For λ in R + we define

f λ = f -1 ({λ}), f + λ = f -1 (]λ, ∞[), f - λ = f -1 ([0, λ[).
Since the study of W * requires only consideration of uniform populations, we will often write i instead of (i). For instance, V (i, j) stands for V ((i), (j)). For any graph g over A we define

V (g) = (i→j)∈g V (i, j).
If X and Y are two subsets of S, we denote by G X (Y ) the set of Y -graphs over X ∪ Y . For instance, we have G(X) = G S (X) = G S\X (X). We put

W X (Y ) = min g∈G X (Y ) (i 1 →i 2 )∈g V (i 1 , i 2 ) = min g∈G X (Y ) V (g).
We note G * X (Y ) the set of graphs in G X (Y ) which realize this minimum. We define G X (∅) as the union of all x-graphs over X for all x in X:

G X (∅) = x∈X G X {x}.
If g is a graph over E, its restriction g |λ to the level λ is the graph

g |λ = { (i 1 , i 2 ) ∈ g : i 1 ∈ f λ }.
for any graph g in G{i}. Taking the minimum over all g in G{i} we obtain

W (i) ≥ λ∈f (E) W f λ (f + λ ∪ f - λ ) -min j:f (j) =f (i) V (i, j) so that finally (24) min i∈E\f * W (i) ≥ λ∈f (E) W f λ (f + λ ∪ f - λ ) -max i∈E\f * min j:f (j) =f (i)
V (i, j).

We build now a near-optimal graph g.

For each λ in f (E), select a graph g λ in G * f λ (f + λ ). If λ = f (f * ), g λ is an i * -graph over f * for some i * in f * (in this case, f +
λ is empty). We define the graph g as the union of the graphs g λ :

(i 1 → i 2 ) ∈ g ⇐⇒ ∃ λ ∈ f (E) (i 1 → i 2 ) ∈ g λ .
It is easy to see that g is in G{i * }. Furthermore we have by construction

∀λ ∈ f (E) g |λ = g λ whence V (g) = λ∈f (E) V (g |λ ) = λ∈f (E) W f λ (f + λ )
from which we deduce, for the point i * of f * , ( 25)

W (i * ) ≤ λ∈f (E) W f λ (f + λ ).
Putting together the inequalities ( 21), ( 24) and ( 25), we obtain

W (i * ) < min i∈E\f * W (i) where i * is a point of f * , which clearly implies W * ⊂ f * .
The condition (21) depends strongly on the structure of the optimization problem: the essential ingredients which are hidden behind the functionals W and V are the function f and the kernels α and β which determine the asymptotic dynamics of the chain (X l n ). We derive now some stronger but more explicit conditions for the inequality (21) to hold. Clearly it is sufficient to have both

(26) ∀λ ∈ f (E) \ f (f * ) W f λ (f + λ ) ≤ W f λ (f + λ ∪ f - λ ) and (27) W f * (∅) < W f * (E \ f * ) -max i∈E\f * min j:f (j) =f (i) V (i, j).
We next prove two fundamental lemmas describing the behavior of the function V with respect to the size of the population m.

Increasing of fitness

Lemma 9.1. Let R be the smallest integer such that

∀i, j ∈ E ∃ r ≤ R ∃i 0 , • • • , i r ∈ E such that i 0 = i, i r = j, ∀k ∈ {0, • • • , r -1} α(i k , i k+1 ) > 0 i.e. R
is the minimal number of transitions necessary to join two arbitrary points of E through the kernel α (since α is irreducible and E is finite, R is finite). Put

∆ • = max i∈E\f * min max 0≤k<r (f (i)-f (i k )) : i 0 = i, r ≤ R, f (i r ) > f (i), 0≤k<r α(i k , i k+1 ) > 0 ∆ * = max i,j∈f * min max 0≤k<r (f (i) -f (i k )) : i 0 = i, r ≤ R, i r = j, 0≤k<r α(i k , i k+1 ) > .
Let i be any point in E \ f * . With these notations we have

(28) ∀m ∈ N * min j∈E,f (j)>f (i) V (i, j) ≤ aR + c(R -1)∆ • .
Let i, j be two points in f * . We have also

(29) ∀m ∈ N * V (i, j) ≤ aR + c(R -1)∆ * . Corollary 9.2. Let ∆ = max (∆ • , ∆ * ). Then, (30) sup m∈N
* max i∈E\f * min j∈E,f (j)>f (i) V (i, j), max i,j∈f * V (i, j) ≤ aR + c(R -1)∆ , (31) sup m∈N
* max y∈S\(f * ) m min z∈S, b f (z)> b f (y) V (y, z) ≤ aR + c(R -1)∆ • . Remark. Putting ∆ = max { |f (i) -f (j)| : i, j ∈ E },
we have clearly ∆ = max (∆ • , ∆ * ) ≤ ∆. All the right-hand side members of the preceding inequalities are thus smaller than (a + c∆)|E|.

Proof of lemma.

Let i be a point of E \ f * . By the very definition of ∆ • , there exists a sequence i 0 , • • • , i r of points of E such that:

i 0 = i, r ≤ R, f (i r ) > f (i), 0≤k<r α(i k , i k+1 ) > 0 and max 0≤k<r (f (i) -f (i k )) ≤ ∆ • .
Among these sequences, we choose a sequence of minimal length, so that necessarily

∀k ∈ {0, • • • , r -1} f (i k ) ≤ f (i).
We put j = i r and we define

x 0 = (i 0 ), x 1 = (i, • • • , i, i 1 ), • • • , x r-1 = (i, • • • , i, i r-1 ), x r = (j), u 0 = v 0 = x 1 , • • • , u r-2 = v r-2 = x r-1 , u r-1 = v r-1 = (i, • • • , i, i r ).
The path

x 0 → u 0 → v 0 → x 1 → u 1 → • • • → v r-2 → x r-1 → u r-1 → v r-1 → x r
clearly belongs to D r ((i), (j)) and its cost is ar

+ c r-2 k=0 (f (i) -f (i k+1 )). It follows that V (i, j) ≤ V r (i, j) ≤ ar + c r-2 k=0 f (i) -f (i k+1 ) .
from which we obtain immediately

V (i, j) ≤ aR + c(R -1)∆ •
and since f (j) > f (i), inequality (28) is proved. Let i, j be two points in f * . By the very definition of ∆ * , there exists a sequence i 0 , • • • , i r leading from i to j through the kernel α such that r ≤ R and max

0≤k<r (f (i) -f (i k )) ≤ ∆ * .
We build the sequences u k , v k , x k as above and thus obtain a similar inequality for V (i, j) (with ∆ * instead of ∆ • ), which proves inequality (29).

Proof of corollary. Inequality (30) is a direct consequence of lemma 9.1. To prove inequality (31), consider an element y of S \ (f * ) m . Pick up a point i in y. Using inequality (28), we see that there exists a point j in E such that f (j) > f (i) and V (i, j) ≤ aR+c(R-1)∆ • . We build as above a path of length smaller than R which transforms by successive mutations the individual i in j and leaves unchanged all others individuals of y and obtain

V (y, (j)) ≤ aR + c(R -1)∆ • whence min j∈E,f (j)> b f (y) V (y, (j)) ≤ aR + c(R -1)∆ • .
Taking the maximum over y in S \ (f * ) m and then over m in N * yields inequality (31).

Our last condition

Theorem 11.1. Let m be an integer such that

(33) ρm > max i∈E\f * min j∈E,f (j)>f (i) V (i, j) , (34) ρm ≥ max i,j∈f * V (i, j).
For this integer m we have

lim l→∞ µ l (f * ) = lim l→∞ ν l (f * ) = 1.
Remark. Inequality (30) shows that such an integer m always exists: the quantities on the right-hand side are bounded by (a + c∆)|E| which is a constant independent of m. In particular, the conditions are fulfilled as soon as m > (aR + c(R -1)∆ )/ min(a, b/2, cδ).

Proof. We prove that for such an m the preceding sufficient conditions (26) and ( 27) are satisfied. Let λ be in f (E) \ f (f * ) and let g be a graph in

G f λ (f + λ ∪ f - λ ).
This graph may contain a finite number of transitions from

f λ to f - λ : i 1 → i 1 , • • • , i r → i r .
The first inequality (33) implies that for each i k , 1 ≤ k ≤ r, there exists a state j k in f + λ such that V (i k , j k ) < ρm and inequality (32) of lemma 10.1 implies V (i k , j k ) < V (i k , i k ). Let g be the graph obtained by replacing the r arrows i

1 → i 1 , • • • , i r → i r by i 1 → j 1 , • • • , i r → j r . The graph g is in G f λ (f + λ ) and V (g ) ≤ V (g). This construction being valid for any graph in G f λ (f + λ ∪ f - λ ), we have W f λ (f + λ ) ≤ W f λ (f + λ ∪ f - λ ) for all λ in f (E) \ f (f * ) and the first inequality (26) is proved. Now let g be a graph in G f * (E \ f * ). This graph contains a finite number of transitions from f * to E \ f * : i 1 → j 1 , • • • , i r → j r .
Let g be the graph obtained from g by replacing the arrows i 2 → j 2 , • • • , i r → j r by the (r -1) arrows i 2 → i 1 , • • • , i r → i 1 . Inequality (34) implies that V (g ) ≤ V (g). In addition the graph g is still in G f * (E \ f * ) and may be decomposed as the union of a graph of G f * (∅) and the arrow i 1 → j 1 whence

V (g) ≥ V (g ) ≥ W f * (∅) + V (i 1 , j 1 ) ≥ W f * (∅) + ρm.
Taking the minimum over all graphs in G f * (E \ f * ) yields

W f * (E \ f * ) ≥ W f * (∅) + ρm.
Since by (33)

ρm > max i∈E\f * min j∈E,f (j)>f (i) V (i, j) ≥ max i∈E\f * min j∈E,f (j) =f (i) V (i, j) then W f * (E \ f * ) > W f * (∅) + max i∈E\f * min j∈E,f (j) =f (i) V (i, j)
and the second inequality (27) is proved.

12. Existence of a critical population size m * Our fundamental result may be stated as follows.

Theorem 12.1. Fix the space E, the fitness function f , the mutation and crossover kernels α and β and the positive constants a, b, c. There exists a critical population size m * depending upon all these objects such that, when m is greater than m * , the set W * of the minima of the virtual energy W is included in the set f * of the global maxima of the fitness function f . Moreover we have

m * ≤ aR + c(R -1)∆ min a, b/2, cδ .
(We recall that δ = min{ |f (i) -f (j)| : i, j ∈ E, f (i) = f (j) } and ∆ is defined in Lemma 9.1 and Corollary 9.2.

Moreover ∆ ≤ ∆ = max { |f (i) -f (j)| : i, j ∈ E }.)
Of course, we expect that in most practical situations, the critical size m * is much smaller than the cardinality of the search space E, which is a basic requirement for a reasonable algorithm. In fact, considering the upper bound given in theorem 12.1 with b = ∞, we obtain that m * is less than R (a + c∆ )/ min(a, cδ). The quantities ∆ and δ depend only on the values of the fitness function. In a concrete situation, δ can be chosen as the level of precision that we require from the genetic algorithm to feel. The coefficient R can be thought of as the diameter of the search space, which is of course much less than its cardinality. In addition, the general upper bound given above is very rough and might be considerably enhanced (for instance by analyzing carefully the lefthand side of (28) and using the condition (33)). In fact, the critical size is related to the difficulty of escaping from a local maximum to go to a better local maximum, rather than to the cardinality of the search space. We will illustrate this on a concrete example on the space {0, 1} N d in the last section of the paper. Let us first examine more simple examples.

Consider the space E = { 0, 1, • • • , N } and let α be a markovian mutation kernel on E such that ∀i ∈ {1, • • • , N -1} α(i, i-1)+α(i, i+1)+α(i, i) = 1, α(i, i-1) > 0, α(i, i+1) > 0, (the possible mutations from i are i → i -1 and i → i + 1), and for the points 0 and N , α(0, 1) + α(0, 0) = 1, α(0, 1) > 0, α(N, N -1) + α(N, N ) = 1, α(N, N -1) > 0.

For the crossover kernel β, we choose the identity matrix over E × E. We consider the fitness function f defined by

∀k ∈ E f (k) = 2 cos k π 2 + k .
This function has local maxima at the points 4k. If we apply directly the rough upper bound given above, we obtain m * ≤ (aN + 2c(N -1))/ min(a, b/2, c). This bound is very bad, because here the diameter of the space is of the same order as the cardinality of the space itself. However, if we compute the lefthand side of (28) and use the condition (33), we get m * ≤ 3(a + c)/ min(a, b/2, c) which does not even depend on N . Taking b large and a = c, we see that m * ≤ 6.

Clearly, the algorithm can't possibly work with one single individual: when m = 1, the selection operator does nothing at all and the fitness function f does not intervene in the transition mechanism. Perhaps two individuals is enough to solve all the optimization problems, i.e., m * = 2?

The answer is no. We give now examples where the critical size m * is arbitrarily large. We consider the space E defined in the previous example, with the same crossover and mutation kernels. Let f be a function defined on E, with values in R * + such that

f (1) = f (2) = • • • = f (N -1), f (0) = f (1) + γ, f (N ) = f (N -1) + Γ,
where 0 < γ < Γ. We have

∀i ∈ {1, • • • , N -1} V (i, i -1) = V (i, i + 1) = a,
(these transitions need only the mutation of one individual: the remainder of the population, driven by the process (X ∞ n ), then comes on the new point)

V (0, 1) = min am, a + cmγ , V (N, N -1) = min am, a + cmΓ .

The transition from 0 to 1, for instance, takes place either trough a general mutation, of cost am, or by a general anti-selection after an individual has mutated to 1, of cost a+cmγ. Finally V (0, N ) = min am + a(N -1), a + cmγ + a(N -1), aN + c(N -1)γ .

There are two possible trajectories to go from 0 to N ; either the whole population go from 0 to 1 (cost V (0, 1)) and then from 1 to N (cost a(N -1)), or an explorer starts alone from 0 and mutates until the point N , while the remainder of the population waits in 0 (cost aN + c(N -1)γ): when the explorer reaches the point N , the whole population throw themselves on N , since γ < Γ. Similarly

V (N, 0) = min am + a(N -1), a + cmΓ + a(N -1), aN + c(N -1)Γ + cm(Γ -γ) .
Notice yet a crucial difference for the passage from N to 0: for the second kind of trajectories, whenever the explorer reaches 0, it is necessary to perform a general anti-selection of cost cm(Γ -γ) to bring everyone in 0. In particular, we check that

lim m→∞ V (0, N ) = aN + c(N -1)γ < ∞ and lim m→∞ V (N, 0) m = min a, c(Γ -γ) > 0.
We evaluate now the virtual energy W .

For i in {1, • • • , N -1}, the i-graph of minimal cost is 0 → 1 • • • → i -1 → i ← i + 1 ← • • • ← N -1 ← N so that W (i) = a(N -2) + V (0, 1) + V (N, N -1).
For the point 0, there are two possible 0-graphs:

N → N -1 → • • • → 1 → 0 et N -1 → • • • → 1 → 0 ← N whence W (0) = a(N -1) + min V (N, N -1), V (N, 0) = a(N -1) + min am, a + cmΓ, aN + c(N -1)Γ + cm(Γ -γ) .
Similarly, we have W (N ) = a(N -1) + min am, a + cmγ, aN + c(N -1)γ .

Since N > 1 and V (0, 1) > a, V (N, N -1) > a, for all i in {1, • • • , N -1}, we have W (0) < W (i), W (N ) < W (i), so that W * ⊂ {0, N }. Suppose am ≤ a + cmγ et am ≤ aN + c(N -1)γ.

For these inequalities to hold, it is enough that a ≤ cγ and m ≤ N . In this situation W (0) = W (N ) = a(N -1) + am and the set W * = {0, N } is not included in f * = {N }. Thus, whenever a ≤ cγ, we have m * > N .

Let us try to explain this phenomenon: the exits from the attraction basins of 0 and N 29 take place with a general mutation, of cost am, which brings the whole population in 1 and N -1. The heights γ and Γ do not intervene and the dynamics can't discriminate between the points 0 and N . The limiting law is then concentrated on both points {0, N }. Whenever the mutation kernel α is not symmetric, we may in addition build examples where f * ∩ W * = ∅ whenever m < m * . Consider the space

E = { 0, 1, • • • , N, N + 1, • • • , N + M -1 } and let f : E → R * + be such that f (1) = f (2) = • • • = f (N -1) = f (N + 1) = • • • = f (N + M -1), f (0) = f (1) + γ, f (N ) = f (N -1) + Γ,
where 0 < γ < Γ. Let α be a markovian mutation kernel satisfying

∀i ∈ {0, • • • , N + M -2} α(i, i + 1) = 1 and α(N + M -1, 0) = 1. We have V (N, 0) = min am + a(M -1), a + cmΓ + a(M -1), aM + c(M -1)Γ + cm(Γ -γ) , V (0, N ) = min am + a(N -1), a + cmγ + a(N -1), aN + c(N -1)γ , and
W (0) = a(N + M -2) + min am, a + cmΓ, aM + c(M -1)Γ + cm(Γ -γ) , W (N ) = a(N + M -2) + min am, a + cmγ, aN + c(N -1)γ .
To have W (0) < W (N ) (and

W * = {0}), it is enough that aM + c(M -1)Γ + cm(Γ -γ) < min am, a + cmγ, aN + c(N -1)γ). Equivalently aM + c(M -1)Γ < m(a -c(Γ -γ)), aM + c(M -1)Γ < a + cm(2γ -Γ), aM + c(M -1)Γ + cm(Γ -γ) < aN + c(N -1)γ.
We choose γ, Γ such that 0 < Γ/2 < γ < Γ and a, c such that a > c(Γ -γ).

Fix an integer M . The first two inequalities are fulfilled for m sufficiently large. Fix such a value of m. The parameter N appears only in the third inequality. All other parameters being fixed, this inequality is satisfied for N sufficiently large.

We have thus built a space E, a mutation kernel α, a function f : E → R * + and a set of parameters (a, c, m) (where the population size m may be arbitrarily large) such that W * ∩ f * = ∅.

Both hypothesis H p and H q are clearly fulfilled. We put as before

∆ • = max i∈E\f * min max 0≤k<r (f (i)-f (i k )) : i 0 = i, r ≤ R, f (i r ) > f (i), 0≤k<r H(i k , i k+1 ) > 0 ∆ * = max i,j∈f * min max 0≤k<r (f (i) -f (i k )) : i 0 = i, r ≤ R, i r = j, 0≤k<r H(i k , i k+1 ) > 0 and let ∆ = max (∆ • , ∆ * ).
Within this framework, all convergence results proved in part IV do apply. We restate and comment only the rough results, which take a more precise form in this context.

The homogeneous case

Theorem 2.1. (Convergence of the homogeneous simple genetic algorithm) If

(35) m > aN + c(N -1)∆ min (a, b/2, cδ) then ∀x ∈ E m lim l→∞ lim n→∞ P [X l n ] ⊂ f * /X l 0 = x = 1.
Our condition shows that the optimization problem may be solved with a sufficiently large population size. In addition, we are completely free for the choice of the parameters a, b and c: if we take a very large c, so that min (a, b/2) ≤ cδ, we see that δ does not intervene any more in the above condition. The parameter ∆ describes the difficulty for an individual to travel from one point to a better point. If ∆ is fixed and N becomes large, a population size of order CN for some constant C will always suffice to handle the problem. Finally, let us remark that our condition on the population size m is very rough: neither the fine structure of the optimization problem nor the possibility of using crossovers to travel between two populations have been taken into account to obtain it. The crucial point is that the perturbation mechanism allows the process to visit all the space, even when the random perturbations are very small. The role of the crossover is thus not fundamental (the algorithm without crossover corresponds to the case b = ∞): however this operator creates a lot of new possible transitions for the chain (X n ) and thus certainly decreases the values of the functions V and V as well as the optimal population size.

then ∀x ∈ E m lim n→∞ P ([X n ] ⊂ f * /X 0 = x) = 1.
The proofs of these results involve technical large deviations estimates. To avoid lengthy developments, we do not reproduce them here. To ensure the convergence toward f * in the inhomogeneous case, we must first ensure that the homogeneous algorithm converges (condition (35)) and that the rate of increasing of the sequence l(n) is carefully adapted to avoid the process to be trapped in any sub-optimal population.

The next very important issue is the speed of convergence of the algorithm. This question is addressed in [START_REF] Cerf | The dynamics of mutation-selection algorithms with large population sizes[END_REF] using the very technical tools developed by Catoni and Trouvé for analyzing generalized simulated annealing processes.

The ground state of the Ising model

We consider a finite box Λ = {1 • • • N } d of side length N in the d-dimensional integer lattice Z d . A point of this box is called a site. Sites will be denoted by the letters x, y. We wrap this box into a torus and we define a neighbourhood relation on Λ by: x ∼ y if all the coordinates of x and y are equal except one which differs by 1 or N -1. At each site x of Λ there is a spin taking the values -1 or +1. The set of all possible spins configurations is X = {-1, +1} Λ . Configurations of spins will be denoted by the letters η, σ. The value of the spin at site x for a configuration σ is denoted by σ(x).

The energy of a configuration σ is where h > 0 is the external magnetic field. We consider the situation where N is large and h is small, with N h > 2d. A ground state of the Ising model is a configuration realizing the global minimum of the energy. Because h is small, in order to minimize the energy, we should first try to minimize the first term appearing in the energy. The effect of this term is to make neighbouring sites have the same sign (this is the ferromagnetic interaction). The effect of the second term is to make the spins choose the sign of the magnetic field. Hence the unique ground state is the configuration with all spins up. There is another very stable configuration, the one with all spins down. We use the mutation kernel of the simple genetic algorithm defined in paragraph V.1. For any configurations σ, η, we set Hence the mutation consists in making one spin-flip at any site of the configuration. The energy landscape associated with this energy and this mutation kernel is very complicated and full of local minima. For instance any configuration which is the union of parallelepipedic droplets of + in a sea of -such that any two droplets are at a distance larger than two is a local minimum. Some features of this energy landscape have been analyzed in dimension two and three in order to deal with metastability questions associated to the Metropolis dynamics [START_REF] Ben Arous-R | Metastability of the three dimensional Ising model on a torus at very low temperatures[END_REF][START_REF] Neves | Critical droplets and metastability for a Glauber dynamics at very low temperatures[END_REF]. Partial results are also available in arbitrary dimension [START_REF] Neves | A discrete variational problem related to Ising droplets at low temperatures[END_REF]. We will rely on these results to give upper bounds on the critical population size m * associated to the problem of minimizing the energy E. The critical population size of the simple genetic algorithm running on the space X = {-1, +1} N d , with the mutation kernel defined in (37), with no crossover, and with the fitness function -E defined in (36) satisfies

m * ≤ N d a + c C(d) h d-1 (min(a, ch)) -1
where C(d) is a constant depending only on the dimension. Thus the critical size is less than a polynomial in N , while the size of the state space is 2 N d . In fact we have used here the rough upper bound given in (28). A more careful analysis would yield

m * ≤ 2d h d a + c C(d) h d-1 (min(a, ch)) -1
which is even independent of N ! Notice however that the length N has to satisfy N h > 2d. This phenomenon illustrates the fact that the critical size is related to a fitness barrier the population has to overcome to find the true global minimum, rather than to the cardinality of the search space. Here the fitness barrier corresponds to the energy necessary to escape from the metastable state (all spins down) to reach the ground state (all spins up). This energy barrier is the energy of a critical droplet (see [START_REF] Ben Arous-R | Metastability of the three dimensional Ising model on a torus at very low temperatures[END_REF][START_REF] Neves | A discrete variational problem related to Ising droplets at low temperatures[END_REF][START_REF] Neves | Critical droplets and metastability for a Glauber dynamics at very low temperatures[END_REF]).

if η = σ x for some site x 1 -

 1 N d l -a if η = σ 0 otherwise 34where σ x is the configuration σ flipped at x i.e.σ x (y) = +σ(y) if y = x -σ(y) if y = x

is satisfied, then the set W * of all points which may be settled ultimately by the homogeneous algorithm with positive probability is included in the set f * of the global maxima of the fitness function i.e.

∀x ∈ E m lim l→∞

Proof. Let g be a graph over E. We may decompose the sum V (g) in the following way:

Suppose now that g is in G{i} for some i in

Let j be any point outside f θ ; then

and the inequality being valid for all j outside f θ , we have

Summing up inequalities ( 22) and (23) yields

10. Decreasing of fitness Lemma 10.1. Let ρ = min(a, b/2, cδ) where δ is defined by [START_REF] Neves | Critical droplets and metastability for a Glauber dynamics at very low temperatures[END_REF]. Then

Remark. The above inequality, being true for the function V , is a fortiori true for the function V .

Proof. We show by induction on r that V r satisfies this inequality. Consider first the case r = 0. Either z 1 = z 2 and the right-hand member vanishes or

). We consider two situations:

) then all elements of the set

} have been destroyed during the transition between z 1 and v 0 whence

The induction hypothesis at rank r implies

and the desired inequality is true in both cases for the path under consideration. Taking the minimum over all elements of D r+1 (z 1 , z 2 ) yields the result at rank r + 1 and the induction is completed.

V. Application to the classical simple genetic algorithm

We specialize our general model in order to apply our results to

The simple genetic algorithm

We take E = {0, 1} N for some integer N . A point i of E is a word of length N over the alphabet {0, 1} and is noted i = i 1 • • • i N where i k ∈ {0, 1}. The Hamming distance H(i, j) between two points i, j of E is the number of letters where i and j differ:

The mutation kernel p l is defined by

The associated kernel α is

It is irreducible: the minimal number of transitions necessary to join two arbitrary points of E through the kernel α is R = N . In order to build the crossover operator, we define now a cutting operator

The kernel β is then defined to be

and the kernel q l is q l ((i, j), (i , j )) = β((i, j), (i , j )) l -b if (i, j) = (i , j ) and q l ((i, j), (i, j)) = 1 -

β((i, j), (i , j )) l -b .

The inhomogeneous case

In practice, one will not wait for the Markov chain to reach equilibrium before decreasing the perturbations. The idea is then to decrease the intensity of the perturbations as time goes, and one wishes do to this in the most efficient way to obtain the same limiting law. From now onwards we assume that l and n increase simultaneously. More precisely, we take l as an increasing function of n with lim n→∞ l(n) = ∞. The Markov chains of our model become thus time inhomogeneous and the transition probabilities depend on the time via the function l(n). We suppress the index l in our notations; for instance X n stands for X l n . Recall that the convergence exponent [10, chapter 3] of the increasing sequence l(n) is defined as the unique non-negative real number λ having the following property: the series

converges for θ > λ and diverges for θ < λ. We denote by T 1 , • • • , T n , • • • the instants of the successive visits of the chain (X n ) in the set S of equi-fitness populations i.e. T n = inf{ k : k > T n-1 , X k ∈ S }. We are mostly interested in the behavior of the chain (X T n ).

Theorem 3.1. (Convergence of the inhomogeneous simple genetic algorithm) 1) For the chain (X T n ) to be trapped in f * after a finite number of transitions, i.e., to have

the convergence exponent of the sequence l(n) must be a positive real number; that is, there must exist two positive real numbers θ 1 and θ 2 such that

2) If the convergence exponent λ of the sequence l(n) and the population size m satisfy the inequalities aN + c(N -1)∆ < λ < min (a, b/2, cδ) m then, with probability one, the chain (X T n ) is trapped in f * after a finite number of transitions, i.e.,

∀x ∈ E m P (∃N ∀n ≥ N [X T n ] ⊂ f * /X 0 = x) = 1.

3) Suppose there exists a real number t strictly greater than one such that for all r in N, the sequences l( tn + r) and l(n) are logarithmically equivalent. If the convergence exponent λ of the sequence l(n) and the population size m both satisfy the inequalities aN + c(N -1)∆ < min (a, b/2, cδ) m ≤ λ


Amanda Cigier

François Varray
email: francois.varray@creatis.insa-lyon.fr

Damien Garcia
email: damien.garcia@inserm.fr

SIMUS: An open-source simulator for medical ultrasound imaging. Part II: Comparison with four simulators

Background and Objective: Computational ultrasound imaging has become a well-established methodology in the ultrasound community. In the accompanying paper (part I), we described a new ultrasound simulator (SIMUS) for MATLAB, which belongs to the Matlab UltraSound Toolbox (MUST). SIMUS can generate pressure fields and radiofrequency RF signals for simulations in medical ultrasound imaging. It works in a harmonic domain and uses far-field and paraxial linear equations. Methods: In this article (part II), we illustrate how SIMUS compares with other ultrasound simulators (Field II, k-Wave, FOCUS, and Verasonics) for a homogeneous medium. We designed different transmit sequences (focused, planar, and diverging wavefronts) and calculated the corresponding 2-D and 3-D (with elevation focusing) RMS pressure fields. Results: SIMUS produced pressure fields similar to those of Field II, FOCUS, and k-Wave. The acoustic fields provided by the Verasonics simulator were significantly different from those of SIMUS and k-Wave, although the overall appearance remained consistent. Conclusion: Our simulations tend to demonstrate that SIMUS is reliable and can be used for realistic medical ultrasound simulations.

Introduction

Computational ultrasound imaging [START_REF] Verweij | Simulation of ultrasound fields[END_REF] has become popular since 20 0 0 after the availability of the Field II [START_REF] Jensen | FIELD: a program for simulating ultrasound systems[END_REF] ultrasound simulator ( Fig. 1 of Part I, [START_REF] Garcia | SIMUS: an open-source simulator for medical ultrasound imaging. Part I: theory & examples[END_REF] ). The recent growing interest in computational ultrasound imaging is probably related to the accessibility of raw signals, driven by the emergence of ultrasound open platforms [START_REF] Tortoli | ULA-OP: an advanced open platform for ultrasound research[END_REF][START_REF] Jensen | SARUS: a synthetic aperture real-time ultrasound system[END_REF][START_REF] Boni | Ultrasound open platforms for next-generation imaging technique development[END_REF][START_REF] Garcia | Make the most of MUST, an open-source Matlab UltraSound Toolbox[END_REF] dedicated to research (e.g. SARUS, UlaOP, Ultrasonix, Verasonics). The interplay between simulations and experimentations has fostered the development and validation of new techniques, both in transmit (e.g. coding, sparse arrays, transmission design) and in post-processing of ultrasound signals (e.g. beamforming, filtering, motion estimation). Access to realistic ultrasound simulations also enables teams not equipped with ultrasound scanners to participate and innovate in the processing of ultrasound signals. The potential offered by ultrasound simulations has motivated some investigators to address configurations not covered by Field II. Other simulators were then created, some of which have been mentioned in Part I. Among the open-source simulators, k-Wave quickly gained popularity and became another reference tool in computational ultrasound imaging. Based on the number of annual citations ( Fig. 1 in Part I, [START_REF] Garcia | SIMUS: an open-source simulator for medical ultrasound imaging. Part I: theory & examples[END_REF] ), k-Wave is now as widely adopted as Field II by the ultrasound imaging community. Although a substantial number of software packages have been developed to broaden the spectrum of simulated ultrasound conditions (a non-exhaustive list is given on the k-Wave website 1 ), Field II and k-Wave are the most widely employed.

In the accompanying paper (part I), we introduced a new ultrasound simulator (SIMUS) that works in the harmonic domain and is written in MATLAB language [START_REF] Garcia | SIMUS: an open-source simulator for medical ultrasound imaging. Part I: theory & examples[END_REF] . SIMUS is part of the opensource MUST Matlab UltraSound Toolbox 2 that proposes codes for ultrasound experimentations and simulations [START_REF] Garcia | Make the most of MUST, an open-source Matlab UltraSound Toolbox[END_REF] . The SIMUS program is based on far-field and paraxial equations and simulates radiofrequency RF signals for uniform linear or convex probes. SIMUS calls the function PFIELD that returns one-way (transmit) and twoway (transmit + receive) acoustic pressure fields. The theoretical aspect is described in detail in Part I. In particular, it is explained how the splitting of the transducer elements enables the use of farfield expressions. The RF signals can be analyzed with the opensource functions of the MUST toolbox [START_REF] Garcia | Make the most of MUST, an open-source Matlab UltraSound Toolbox[END_REF] (beamforming [START_REF] Perrot | So you think you can DAS? A viewpoint on delay-and-sum beamforming[END_REF] , wall II). The maps represent the RMS (root mean square) acoustic pressure fields. They are 5-cm high.

filtering, Doppler [START_REF] Madiena | Color and vector flow imaging in parallel ultrasound with sub-Nyquist sampling[END_REF] , speckle tracking [START_REF] Perrot | Back to basics in ultrasound velocimetry: tracking speckles by using a standard PIV algorithm[END_REF] , etc.) to produce realistic ultrasound images. Some examples are provided in the accompanying paper and on the MUST website. The objective of this part-II paper was to compare the one-way acoustic pressure fields returned by PFIELD with those simulated by Field II, k-Wave, and FOCUS [START_REF] Mcgough | Rapid calculations of time-harmonic nearfield pressures produced by rectangular pistons[END_REF] . We also included the Verasonics Research Ultrasound Simulator in the comparison protocol. Although the latter is a recent and therefore little cited [START_REF] Huang | Simultaneous noise suppression and incoherent artifact reduction in ultrafast ultrasound vascular imaging[END_REF] software, it could quickly gain popularity as the Verasonics research scanners are widely used by the ultrasound community. We consider it important that users be aware of how the Verasonics simulator differs from others. Since the theoretical background used in PFIELD/SIMUS is relatively close to that of Field II, except for the fact that SIMUS operates in the frequency domain, we also compared two-way PSFs (point-spread functions) generated by these two software packages. Before describing the comparison protocol and the results, we briefly review the theoretical and numerical specificities of the four simulators Field II, k-Wave, FOCUS, and Verasonics.

Field II, k-Wave, FOCUS, and Verasonics

Each ultrasound simulator has its own purpose and/or numerical specificity. Opting for one rather than the other depends essentially on the targeted goals, but also the ease of use and the realism of the simulations. The five ultrasound simulators that were compared (PFIELD, Field II, k-Wave, FOCUS, Verasonics) can all provide one-way acoustic pressure fields. In this study, we compared RMS (root mean square) acoustic pressures for standard emission sequences. We also compared two-way PSFs (SIMUS vs. Field II) after beamforming. We made two-way comparisons against Field II only, since Field II also uses a model based on point scatterers. The way PFIELD and SIMUS operate is explained in the part-I paper. Table 1 summarizes its features and those of the other three simulators.

Field II

J.A. Jensen introduced the MATLAB simulator Field in 1996 [START_REF] Jensen | FIELD: a program for simulating ultrasound systems[END_REF] , which then became Field II, after a previous work on the simulation and validation of acoustic pressure fields generated by medical ultrasound transducers [START_REF] Jensen | Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers[END_REF] . Field II works in the time domain and uses linear acoustics. The acoustic pressure field is calculated by dividing the transducer into small rectangular sub-elements and then summing their individual responses. For a given arbitrary point in the region of interest, Field II calculates the far-field spatial impulse response of a piston-like sub-element [START_REF] Stepanishen | Transient radiation from pistons in an infinite planar baffle[END_REF] . The acoustic pressure field at the same point is then deduced by convolving this local response with the excitation function of the transducer. The received responses are calculated from the acoustic reciprocity principle, by assuming that the point scatterers behave as monopole sources and do not interact with each other (weak scattering).

k-Wave

B.E. Treeby and B.T. Cox introduced the MATLAB toolbox k-Wave in 2010 [START_REF] Treeby | k-Wave: MATLAB toolbox for the simulation and reconstruction of photoacoustic wave fields[END_REF] , after a previous work on acoustic propagation models in heterogeneous media [START_REF] Cox | k-space propagation models for acoustically heterogeneous media: application to biomedical photoacoustics[END_REF] . The toolbox k-Wave works in the k-space (wavenumber domain). The acoustic field is calculated at the nodes of regular mesh grids using both finite differences (over time) and a pseudo-spectral [START_REF] Trefethen | Spectral methods in MATLAB[END_REF] method (in space). Point monopole sources and sensors are assigned to the grid points that belong to each transducer element. The received responses naturally derive from the numerical acoustic equations included in the software: backscattering occurs whenever an incident wave encounters an acoustic impedance gradient. In contrast to Field II and PFIELD, k-Wave can implicitly handle non-linear propagations and heterogeneous media [START_REF] Treeby | Modeling nonlinear ultrasound propagation in heterogeneous media with power law absorption using a kspace pseudospectral method[END_REF] . Note, however, that some non-linear aspects were covered with add-ons to Field II [ 19 , 20 ] by considering only the fundamental and second harmonics and using the angular spectrum approach.

FOCUS

R.J.

McGough introduced FOCUS to provide accurate pressure fields in the very near field where solving the acoustic equation can be a numerical issue. To this end, FOCUS uses the fast nearfield method developed by the author [START_REF] Chen | A 2D fast near-field method for calculating near-field pressures generated by apodized rectangular pistons[END_REF] . The fast near-field method reduces a double integral -which is present in the impulse response approach -to a rapidly converging single integral and removes a singularity. The pressure fields are calculated on a grid by the angular spectrum approach. FOCUS can consider different transducer element geometries, layered media, and nonlinearity. Although one proceeding presented a comparison with Field II on a two-way example [START_REF] Zhu | A comparison of ultrasound image simulations with FOCUS and Field II[END_REF] , the FOCUS software does not yet offer the possibility to simulate pulse-echo scenarios. The FO-CUS software is not fully open as calculations are performed via compiled codes written in C ++ and called by MATLAB functions.

Verasonics

The Verasonics Vantage ultrasound scanners are probably the most widely used research scanners for ultrasound medical imaging. Their flexibility and accessibility to raw ultrasound data make them well suited for 2-D and 3-D ultrasound research. Verasonics recently introduced its ultrasound simulator software. It is likely to become increasingly popular as it is integrated into the Vantage platform, which is why we have chosen to include this simulator in our comparison panel. To date, there is no theoretical document that describes the Verasonics simulator. According to the Verasonics website [START_REF] Vega | Ultrasound simulator software[END_REF] , the Verasonics simulator simulates backscattered echoes from a set of point targets with a given reflectivity. The contributions of each emitting element are summed by assuming linear acoustic propagation. In the Verasonics simulator, the transmit waveform, the transmit and receive apodizations, as well as the directivities of the elements are considered to simulate RF signals. The website does not specify whether the simulator applies 2-D or 3-D equations. Since the values of the element height and elevation focus did not modify the acoustic pressures, we postulated that the Verasonics simulator is 2-D.

Method: comparison of the four simulators

One-way configurations: acoustic pressure fields

We compared the five ultrasound simulators PFIELD, Field II, k-Wave, FOCUS, and Verasonics, by simulating acoustic pressure fields with four different short-pulse emission patterns: plane waves with a tilt angle of 0 °or 10 °, diverging wave, and focused wave. The properties of the transducers and transmits are summarized in Table 2 . We tested either 2-D or 3-D conditions or both, depending on the simulator characteristics. Two ultrasound transducers were simulated, a linear and a phased array. Their properties (Table II) were those of the L11-5v and P4-2v Verasonics transducers. The medium was homogeneous and non-attenuating, with a speed of sound of 1540 m/s. Note that the current versions of PFIELD and SIMUS cannot handle inhomogeneous media but can consider attenuation. The excitation pulse length was of one wavelength (see Fig. 8 in part I, [START_REF] Garcia | SIMUS: an open-source simulator for medical ultrasound imaging. Part I: theory & examples[END_REF] ). To ensure similar transmit conditions between the software packages, the pulse generated by PFIELD was used as pressure input to the other simulators. We simulated RMS (root mean square) acoustic pressure fields and compared them on a decibel dB scale. All pressure fields were normalized so that their maxima were equal to one. The pressure points were located on Cartesian grids with steps of approximately half a wavelength in the x -and z -directions, and one sixth of a wavelength in the y -direction. The transducer elements were split in both x -and y -directions for Field II and k-Wave (to obtain 2 × 17 sub-elements). The elements were not split in the xdirection in FOCUS since a near-field method is used [START_REF] Mcgough | Rapid calculations of time-harmonic nearfield pressures produced by rectangular pistons[END_REF] . It is re-called that dividing the elements in the y -direction enables a focus in elevation (see Part I for more details). On the other hand, the elements were not split in the y -direction with PFIELD since it applies continuous delays along the elevation direction (see Equation (4) in part I). To compare PFIELD with the other simulators, we calculated the absolute pressure differences in the dB scale as 20 log (| P (PFIELD) -P (other)|), after normalization. To obtain a global comparison, we pooled all simulated pressure field values from the same simulator (for 2-D and 3-D, respectively), and calculated the inter-simulator correlation coefficients. To reduce biased trends induced by the lowest values, only acoustic pressures above -20 dB were taken into account in the calculation of the correlation coefficients.

Two-way configurations: PSF

SIMUS simulates radio-frequency RF signals from the function PFIELD. We also compared SIMUS and Field II in two-way configurations by calculating two-way PSFs (point-spread functions). A two-way PSF was defined as the real envelope returned by a single scatterer after beamforming. We tested two ultrasound 3-D conditions: an unsteered plane wave with a L11-5v linear array (Table II, leftmost column) and a focused wave with a P4-2v phased array focusing at 6 cm deep (Table II, rightmost column). To investigate the two-way PSFs, one scatterer was located at (0 cm; 2 cm) and (0 cm; 6 cm) in the plane-and focused-wave configurations, respectively. The simulated RF signals were I/Q demodulated and beamformed with a delay-and-sum (DAS) [START_REF] Perrot | So you think you can DAS? A viewpoint on delay-and-sum beamforming[END_REF] . During beamforming, we used a receive f-number of 1.5 with the linear array, while a full aperture was used with the phased array. We calculated the absolute differences of the PSF in the dB scale as 20 log (|PSF(SIMUS) -PSF(Field II)|).

Results

Fig. 1 to Fig. 7 show the simulated acoustic pressure maps and the absolute differences with respect to PFIELD. We observed a good concordance between PFIELD, Field II, k-Wave, and FOCUS with the 3-D simulations ( r 2 ≥ 0.93). There was also a good concordance between PFIELD and k-Wave under 2-D conditions. The 2-D acoustic pressure fields simulated by the Verasonics ultrasound simulator were significantly different from those generated by PFIELD and k-Wave. In 3-D, when considering all simulations, the squares of the inter-simulator correlation coefficients were all greater than 0.97 (PFIELD vs. Field II vs. k-Wave, Fig. 8 ). The correlation between PFIELD and k-Wave was lesser in 2-D ( r 2 = 0.89), as side lobes were more visible with PFIELD than with k-Wave ( Fig. 4 , bottom). The 2-D correlation coefficients associated with the Verasonics simulator were the smallest ( r 2 < 0.77). The two-way PSFs returned by SIMUS and Field II are displayed in Fig. 9 and Fig. 10 . II). The maps represent the RMS (root mean square) acoustic pressure fields. They are 10-cm high.

They were very similar in both plane-wave and focused-wave configurations.

The computing times were significantly different from one simulator to another. It should be noted that the notion of calculation time is somewhat subjective. Each software has its specific applications, and proposes its own numerical tricks to reduce the computation time, generally at the expense of the accuracy. In this paper, we have focused on accuracy. In addition, each software can be run in parallel, with different algorithmic strategies. An exhaustive comparative study would require testing several test cases, each with several tuning parameters and operating systems. As an indication only (and under the request of reviewers), the computing times on the same personal laptop (non-parallel computations) were for the "plane wave" sequence ( Fig. 1 ): PFIELD: 260 s, Field II: 480 s, k-Wave: 11600 s (380 s in parallel), FOCUS: 1700 s; for the "diverging wave" sequence ( Fig. 6 ): PFIELD: 15 s, Field II: 80 s, k-Wave: 73,60 0 s (170 0s in parallel), FOCUS: 210 0 s.

Discussion

The new ultrasound simulator SIMUS generated acoustic pressure fields similar to those simulated by the other software packages Field II, k-Wave, and FOCUS, especially in 3-D ( r 2 > 0.92, Fig. 8 ). Like these popular ultrasound simulators, SIMUS can therefore be considered to produce synthetic medical ultrasound signals, under the conditions described in parts I and II (linearity, homogeneous medium, etc.). Under 2-D conditions, although the radiation patterns appeared to be identical, we observed some differences with the Verasonics simulator ( r 2 < 0.77, Fig. 8 ). Regarding the differences found with the Verasonics commercial simulator, it is difficult to give a justification because, unlike Field II, k-Wave, and FOCUS, the theoretical and numerical methods are not publicly available. The best match was with Field II. This is unsurprising because we founded SIMUS based on the same assumptions as Field II, namely: linearity, homogeneity, and far-field equations with the use of sub-elements. Like Field II, SIMUS uses weak and monopole scattering. It follows that the two-way PSFs obtained with SIMUS and Field II were similar. The main numerical difference between PFIELD/SIMUS and Field II is the computational domain: temporal with Field II, frequency with SIMUS. The frequency domain is indeed appropriate when dealing simulations with signals with bandpass signals. It can also easily consider acoustical features that depend on the frequency, such as the directivity of the elements, attenuation in the tissues, or Rayleigh scattering, for example. Another difference appears with regard to elevation focusing. In SIMUS, to simulate focusing in the elevation direction, we adjusted the velocity delays of the pistons along their height (see Fig. 4 in the accompanying part-I paper), as suggested in [START_REF] Schmerr | Fundamentals of Ultrasonic Phased Arrays[END_REF]. In Field II, the elements are split into sub-elements in the height direction. A transmission delay is associated with each of these sub-elements to provide focusing in elevation. A similar strategy is used in k-Wave.

It is important to note that the correlation coefficients reported in Fig. 8 characterize the inter-simulator similarities that we observed with a large and deep medium. The dimensions of the medium we chose are those typically encountered in B-mode or Doppler imaging. A smaller coefficient between two simulators means that they give substantially different results under these particular conditions. It does not imply that one of them is "inferior". It was difficult to compare the simulators under optimal conditions given that they each have distinct numerical characteristics. To achieve comparisons that were as fair as possible, we chose sufficiently fine grids. This strategy did not allow us to take into account the very near field, for which FOCUS and k-Wave are expected to be better adapted than Field II and PFIELD. A much finer grid (and much longer computation time) would have been necessary to examine the nearby region. In most medical imaging cases, however, it is not needed to achieve high-accuracy simulations at the near-field because the thicknesses of the acoustic matching layer and lens ( ∼1 mm for a linear transducer) [START_REF] Mckeighen | Design Guidelines For Medical Ultrasonic Arrays[END_REF] , and that of the skin ( ∼0.5 mm), separate the scanned region by several wavelengths from the surface of the transducers. As shown in the appendix of part I (Figures 14 and15), the radiation fields from PFIELD approximate well those calculated by the Rayleigh-Sommerfeld integral beyond a distance of two wavelengths (relative errors < 3%). These results confirm that the far-field equation can be used when slicing the transducer elements [START_REF] Schmerr | Fundamentals of Ultrasonic Phased Arrays[END_REF] .

When using simus and PFIELD?

SIMUS and PFIELD are part of the MATLAB Ultrasound Toolbox MUST toolbox [START_REF] Garcia | Make the most of MUST, an open-source Matlab UltraSound Toolbox[END_REF] . This toolbox, as illustrated by the numerous examples available on the website, 3 allows students and researchers to quickly learn experimental ultrasound imaging: acquisitions, demodulation, beamforming, wall filtering, color Doppler, vector Doppler, speckle tracking… SIMUS and PFIELD complete the series by offering to perform realistic ultrasound simulations for pedagogical and research purposes. In contrast to Field II, the MAT-LAB codes are open and highly documented. The default functions allow many tasks to be performed, such as those illustrated by the examples in the Part-I paper. An advanced user can modify the codes or the programming language according to her needs. PFIELD and SIMUS have been already used in several published works. For example, far-field patterns were generated with PFIELD to design diverging-wave transmit sequences for high-frame-rate motion-compensated echocardiography in human hearts [START_REF] Porée | Highframe-rate echocardiography using coherent compounding with Doppler-based motion-compensation[END_REF] . Vec-tor flow imaging with vector Doppler was simulated in a carotid bifurcation by coupling SIMUS with an SPH (smoothed particle hydrodynamics) flow model [START_REF] Shahriari | Meshfree simulations of ultrasound vector flow imaging using smoothed particle hydrodynamics[END_REF] . In our team, SIMUS has also recently been applied in the generation of simulated moving ultrasound phantoms and left ventricles to train convolutional neural networks for motion estimation in echocardiography [ 28 , 29 ]. We also simulated clinical-type Doppler-echo cineloops with the future objective to train deep learning algorithms for intracardiac flow imaging [START_REF] Sun | A pipeline for the generation of synthetic cardiac color Doppler[END_REF] . In the same vein, Milecki et al . trained a convolutional neural network by using SIMUS simulations in a synthetic murine brain microvascular network for deep-learning-based ultrasound localization microscopy [START_REF] Milecki | A deep learning framework for spatiotemporal ultrasound localization microscopy[END_REF] .

What simus cannot do, what it can

Like Field II and k-Wave, SIMUS can be used to produce realistic radiofrequency RF signals for pulse-echo ultrasound imaging. Only linear acoustics into homogeneous media, however, is considered in the current version of SIMUS. The homogeneity hypothesis means that both density and speed of sound are constant. Aberrations due to acoustic impedance gradients thus cannot be simulated at the time of publication of this article. The linearity assumption simplifies the problem since it allows the summation of the time-harmonic solutions, but this prevents one from addressing non-linear properties as in harmonic imaging. SIMUS also considers only weak (single) scattering to avoid computational overload of multiple scattering. The multiple reverberation artifacts, which are sometimes observed in clinical ultrasound, therefore cannot be simulated. In contrast, the simulator k-Wave can account for nonlinearity, acoustic heterogeneities, and multiple scattering. However, because k-Wave uses a mesh-based method, the locations of the scatterers are linked to the geometry of the grid. This represents a major difference with SIMUS (or Field II), which is a particle-based software. In SIMUS, the speckles of the B-mode images result from the pseudorandom positions of the scatterers, while they are due to local changes in impedance in k-Wave. In addition to numerical simplifications, the use of a particle-based model readily allows the motion of individual particles. This makes it an ideal simulation method for tissue displacements and the investigation of ultrasound techniques for motion analysis (e.g. speckle tracking [START_REF] Evain | A pilot study on convolutional neural networks for motion estimation from ultrasound images[END_REF] , Doppler [START_REF] Mehregan | Doppler vortography: a color Doppler approach to quantification of intraventricular blood flow vortices[END_REF] , vector flow imaging [START_REF] Shahriari | Meshfree simulations of ultrasound vector flow imaging using smoothed particle hydrodynamics[END_REF] ).

Although SIMUS is able to consider near fields through element splitting (see accompanying paper, part I, [START_REF] Garcia | SIMUS: an open-source simulator for medical ultrasound imaging. Part I: theory & examples[END_REF] ), the singularity related to the 1/r term present in the equations does not allow for evaluations near the surface of the transducer. A free software such as FOCUS removes this singularity and offers high accuracy in the very near field regions [START_REF] Chen | A 2D fast near-field method for calculating near-field pressures generated by apodized rectangular pistons[END_REF] . As explained in a previous paragraph, this aspect is not a limitation for ultrasound medical imaging for which this specific region is generally not investigated.

Potential future perspectives

In creating SIMUS, our intention was not to be "better" or more innovative than other simulators. As we mentioned earlier, each ultrasound simulator has its own specificity. Our first goal was to release an open-source MATLAB toolbox (MUST) for students and researchers in ultrasound imaging [START_REF] Garcia | Make the most of MUST, an open-source Matlab UltraSound Toolbox[END_REF] . To make it self-contained, we intended this toolbox to include a reliable, open-code ultrasound simulator that is sufficiently documented to make users aware of its capabilities and limitations. From the comparative results described in the present work, it appears that SIMUS can be used for medical ultrasound research. Since the MUST toolbox is available online, PFIELD and SIMUS have been used or tested by researchers and students on several continents ( Fig. 11 ). The present SIMUS version is for 1-D uniform (convex or linear) arrays. Depending on the success of SIMUS, we will incorporate the option to simulate matrix arrays for volume ultrasound imaging. A partition of the elements in the Y -direction of elevation will then be necessary to enable the use of the far field equations (see accompanying paper, part I). Similarly, a future version will allow the positioning and orientation of the elements at the user's convenience. Such an update would allow completing the medical imaging spectrum by considering, for example, sparse arrays or ultrasound computer tomography. It is understood that these updates will significantly increase the number of element-scatterer interactions. Some numerical strategies will have to be implemented to reduce the computation time. As we are interested in the simulation of cardiac images (echocardiography), we may also plan to integrate an add-on for the calculation of the second harmonic using a quasi-linear approximation of the nonlinear wave propagation [ 33 , 34 ]. Finally, by using a geometric ray-tracing method based on the eikonal approximation, generating non-homogeneous media would also be conceivable. This approach could be applied to derive artifacts due to aberrations, or to investigate quantitative ultrasonography methods such as the local measurement of the propagation speed.

As artificial intelligence has gained momentum in medical imaging, ultrasound simulations have found a new interest through deep learning. In medical imaging, whether ultrasound or other, the main challenge is to train neural networks with hundreds or thousands of images, which are sometimes difficult to access. Rather than being based on patient data, training can be done through patient-like simulations whose number is limited by computational time only. The key task is then to obtain simulations that are realistic enough to fool the network. "Realistic" does not mean "exact". A simulator is sufficiently "realistic" if images or signals similar to those obtained in a patient (e.g., RF or I/Q signals, Bmode or color Doppler images) can be simulated. It is understood that the subjective degree of realism depends on the intended goal. The network will then be functional when dealing with real clinical images. If simulated signals or images can be used to train neural networks with limited effort, we anticipate that ultrasound simulators still have a bright future ahead of them.

Conclusion

The open-source MATLAB simulator SIMUS is efficient and easy to use. Like Field II, k-Wave, and FOCUS, SIMUS could contribute to the growing popularization of computational ultrasound imaging.

Fig. 1 .

 1 Fig. 1. Simulations of plane waves (tilt angles = 0 °, 10 °) transmitted by a 128element linear array with PFIELD (see TableII). The maps represent the RMS (root mean square) acoustic pressure fields. They are 5-cm high.

Fig. 2 .

 2 Fig. 2. 3-D simulations: comparison between PFIELD, Field II, k-Wave, and FOCUS for the 0 °-tilt plane wave (see Fig. 1 , left). The three planes are those displayed in Fig. 1 , left. From top to bottom: azimuthal plane, elevation plane, focal plane. The green-to-red maps represent the absolute differences with respect to PFIELD.

Fig. 3 .

 3 Fig. 3. 2-D simulations: comparison between PFIELD, Verasonics, and k-Wave for the 0 °-tilt plane wave. The second row represents the absolute differences with respect to PFIELD.

Fig. 4 .

 4 Fig. 4. 3-D (top) and 2-D (bottom) simulations of a 10 °-tilt plane wave. A 3-D overview is given in Fig. 1 , right. The second rows represent the absolute differences with respect to PFIELD.

Fig. 5 .

 5 Fig. 5. Simulations of diverging and focused waves transmitted by a 64-element phased array with PFIELD (see also TableII). The maps represent the RMS (root mean square) acoustic pressure fields. They are 10-cm high.

Fig. 6 .

 6 Fig. 6. 3-D (top) and 2-D (bottom) simulations of a diverging wave. A 3-D overview is given in Fig. 5 , left. The second rows represent the absolute differences with respect to PFIELD.

Fig. 7 .

 7 Fig. 7. 3-D (top) and 2-D (bottom) simulations of a focused wave. A 3-D overview is given in Fig. 5 , right. The second rows represent the absolute differences with respect to PFIELD.

Fig. 8 .

 8 Fig. 8. Overall comparison between the five simulators after pooling all simulated acoustic pressures. The numbers are the squared correlation coefficients ( r 2 ).

Fig. 9 .

 9 Fig. 9. Comparison of two-way PSFs (SIMUS vs. Field II) in a focused-wave configuration with a simulated phased array.

Fig. 10 .

 10 Fig. 10. Comparison of two-way PSFs (SIMUS vs. Field II) in a plane-wave configuration with a simulated linear array.

Fig. 11 .

 11 Fig. 11. Cities where SIMUS has been used (2022/03).

Table 1

 1 Properties of the Five Simulators.

		PFIELD	Field II	k-Wave	FOCUS	Verasonics
	method	meshfree	meshfree	mesh	mesh	meshfree
	acoustics	linear	linear	nonlinear	nonlinear	linear
	time domain	harmonic	temporal	both	both	temporal
	space domain	x,y,z	x,y,z	x,y,z + k	x,y,z + k	x,y
	dimensionality	2D, 3D	3D	2D, 3D	3D	2D
	scattering	weak	weak	multiple	-	w e a k
	medium	homogeneous	homogeneous	heterogeneous	layered	homogeneous
	open-source	yes	no	yes	no	no
	free	yes	yes	yes	yes	no

Table 2

 2 Simulated Conditions (one-way).

		0 °-plane wave	10 °-plane wave	diverging wave	focused wave
	transducer	L11-5v (Verasonics)	P4-2v (Verasonics)
	center frequency	7.6 MHz	2.7 MHz
	elements	128	64
	kerf width	30 μm	50 μm
	element width	0.27 mm	0.25 mm
	pitch	0.3 mm	0.3 mm
	element height	5 mm (if 3-D)	1.4 cm (if 3D)
	elevation focus	1.8 cm (if 3-D)	6 cm (if 3D)
	two-way bandwidth	77%	74%
	TX focus	∞		(0, -1.64) cm	(3,5.2 cm)
	TX apodization	boxcar	Tukey	boxcar

https://www.biomecardio.com/MUST

Acknowledgment

This work was supported in part by the LABEX CeLyA ( ANR-10-LABX-0060 ) of Université de Lyon , within the program "Investissements d'Avenir" ( ANR-16-IDEX-0 0 05 ) operated by the French National Research Agency.