Thierry E Huillet

PARETO GENEALOGIES ARISING FROM A POISSON BRANCHING EVOLUTION MODEL WITH SELECTION

Keywords: Pareto coalescents, scaling limits, Poisson-Dirichlet, Kingman and Beta coalescents, Poisson Point Process, evolution model including selection

Pareto genealogies arising from a Poisson branching evolution model with selection

Introduction

We first investigate discrete-time finite coalescents derived from sampling from N i.i.d. Pareto(α) random variables, normalized by their sum, with α > 0. We include size-biasing on total length effects involving a parameter β < α. Depending on the range of α, we derive the large N limit coalescents structure: The case α ∈ [0, 1) leads to a discrete-time Poisson-Dirichlet(α, -β) Ξ-coalescent (with no time-scaling). The case α = 1 gives rise to a continuous-time beta(1, 1β) Λ-coalescent, involving a logarithmic time scaling log N . The case α ∈ (1, 2) leads to a continuous-time Beta(2α, αβ) Λ-coalescent, involving a power-law time scaling according to N α-1 . The range α ≥ 2 gives rise to the standard Kingman coalescent (with time scaling N if α > 2 and N/ log N in the critical case α = 2). We give for each case the exact speeds of convergence. We establish a loose link with a Generalized Central Limit Theorem for stable random variables and we briefly recall the main statistical features akin to general Λ-coalescents.

We indicate that the above special classes of coalescent processes (and their scaling limits) may be viewed as the genealogical processes of some forward in time evolving branching population models including selection effects. These models are closely related in spirit to the additive exponential model discussed in Brunet et al [START_REF] Brunet | Effect of selection on ancestry: an exactly soluble case and its phenomenological generalization[END_REF], [START_REF] Brunet | Genealogies in simple models of evolution[END_REF].

In the models we first consider, the size N of the population is kept constant over the generations. Each alive individual is assigned some positive fitness x > 0. In each generation and for each of the N offspring alive independently, the reproduction step is based on a fitness-dependent Poisson Point Process (PPP) with scaling power-law(α) intensity; this procedure assigns new fitnesses to the (potentially infinitely many) individuals of the next generation, in a multiplicative way. We call f (x) = x α the output fitness of an individual with fitness x. The selection step then consists of sorting out the N fittest individuals issued from the reproduction step 1 . The process is iterated independently over the subsequent generations. To the first large N approximation, the logarithm of the mean output fitness within each generation k, scaled by the generation number k, is shown to shift to the right at speed v N = log log N as k → ∞.

While adopting a sampling point of view based on the intensity of the PPP to compute the coalescence probabilities that some offspring is the one of a parental individual with given fitness, it is shown that the genealogy of the branching model with selection is in the domain of attraction of the beta(1, 1β) coalescent in the large N limit (Bolthausen-Sznitman coalescent if β = 0). It is also shown that the full class of the Pareto-coalescents discussed earlier can be obtained while considering the PPP which is the output image of the original one, given by the output map f (x) = x α in fitness space. In this setup, the large N limit computations of the coalescence and merging probabilities are based not on the fitnesses but on the output deformed fitnesses.

2.

Coalescents derived from Pareto-Sampling 2.1. Pareto sampling and coalescents. Let X 1 , .., X N be N i.i.d. Pareto(α) random variables (rvs) with P (X 1 > x) = x -α , α > 0 and x ≥ 1. Let F X1 (x) = 1 -P (X 1 > x) denote its probability distribution function (pdf). The density of X 1 is f X1 (x) = αx -(α+1) . Let S n := X n / N 1 X n , n = 1, ..., N define a random partition of the unit interval, upon normalizing the Xs by their sum. The S n s are identically distributed but not independent of course as they sum up to 1; by doing so, the unit interval [0, 1] is thus broken into N random pieces (subintervals or segments) of sizes S n , n = 1, ..., N.

By sampling the Ss, we mean that we draw independently i uniform random variables on the unit interval with i ≤ N , looking at the subintervals which are being hit in the process. From this procedure, for instance, the probability that the i-sample hits any one of the S n s only once is [START_REF] Bérard | Brunet-Derrida behavior of branching-selection particle systems on the line[END_REF] P

(N ) i,1 = E N n=1 S i n = N E S i 1 .
Let Σ N := N n=1 X n denote the partial sum of the Xs. For the values of β for which E Σ β N exists, we can size-bias the latter probability by the total length Σ N and consider instead 2

P (N ) i,1 = E Σ β N N n=1 S i n E Σ β N = N E Σ β N S i 1 E Σ β N .
The latter event under consideration corresponds to an i to 1 merger of some Markov coalescent process where i particles are identified to a single one (share the same ancestor) whenever the i-sample hits the same subinterval of the unit partition.

In this setup, P

i,1 is therefore the entry (i, 1) of its one-step transition matrix. The quantity c N := P (N ) 2,1 , which is the probability that two individuals chosen at random out of N share the same common ancestor, is called the coalescence probability.

Similarly we can define a i to j merger (j ≤ i) by considering the event that the i particles hit any size-j subset of the segments S constituting the partition of unity. We get ( 2)

P (N ) i,j = N j * i1+...+ij =i i i 1 ...i j E j l=1 S i l l = N j j l=1 (-1) j-l j l E (S 1 + ... + S l ) i ,
where the star-sum in (2) runs over the i l ≥ 1. The quantity

E j l=1 S i l l
is the probability of a (i 1 , ..., i j ) -merger from i to j. Using the same abuse of notation, we shall also write the size-biased version of the latter probability as

(3)

P (N ) i,j = N j * i1+...+ij =i i i 1 ...i j E Σ β N j l=1 S i l l E Σ β N .
Unless stated otherwise, whenever we speak in the sequel of P (N ) i,j , we mean (3) and not [START_REF] Bertoin | Lévy processes with no negative jumps and branching processes[END_REF].

Clearly,

i j=1 P (N ) i,j
= 1 and so the matrix P (N ) with entries P (N ) i,j , i = 1, ..., N and j = 1, ..., i, is a N × N lower-triangular stochastic matrix corresponding to some Markov discrete-time-k coalescent (pure death) process, say x (N ) k , with finite state-space 3 and state {1} absorbing. Let us first investigate the expression of the size-biased probability P (N ) i,1 , showing that its large N estimate depends on the understanding of the β-moments of Σ N . 2 We abusively use the same notation P (N) i,1 in the size-biased setup as in (1) (corresponding to β = 0), to avoid overburden notations.

3 A 'true' coalescent process takes values in the set of equivalence relations or partitions on {1, ..., N } and we rather deal here and throughout with its block-counting counterpart.

Proposition 1. When -∞ < β < α < 2, the size-biased probability of an i to 1 merger reads

(4) P (N ) i,1 = N α E Σ β-α N -1 E Σ β N Γ (i -α) Γ (α -β) Γ (i -β) , i ≥ 2. Proof: We have 1/S 1 = 1 + Σ ′ N -1 /X 1 where Σ ′ N -1 := N n=2 X n d = N -1 n=1 X n =: Σ N -1 . By conditioning on X 1 , with f ΣN-1 the density of Σ N -1 , we get P (N ) i,1 = N E Σ β N ∞ 1 dx • f X1 (x) x β+1 ∞ 1+ N -1 x z β-i f ΣN-1 (x (z -1)) dz.
Reversing the integration and after two changes of variables

P (N ) i,1 = N E Σ β N ∞ 1 dz • z β-i ∞ N -1 z-1 x β+1 f X1 (x) f ΣN-1 (x (z -1)) dx = N α E Σ β N ∞ 1 dz • z β-i (z -1) α-β-1 ∞ N -1 s β-α f ΣN-1 (s) ds = N α E Σ β-α N -1 E Σ β N 1 0 u i-2 u 1-α (1 -u) α-β-1 du.
When β < α < 2, the latter integral is (upon adequate normalization by a beta function term B (2α, αβ)) identified with the order i -2 moment of a beta(2α, αβ) rv. In this parameter range, for i ≥ 2

1 0 u i-2 u 1-α (1 -u) α-β-1 du = Γ (i -α) Γ (α -β) Γ (i -β) =: B (i -α, α -β)
are well-defined. ⋄

The Pareto rv X 1 has power-law tails with index α. Therefore E X β 1 only exists

when β < α, with E X β 1 = α/ (α -β) .
Clearly also, the tails of the sum Σ N obey

P (Σ N > s) ∼ s→∞ N P (X 1 > s) = N s -α . Because N P (X 1 > s) ∼ s→∞ P (M N > s) = 1-(1 -s -α )
N where M N = max (X 1 , .., X N ), this means that for large s, the event Σ N > s is essentially determined by the event M N > s. Note that, as a result, Σ N has the same tail index as X 1 , indicating that the β-moment of Σ N only exists for β < α.

We will show below from this, that large N estimates of E Σ β N can be obtained.

As a result, for example, based on (4), whenever 1 < α < 2 and β < α, it will be checked that

c N := P (N ) 2,1 = N α E Σ β-α N -1 E Σ β N B (2 -α, α -β) ∝ N -(α-1) → N →∞ 0.
We will show that, in that case, for each i ≥ 2, the limits lim

N →∞ c -1 N P (N ) i,1 exist, with c -1 N P (N ) i,1 → N →∞ 1 B (2 -α, α -β) Γ (α -β) Γ (i -α) Γ (i -β) = 1 0 u i-2 Λ (du) > 0,
where Λ (du) has the density on [0, 1]:

u 1-α (1 -u) α-β-1 /B (2 -α, α -β), which is a beta(2 -α, α -β) density.
Because in the large N limit, simultaneous multiple collisions will be seen to be negligible, we conclude, using similar arguments to the ones in [START_REF] Sagitov | The general coalescent with asynchronous mergers of ancestral lines[END_REF] and [START_REF] Möhle | A classification of coalescent processes for haploid exchangeable population models[END_REF], that, in the range 1 < α < 2 and β < α, a time-scaled version of the finite discrete-timek coalescent x (N ) k arising from size-biased sampling out of Pareto(α) partitions converges weakly (as

N → ∞) to a continuous-time-t Λ-coalescent x t where Λ is a beta(2 -α, α -β) probability measure on [0, 1] ; namely x t d = lim N →∞ x (N ) [t/cN ] 4 .
In other words, if t is continuous-time, the appropriate scaling is k → [t/c N ] (the integral part of t/c N ), showing that time t should be measured in units of N e := c -1 N (the effective population size). The limiting process x t is a continuous-time pure death Markov process on N = {1, 2, ...}, absorbed at state {1} and with transition rates from i to j given by

c -1 N P (N ) i,j → N →∞ λ i,j := i j -1 1 0 u i-j-1 (1 -u) j-1 Λ (du) , 1 ≤ j < i.
The rate terms λ i,j may also be written as

λ i,j := j-1 l=0 (-1) l l! i j -l -1 q l+i-j+1 ,
where q l := 1 0 u l-2 Λ (du) are the l-moments of u -2 Λ (du) and also the rates of l to 1 mergers λ l,1 .

2.2. Generalities on Λ-coalescents. Whenever one gets a family of rates λ i,j as above for some finite probability measure Λ on [0, 1] , one speaks of continuoustime-t Λ-coalescents (see [START_REF] Pitman | Coalescents with multiple collisions[END_REF] for a precise definition). These Λ-coalescents are non-increasing pure death Markov processes, say x t , on the state-space N, and the λ i,j s are the rates at which i to j < i mergers for x t occur. The state {1} is absorbing. In such processes, multiple collisions of any order (when 1 ≤ j < i) can occur, but never simultaneously. The total death rate at which some merger occurs, starting from state i, is λ i := i-1 j=1 λ i,j . One can check that when Λ = δ 0 (corresponding to the Kingman coalescent), λ i,j = 0 only when j = i -1 with

λ i = λ i,i-1 = i 2 .
The general expression of λ i is When Λ ({0}) = 0 (excluding thereby the Kingman coalescent), the precise dynamics of x t when started at i is given by x 0 = i and [START_REF] Birkner | Computing likelihoods for coalescents with multiple collisions in the infinitely many sites model[END_REF] (5)

λ i = [0,1] u -2 1 -(1 -u) i -iu (1 -u) i-1 Λ (du) .
x tx 0 = -

(0,t]×(0,1] B x s-, u -1 B(xs -,u)>0 N (ds × du) = - (0,t]×(0,1] B x s-, u -1 + N (ds × du) .
Here, x + = max (x, 0), N is a random Poisson measure on [0, ∞)× (0, 1] with intensity ds × 1 u 2 Λ (du) and B x s-, u d ∼ bin x s-, u is a binomial rv with parameters x s-, u . As a result, with [START_REF] Brunet | Effect of selection on ancestry: an exactly soluble case and its phenomenological generalization[END_REF] r (i) :=

(0,1] ui -1 + (1 -u) i u -2 Λ (du) , i > 0,
upon taking the expectation in [START_REF] Brunet | Noisy traveling waves: effect of selection on genealogies[END_REF], it holds that

E dx t | x t-= -r x t-dt.
From this, the quantity r (i) is the rate at which size i blocks are being lost as time passes by. Clearly r is also

r (i) = iλ i - i-1 j=1 jλ i,j = i-1 j=1 (i -j) λ i,j .
Consequently, the reciprocal function 1/r (i) of the rate r (i) interprets as the expected time spent by x t in a state with i lineages and therefore x0=i j=2 1/r (j) will give a large i estimate of the expected time to the most recent common ancestor (the height of the coalescent tree):

τ i,1 := inf (t ∈ R + : x t = 1 | x 0 = i) .
There are lots of detailed studies in the literature (see a precise partial list below) on other functionals of x t such as the total branch length L i of the Λ-coalescent, its total external branch length L e i , the length l i of its external branch (the time till first collision of a branch chosen at random out of i), the number of collisions C i till time to most recent common ancestor,... All these functionals obey some distributional identities which prove useful to obtain some insight on their limit laws as i → ∞. Given x 0 = i, all involve the number U i of singletons taking part in the first collision occurring at time T i d ∼ exp(λ i ) , giving x Ti -1 = i -U i singletons not participating to the first collision (with T i , x Ti independent). The quantity U i is important in itself because, due to

P (U i = j) = λ i,i-j+1 /λ i , j = 2, .., i r (i) = λ i (E (U i ) -1)
where

E (U i ) = i λ i (0,1] 1 -u -(1 -u) i u (1 -u) Λ (du) .
Provided Λ has no atom at point {1} , the condition ∞ i=2 1/r (i) < ∞ is the necessary and sufficient condition for x t to come down from infinity, [START_REF] Schweinsberg | A necessary and sufficient condition for the Λ-coalescent to come down from infinity[END_REF].

Clearly indeed 5 , τ i,1 d = T i +τ xT i ,1 , L i d = iT i +L xT i , L e i d = iT i +L e xT i -1 , C i d = 1+C xT i and l i d = T i + B i l xT i -1
where B i is a Bernoulli rv, given by: P

(B i = 1 | x Ti ) = (x Ti -1) /i with B i l xT i -1 independent of T i .
Famous examples that we shall deal with in the sequel, include Λ-coalescents for which: -

-Λ (du) = B (a, b) u a-1 (1 -u) b-1 1 [0,1] (u) du,
Λ (du) = B (2 -α, α) u 1-α (1 -u) α-1 1 [0,1] (u) du, (α ∈ (0, 1) ∪ (1, 2)); this is the beta(2 -α, α) coalescent.
-(Lebesgue) Λ (du) = 1 [0,1] (u) du : this is the Bolthausen-Sznitman coalescent or beta(1, 1) coalescent.

-Λ (du) = δ 0 : we get the Kingman coalescent where only binary mergers can occur (j = i -1) one at a time, [START_REF] Kingman | The coalescent[END_REF].

-Λ (du) = δ 1 : we get the star-shaped coalescent involving a single big collision. ⋄ Using the above distributional identities, it can be shown for instance that, for the Kingman coalescent and for large i, to leading order of magnitude, rough estimates are:

τ i,1 ∼ 2 (1 -1/i) [40], L i ∼ 2 log i [12], L e i ∼ 2 [22], C i = i -1 and l i ∼ 1/i [9].
For the Bolthausen-Sznitman coalescent:

τ i,1 ∼ log log i [17], L i ∼ i/ log i [12], C i ∼ i/ log i [21] and l i ∼ 1/ log i [15].
For the beta(2α, α) coalescent with 0 < α < 1, L i ∼ i [START_REF] Möhle | On the number of segregating sites for populations with large family sizes[END_REF] and L i /L e i converges in probability to 1 [START_REF] Möhle | Asymptotic results for coalescent processes without proper frequencies and applications to the two-parameter Poisson-Dirichlet coalescent[END_REF] and l i ∼ O (1) [START_REF] Gnedin | On asymptotics of exchangeable coalescents with multiple collisions[END_REF].

For the beta(2 -α, α) coalescent with 1 < α < 2, L i ∼ 1/i α-2 [24], L e i ∼ 1/i α-2 [11] and l i ∼ 1/i α-1 [10].
This information is important to grasp the general shape of the coalescent trees in each case.

GCLT for Pareto sums

In this Section, we first sketch a loose connection of the large N estimate of c N with the Generalized Central Limit Theorem for random variables in the domain of attraction of stable laws, [START_REF] Uchaikin | Chance and stability. Stable distributions and their applications[END_REF].

Let Σ N = N n=1 X n be the partial sums of the i.i.d. Xs with Pareto(α) distribu- tions, α > 0, on (1, ∞). Let S α be skewed α-stable rvs with E e -λSα = e -λ α if α ∈ (0, 1) , λ ≥ 0 5 d
= means equality in distribution between random variables.

the Laplace-Stieltjes transform (LST) of a one-sided α-stable rv,

E e iλS1 = e -|λ|(1+isign(λ) 2 π log|λ|) the characteristic function (c.f.) of a skewed 1-stable Cauchy rv on R E e iλSα = e -|λ| α (1-isgn(λ) tan( πα 2 )) if α ∈ (1, 2) , λ ∈ R
the characteristic function (c.f.) of a skewed α-stable Cauchy rv on R and

E e iλSα = e -λ 2 /2 if α ≥ 2, λ ∈ R the c.f. of a standard normal rv on R.
The following Generalized Central Limit Theorem (GCLT) then holds Theorem 2. ( [START_REF] Uchaikin | Chance and stability. Stable distributions and their applications[END_REF], [START_REF] Zaliapin | Approximating the distribution of Pareto sums[END_REF]): Let Σ N denote the partial sum sequence of N i.i.d. Pareto(α) random variables. Then,

(7) Σ N -a N b N d → N →∞ S α
where, with

C α = Γ (1 -α) cos πα 2 1/α if α ∈ (0, 2) \ {1} , C 1 = π 2 , C α = α α -2 - α α -1 2 1/2 if α > 2, b N is given by (i) b N = C α N max(1/α,1/2) if α = 2 (ii) b N = (N log N ) 1/2 if α = 2
and, with γ the Euler constant and E (X 1 ) = µ := α/ (α -1), a N is given by

(i) a N = 0 if α ∈ (0, 1) (ii) a N = N π 2 2 ∞ 1 sin 2x πN dF X1 (x) ∼ N log N + N 1 -γ -log 2 π if α = 1 (iii) a N = N µ if α ∈ (1, ∞) .
When α ≤ 1, the characteristic values of Σ N can be guessed to be what they are claimed to be (N

1/α if α < 1 and N log N if α = 1) while estimating N mN 1 xf X1 (x)
dx where m N is the mode of M N which is seen to grow like N 1/α (see [START_REF] Bouchaud | Anomalous diffusion in disordered media: statistical mechanisms, models and physical applications[END_REF]) and similarly for the fluctuation scaling term b N for α ≤ 2 (of order

N 1/α if α < 2 and (N log N ) 1/2 if α = 2).
From these rough estimates, we would conclude that with c N ∝ N

E(Σ β-α N -1 ) E(Σ β N ) , up to the leading order in N c N ∝ N C α N 1/α β-α C α N 1/α β = O (1) if α ∈ (0, 1) c N ∝ N (N log N ) β-1 (N log N ) β ∼ 1 log N if α = 1 c N ∝ N (µN ) β-α (µN ) β ∼ µ -α N -(α-1) if α ∈ (1, 2) .
Depending on the values of α, we therefore anticipate

• α ∈ (0, 1): Because in that case c N is asymptotic to a constant, this suggests a limiting discrete-time coalescent. We will show below that it is not a discrete-time Λ-coalescent, rather it is a Ξ-coalescent of the Poisson-Dirichlet type with two parameters (α, -β). Ξ-coalescents were first introduced in [START_REF] Möhle | A classification of coalescent processes for haploid exchangeable population models[END_REF] and further studied in [START_REF] Schweinsberg | Coalescents with simultaneous multiple collisions[END_REF]. In sharp contrast with Λ-coalescents, multiple collisions can occur simultaneously at the same transition time. In their block-counting version, they are characterized by the set of numbers φ j (i 1 , .., i j ) defining the probabilities of an (i 1 , .., i j ) -merger (i 1 + .. + i j = i), resulting when the Ξ-coalescent is discrete, in an i to j ≤ i transition with probability

P i,j = 1 j! * i1+...+ij =i i i1...ij φ j (i 1 , .
., i j ). The φ j (i 1 , .., i j ) can be written as

φ j (i 1 , .., i j ) = ∆j j l=1 u i l -2 l Λ j (du 1 , .., du j ) ,
for some finite measures Λ j with density on the (j + 1)simplex

∆ j = (u 1 , .., u j ) ∈ [0, 1] j : u 1 + ... + u j ≤ 1 .
The set of measures Λ j , j ≥ 1, (characterized by their moments φ j ), with values over the simplices ∆ j , completely characterize the Ξ-coalescent, [START_REF] Möhle | A classification of coalescent processes for haploid exchangeable population models[END_REF].

In the simplest cases, with u, u := l≥1 u 2 l , it is also convenient (see [START_REF] Schweinsberg | Coalescents with simultaneous multiple collisions[END_REF]) to rewrite the φ j s as

φ j (i 1 , .., i j ) = ∆ * k1,...,kj all distinct j l=1 u i l k l Ξ (du) u, u ,
where ∆ = u:= (u 1 , .., u l , ..) : u 1 ≥ .. ≥ u l ≥ .. ≥ 0 : l≥1 u l ≤ 1 and Ξ a finite measure concentrated on the subset ∆ * of ∆ consisting of those us exactly summing to 1. Letting ν (du) := Ξ (du) / u, u , the measure ν on the infinite simplex ∆ is such that ν (∆) < ∞.

• α = 1: A Λ-coalescent with logarithmic effective population size where Λ is a beta(1, 1β) probability measure with β < 1 as in [START_REF] Brunet | Genealogies in simple models of evolution[END_REF] (reducing to the Bolthausen-Sznitman coalescent when avoiding size-biasing corresponding to β = 0).

• α ∈ (1, 2): A Λ-coalescent with N e ∝ N α-1
where Λ is a beta(2α, αβ) measure with β < α. In the latter case, β = 0 leads to the standard beta(2α, α) coalescent.

• α ≥ 2 : In the range α > 2, (respectively α = 2), N e ∝ N (respectively N e ∝ N/ log N ) and the obtained scaled limiting coalescent is a Kingman coalescent, as a result of Λ approaching δ 0 . In this latter case, only binary mergers occur in the large N limit.

Remark: The Kingman coalescent also occurs when dealing with sampling from some alternative random partition. For instance, would sampling be defined from a random partition of unity given by S n := X n / N 1 X n , n = 1, ..., N where X 1 now obeys the following gamma(θ) density:

f X1 (x) = Γ (θ) -1 x θ-1 e -x , θ, x > 0, then the law of Σ N = X 1 + ... + X N is f ΣN (x) = Γ (N θ) -1 x N θ-1 e -x
, independent of S 1 and

P (N ) i,1 = E Σ β N S i n E Σ β N = N E Σ β N S i Γ (i + θ) Γ (N θ + i) ∼ N -(i-1) Γ (i + θ) Γ (θ) θ -i . Thus c N = P (N ) 2,1 = 1 N 1+θ θ → N →∞ 0 together with d N = P (N ) 3,1 = 1 N 2 (1+θ)(2+θ) θ 2
. Because triple mergers are asymptotically negligible compared to binary ones (d N /c N → 0), the time-scaled limiting coalescent using

c N = 1 N 1+θ θ ∝ N -1 is a Kingman coalescent. The prefactor 1+θ θ appearing in front of c N is the ratio ρ/µ 2 where ρ := E X 2 1 = θ (θ + 1) and µ := E (X 1 ) = θ and c N is independent of β. 4. Large N asymptotic estimation of E Σ β N

and consequences

In this Section, we compute the asymptotic behavior of the β-moments of Σ N in the cases α ∈ (0, 1) , α ∈ (1, 2) and α = 1, and α ≥ 2, making the previous conclusions based on the GCLT consistent 6 .

• We start with the case α ∈ (0, 1) .

Theorem 3. When α ∈ (0, 1) , as N → ∞, x k d = lim N →∞ x (N ) k
exists and is a discrete-time Poisson-Dirichlet(α, -β) Ξ-coalescent.

Proof: For the values of β for which it makes sense, we have

E Σ β N = 1 Γ (-β) ∞ 0 dλ•λ -β-1 E e -λΣ N = 1 Γ (-β) ∞ 0 dλ•λ -β-1 E e -λX1 N .
When N is large, only the small λ approximation of E e -λX1 to the latter integral contributes. For small λ, we have ( 8)

E e -λX1 = α ∞ 1 x -(α+1) e -λx dx = 1 -α ∞ 1 x -(α+1) 1 -e -λx dx ∼ 1 -α ∞ 0 x -(α+1) 1 -e -λx dx ∼ 1 -Γ (1 -α) λ α ∼ e -Γ(1-α)λ α .
6

The technique we use is inspired from the one used in [START_REF] Brunet | Genealogies in simple models of evolution[END_REF] in a particular case. The author is indebted to B. Derrida for pointing this out to him.

Note that, when λ is small

E X i 1 e -λX1 = (-1) i d i dλ i E e -λX1 ∼ αΓ (i -α) λ α-i . Thus, E Σ β N ∼ 1 Γ (-β) ∞ 0 dλ • λ -β-1 e -N Γ(1-α)λ α .
With the change of variables u = N Γ (1α) λ α , with β < α, we get

E Σ β N ∼ N β/α αΓ (-β) Γ (1 -α) β/α ∞ 0 du•u -β/α-1 e -u = N β/α Γ (1 -α) β/α Γ (1 -β/α) Γ (1 -β) .
Finally, using the above large N estimate of E Σ β N , the identity Γ (x + 1) = xΓ (x) and ( 4) with i = 2, when β < α < 1, we obtain

= Γ (1 -β) Γ (α -β) Γ (1 -α) Γ (2 -α) Γ (α -β) Γ (2 -β) = 1 -α 1 -β =: c ∞ .
Thus the coalescence probability c N converges to c ∞ ∈ (0, 1) . When α ∈ (0, 1), using again (4), we obtain more generally

P (N ) i,1 → N →∞ P i,1 = Γ (1 -β) Γ (1 -α) Γ (i -α) Γ (i -β) ,
which are the probabilities to merge all i particles in one step in the limiting discretetime coalescent.

To derive the full transition probabilities of the limiting discrete-time-k coalescent

x k d = lim N →∞ x (N )
k , recalling that

P (N ) i,j = N j * i1+...+ij =i i i 1 ...i j E Σ β N j l=1 S i l l E Σ β N , we use (9) E Σ β N j l=1 S i l l = 1 Γ (i -β) ∞ 0 dλ • λ i-β-1 j l=1 E X i l l e -λX l E e -λX1 N -j ∼ α j Γ (i -β) j l=1 Γ (i l -α) ∞ 0 dλ • λ i-β-1 λ αj-i e -N Γ(1-α)λ α .

Performing again the change of variables

u = N Γ (1 -α) λ α , with β < α, we get E Σ β N j l=1 S i l l ∼ α j-1 N β/α-j Γ (1 -α) β/α-1 Γ (j -β/α) Γ (i -β) j l=1 Γ (i l -α) .
Using N j ∼ N j /j! for large N , we finally get P These are the full transition probabilities of the limiting discrete-time-k coalescent x k in the regime α ∈ (0, 1) and β < α (satisfying i j=1 P i,j = 1). Note that P i,1 are the probabilities obtained previously and that the diagonal terms (the eigenvalues of P ) read

P i,i = α i-1 Γ (1 -β) Γ (1 -β/α) Γ (i -β/α) Γ (i -β) .
Clearly this discrete-time coalescent is not a discrete Λ-coalescent as simultaneous multiple collisions can occur (it is a Ξ-coalescent). Clearly P i,j is also

P i,j = 1 j! j l=1 Γ ((l -1) α + 1 -β) Γ (1 -α) Γ (lα -β) .
Thus ( 10) is also [START_REF] Dhersin | Asymptotic behavior of the total length of external branches for Beta-coalescents[END_REF])

P i,j = c j,α,β i! j! Γ (αj -β) Γ (i -β) * i1+...+ij=i j l=1 Γ (i l -α) i l ! .
Defining the finite Dirichlet measures Λ j with density on the (j + 1)simplex ∆ j given by:

Λ j (du 1 , .., du j ) = c j,α,β c ∞ j l=1 u 1-α l du l 1 - j l=1 u l αj-β-1
, we get

φ j (i 1 , .., i j ) = c ∞ ∆j j l=1 u i l -2 l Λ j (du 1 , .., du j ) .
The set of finite Dirichlet measures Λ j with parameters

(θ 1 = 2 -α, ..., θ j = 2 -α, θ j+1 = αj -β)
on the simplices ∆ j completely characterize this limiting discrete-time coalescent. Note that Λ 1 is a beta(2α, αβ) probability measure.

One may rewrite the φ j s as (see [START_REF] Schweinsberg | Coalescents with simultaneous multiple collisions[END_REF] and [START_REF] Möhle | Asymptotic results for coalescent processes without proper frequencies and applications to the two-parameter Poisson-Dirichlet coalescent[END_REF])

φ j (i 1 , .., i j ) = ∆ * k1,...,kj all distinct j l=1 u i l k l Ξ (du) u, u ,
where ∆ = u:= (u 1 , .., u l , ..) : u 1 ≥ .. ≥ u l ≥ .. ≥ 0 : l≥1 u l ≤ 1 and Ξ a measure on ∆. Letting ν (du) := Ξ (du) / u, u , the measure ν on the infinite simplex ∆ can be identified (see [START_REF] Möhle | Asymptotic results for coalescent processes without proper frequencies and applications to the two-parameter Poisson-Dirichlet coalescent[END_REF]) to the two-parameter Poisson-Dirichlet(α, -β) measure, with α ∈ [0, 1) and β < α. It holds that ν (∆) = 1. Poisson-Dirichlet measures enjoy many remarkable properties including a stick-breaking property, [START_REF] Pitman | The two parameter Poisson-Dirichlet distribution derived from a stable subordinator[END_REF]. ⋄ Remarks:

(i) From ( 11), the limiting situation α = 0 also makes sense, leading to the oneparameter Poisson-Dirichlet(0, -β) measure with β < α = 0. In this case, from

P i,j = (-β) j Γ (-β) Γ (i -β) s i,j , (11)
where s i,j

:= i! j! * i1+...+ij =i j l=1 1
i l are the absolute first kind Stirling numbers.

(ii) Finally, avoiding size-biasing (β = 0) gives rise to the discrete Poisson-Dirichlet(α, 0) coalescent with one-parameter α ∈ (0, 1) , appearing in [START_REF] Schweinsberg | Coalescent processes obtained from supercritical Galton-Watson processes[END_REF].

• The case α ∈ (1, 2) . Lemma 4. When α ∈ (1, 2) , with µ := α α-1 , for all β < α (12) c N := P (N ) 2,1 ∼ αµ -α B (2 -α, α -β) N -(α-1) → N →∞ 0.
Proof: In this case, defining a := α -1 ∈ (0, 1) , after an integration by parts and using the previous estimate (8) substituting a to α

E e -λX1 = e -λ - λ a a ∞ 1 x -(a+1) e -λx dx ∼ e -λ - λ a (1 -Γ (1 -a) λ a ) .
Thus, for small λ, with µ := α α-1

(13) E e -λX1 ∼ 1 -λµ -Γ (1 -α) λ α ∼ e -λµ .
Thus, consistently with the GCLT approach, to the dominant order 1) . ⋄ This suggests that Proposition 5. When α ∈ (1, 2) , upon scaling time using an effective population size N e = c -1 N (with c N as in ( 12)), we obtain a limiting continuous-time-t Λ-coalescent:

E Σ β N = 1 Γ (-β) ∞ 0 dλ • λ -β-1 E e -λX1 N ∼ (N µ) β , so that c N := N α E Σ β-α N -1 E Σ β N B (2 -α, α -β) ∼ αµ -α B (2 -α, α -β) N -(α-
x t d = lim N →∞ x (N ) [t/cN ] , with Λ a beta(2 -α, α -β) probability mea- sure, β < α. If β = 0, we get the standard beta(2 -α, α) coalescent.
Proof: To confirm this point, we will first evaluate a large N estimate of P (N ) i,1 defined in [START_REF] Bertoin | Lévy processes with no negative jumps and branching processes[END_REF], in the range α ∈ (1, 2) . Using the above small λ estimate of E e -λX1 in the parameter range under concern, for i ≥ 2, we get

E X i 1 e -λX1 = (-1) i d i dλ i E e -λX1 ∼ αΓ (i -α) λ α-i . Thus, using (9) E Σ β N S i 1 = 1 Γ (i -β) ∞ 0 dλ • λ i-β-1 E X i 1 e -λX1 E e -λX1 N -1 ∼ α Γ (i -α) Γ (i -β) ∞ 0 dλ • λ i-β-1 λ α-i e -µλN = α Γ (i -α) Γ (i -β) Γ (α -β) (µN ) α-β .
Finally, we obtain

P (N ) i,1 ∼ N -(α-1) α µ α Γ (α -β) Γ (i -α) Γ (i -β) ,
showing that, with c N = P

(N ) 2,1 ∝ N -(α-1) → 0, for each i, lim N →∞ c -1 N P (N )
i,1 exist and are strictly positive constants. More precisely,

c -1 N P (N ) i,1 → N →∞ φ 1 (i) = 1 B (2 -α, α -β) Γ (α -β) Γ (i -α) Γ (i -β) = 1 0 u i-2 Λ 1 (du) ,
where Λ 1 = Λ is a beta(2α, αβ) probability measure.

To deal with the higher order terms, with m := {l ∈ {1, ..., j} : i l ≥ 2}, assuming 1 ≤ j < i, let us write (3) as

P (N ) i,j = N j j m=1 j m * * i1+...+im=i-j+m i i 1 ...i m E Σ β N m l=1 S i l l j l=m+1 S l E Σ β N ,
where the double-star sum is now over the i l , l = 1, ..., m satisfying i l ≥ 2. Proceeding similarly with higher order terms, it clearly holds that for all j ≥ 2 and i l ≥ 1, l = 1, ..., j, satisfying i l ≥ 2 for at least two l in the list, lim

N →∞ c -1 N E Σ β N j l=1 S i l l /E Σ β N = 0,
so that simultaneous multiple collisions cannot occur in the limit (actually the contribution of these simultaneous multiple collisions terms in P

(N ) i,j is of order O (c m N )
). In the latter expression of P (N ) i,j therefore, only the term corresponding to m = 1 will contribute to the

O (c N ) -order. Because, E X 1 e -λX1 = -d dλ E e -λX1 ∼ µ and E X i1 1 e -λX1 ∼ αΓ (i 1 -α) λ α-i1 with i 1 = i -j + 1 ≥ 2, using (9)
, we get

E Σ β N S i1 1 j l=2 S l E Σ β N = ∞ 0 dλ•λ i-β-1 E X i1 1 e -λX1 j l=2 E X l e -λX l E e -λX1 N -j E Σ β N Γ (i -β) ∼ αµ j-1 Γ (i 1 -α) E Σ β N Γ (i -β) ∞ 0 dλ•λ i-β-1 λ α-i1 e -N µλ = αµ j-1 Γ (i 1 -α) E Σ β N Γ (i -β) Γ (α -β + j -1) (N µ) α-β+j-1 .
This shows, using N j ∼ N j /j!, [START_REF] Drmota | Asymptotic results concerning the total branch length of the Bolthausen-Sznitman coalescent[END_REF] and after some elementary algebra, that for j = 1, ..., i -1,

c -1 N P (N ) i,j → N →∞ λ i,j = i j -1 B (i -j + 1 -α, α -β + j -1) B (2 -α, α -β) = i j -1 1 0 u i-j-1 (1 -u) j-1 Λ (du) , where Λ = Λ 1 d ∼ beta(2-α, α-β). As a result, x t d = lim N →∞ x (N )
[t/cN ] is a continuoustime-t pure death coalescent process on N with infinitesimal transition rates λ i,j . ⋄

• The case α = 1. We view it as a limiting case of the previous analysis deriving from [START_REF] Wakeley | Coalescent processes when the distribution of offspring number among individuals is highly skewed[END_REF] when α → 1 + . Proposition 6. When α = 1, upon scaling time using a logarithmic effective population size N e = c -1 N ∼ log N, the limiting process

x t d = lim N →∞ x (N )
[t/cN ] exists and is a continuous-time-t Λ-coalescent with Λ a beta(1, 1-β) probability measure, β < 1. If β = 0, we get the standard Bolthausen-Sznitman coalescent with Λ uniform.

Proof: Put indeed α = 1 + ε in the previous small-λ estimate ( 13) of E e -λX1 in the parameter range α ∈ (1, 2), with ε > 0 small. Then µ ∼ 1 + 1/ε and because

Γ (1 -α) = Γ (-ε) ∼ 0 + -1/ε, (14) E e -λX1 ∼ 1 -λ - 1 ε λ -λ 1+ε ∼ 1 -λ + λ log λ =: I (λ) .
Thus,

E Σ β N ∼ 1 Γ (-β) ∞ 0 dλ • λ -β-1 I (λ) N .
The leading contribution of I (λ) N is, when λ is small, of order 1/ (N log N ). Putting λ = u/ (N log N )

I (λ) N ∼ e -u 1 + u log N (log u -log log N -1) .
Thus, with β < 1

E Σ β N ∼ (N log N ) β Γ (-β) ∞ 0 du • u -β-1 e -u 1 + u log N (log u -log log N -1) (15) ∼ (N log N ) β 1 + β log N (ψ (1 -β) -log log N -1) ,
where ψ (x) = Γ ′ (x) /Γ (x) is the digamma function. In the latter estimate, we used that differentiating

∞ 0 du•u -β u θ e -u = Γ (1 -β + θ)
with respect to the extra parameter θ and then putting θ = 0 gives

∞ 0 du • u -β log (u) e -u = Γ ′ (1 -β).
Finally, consistently with the GCLT approach

c N = N E Σ β-1 N -1 E Σ β N ∼ 1 log N .
Clearly, with Λ a beta(1, 1β) probability measure and β < α = 1

c -1 N P (N ) i,1 → N →∞ 1 0 u i-2 Λ (du) = Γ (2 -β) Γ (i -1) Γ (i -β) = λ i,1 .
Using similar arguments on higher order terms, showing that simultaneous multiple collisions do not contribute in the limit, one can easily show

(16) c -1 N P (N ) i,j → N →∞ 1 0 u i-j-1 (1 -u) j-1 Λ (du) = B (i -j, j -β) B (1, 1 -β) =: λ i,j , j = 1, ..., i -1, where Λ = Λ 1 d ∼ beta(1, 1 -β).
This confirms that, when α = 1, upon scaling time using an effective logarithmic population size Proof: If α > 2, Σ N is in the domain of attraction of the normal law. As a result, for small λ, with µ := E (X 1 ) = α α-1 and ρ

N e = c -1 N ∼ log N, we get a limiting continuous-time-t Λ-coalescent x t d = lim N →∞ x (N ) [t/cN ] , where is a beta(1, 1 -β) probability measure, β < 1. ⋄ • The case α > 2.
:= E X 2 1 = α α-2 (17) E e -λX1 ∼ 1 -λµ + 1 2 ρλ 2 ∼ e -λµ .
From [START_REF] Goldschmidt | Random recursive trees and the Bolthausen-Sznitman coalescent[END_REF], for small λ, E X 1 e -λX1 = -d dλ E e -λX1 ∼ µ and E X 2 1 e -λX1 = d 2

dλ 2 E e -λX1 ∼ ρ.

We have c N := P

(N ) 2,1 = N E(Σ β N S 2
1 )

E(Σ β N )
with, from (9)

E Σ β N S 2 1 ∼ ρ Γ (2 -β) ∞ 0 dλ • λ 1-β e -λµN = ρ (µN ) β-2 E Σ β N ∼ 1 Γ (-β) ∞ 0 dλ • λ -β-1 e -λµN = (µN ) β .
Thus,

c N = N E Σ β N S 2 1 E Σ β N ∼ ρ µ 2 N -1
goes to 0 as claimed. As in Proposition 5, simultaneous multiple collisions cannot contribute in the limit. Only the term in the expression of P (N ) i,j corresponding to m = 1 will contribute to the O (c N ) -order and so we need to focus on a collision with a single i 1 ≥ 2. Now, from ( 17), E X 1 e -λX1 = -d dλ E e -λX1 ∼ µ and E X i1 1 e -λX1 ∼ ρδ i1,2 with i 1 = ij + 1 ≥ 2. Proceeding again as in Proposition 5 therefore, only effective transitions from i to j < i with j = ii 1 + 1 = i -1 are seen in the limit, corresponding to the binary mergers of a Kingman coalescent. ⋄ • α = 2. To complete the picture, it remains to study the limiting critical case α = 2. Lemma 8. When α = 2, for all values of β, the coalescence probability goes to 0 as N → ∞ like

c N ∼ 1 2 log N N .
Proof: We view the case α = 2 as a limiting case of ( 17) as α → 2 + . When λ is small and for α > 2, we have

E e -λX1 ∼ 1 -λµ + 1 2 ρλ 2 -Γ (1 -α) λ α .
Putting α = 2 + ε, for ε > 0 small, using Γ (-1ε) ∼ 0 + 1/ε, with ρ ∼ 2/ε, when λ is small, we have

(18) E e -λX1 ∼ 1 -2λ + 1 ε λ 2 - 1 ε λ 2+ε ∼ 1 -2λ -λ 2 log λ ∼ e -2λ . Because E X 1 e -λX1 = -d dλ E e -λX1 ∼ 2 and E X 2 1 e -λX1 = d 2 dλ 2 E e -λX1 ∼ -2 log λ -3, we have E Σ β N S 2 1 ∼ 1 Γ (2 -β) ∞ 0 dλ • λ 1-β E X 2 1 e -λX1 E e -λX1 N -1 ∼ ∼ -1 Γ (2 -β) ∞ 0 dλ • λ 1-β (2 log λ + 3) e -2N λ E Σ β N ∼ 1 Γ (-β) ∞ 0 dλ • λ -β-1 e -2λN = (2N ) β .
The Euler integral with a logarithmic term inside appearing in the expression of E Σ β N S 2 1 can be obtained while taking the derivative of ∞ 0 dλ•λ 1-β λ θ e -2N λ with respect to the extra parameter θ and then putting θ = 0 in the obtained expression. Observing therefore Proof: Using ( 18), for small λ, we have

∞ 0 dλ•λ 1-β log (λ) e -2N λ = Γ ′ (2 -β) (2N ) -(2-β) -Γ (2 -β) (2N ) -(2-β) log (2N ) , to the dominant order in N , E Σ β N S 2 1 ∼ 2 (2N ) -(2-β) log N, leading to c N := P (N ) 2,1 = N E Σ β N S 2 1 E Σ β N ∼ 1
E X i 1 e -λX1 = (-1) i d i dλ i E e -λX1 ∼ 2Γ (i -2) λ -(i-2) , i ≥ 3, to which one should add E X 1 e -λX1 ∼ 2 and E X 2 1 e -λX1 ∼ -(2 log λ + 3) . For i ≥ 3, we get E Σ β N S i 1 = 1 Γ (i -β) ∞ 0 dλ • λ i-β-1 E X i 1 e -λX1 E e -λX1 N -1 ∼ 2Γ (i -2) Γ (i -β) ∞ 0 dλ • λ 1-β e -2N λ = 2 B (i -2, 2 -β) (2N ) 2-β .
Thus, for all i ≥ 3, as

N → ∞ c -1 N P (N ) i,1 = N c -1 N E Σ β N S i 1 E Σ β N ∼ 1 log N B (i -2, 2 -β) → 0.
Therefore, due to the extra factor log N appearing in c N , the transitions from i ≥ 3 to 1 cannot be seen in the limit, nor (for the same reason) the transitions involving a single multiple collision with i 1 ≥ 3, nor transitions involving simultaneous multiple collisions of any order. In fact, only the events involving a single multiple collision with i 1 = ij + 1 = 2 (corresponding to transitions from i to j = i -1) will contribute in the limit. Indeed,

E Σ β N S 2 1 j l=2 S l E Σ β N = ∞ 0 dλ•λ i-β-1 E X 2 1 e -λX1 j l=2 E X l e -λX l E e -λX1 N -j E Σ β N Γ (i -β) ∼ -2 j-1 E Σ β N Γ (i -β) ∞ 0 dλ • λ i-β-1 (2 log λ + 3) e -2N λ ∼ 1 2 N -i log N.
This shows, using

P (N ) i,j ∼ N j j 1 i 2 1...1 E Σ β N S 2 1 j l=2 S l E Σ β N with j = i -1, N j ∼ N j /j!,
and after some elementary algebra, that

c -1 N P (N ) i,i-1 → N →∞ λ i,i-1 = i 2 . ⋄
Remarks:

(i) Results of a similar flavor can be found in Schweinsberg's work [START_REF] Schweinsberg | Coalescent processes obtained from supercritical Galton-Watson processes[END_REF]. However, his model and techniques are different from ours because he considers large-N limiting coalescents obtained while sampling without replacement from a discrete super-critical Galton-Watson branching process, assuming the reproduction law of each offspring to exhibit power-law Zipf tails of index α (a in his notations). Note that there is no parameter β in the construction [START_REF] Schweinsberg | Coalescent processes obtained from supercritical Galton-Watson processes[END_REF]. In the same spirit, results can also be found in Huillet-Möhle [START_REF] Huillet | On the extended Moran model and its relation to coalescents with multiple collisions[END_REF] where, following [START_REF] Wakeley | Coalescent processes when the distribution of offspring number among individuals is highly skewed[END_REF], Λ-coalescents are obtained as scaling limits of discrete extended Moran models, the skewed reproduction law of which displaying occasional extreme events with one individual allowed to produce a large amount of offspring. Whenever the reproduction law displays systematic extreme events, discrete coalescents were even shown to emerge in the large N -limit, but in this Moran context, they are only Λ-coalescents with multiple but no simultaneous collisions, [START_REF] Huillet | On the extended Moran model and its relation to coalescents with multiple collisions[END_REF]. This contrasts with the occurrence in our present work of discrete Poisson-Dirichlet Ξ-coalescents. In contrast also with our current work where coalescents are derived from a sampling procedure in the continuum, the coalescents considered in [START_REF] Schweinsberg | Coalescent processes obtained from supercritical Galton-Watson processes[END_REF] and [START_REF] Huillet | On the extended Moran model and its relation to coalescents with multiple collisions[END_REF] were built from discrete reproduction laws at fixed population size N and looking at their scaling limits N → ∞.

We refer to these two works and to [START_REF] Brunet | Genealogies in simple models of evolution[END_REF] for additional background on Λ-coalescent processes adapted to our purposes.

(ii) When α ≥ 2, neither the (large N estimate of the) scaling constant c N nor the limiting (Kingman) coalescent depend on the bias parameter β.

(iii) We observe that the probability c N := P (N )

2,1 that two individuals chosen at random share the same common ancestor, defining the time-scale to derive the large-N limits of the Pareto-coalescents all have the same large-N order of magnitude as the length l N of an external branch chosen at random in the limiting Λ-coalescents x t (α ∈ [1, 2)) or Ξ-coalescent x k (α ∈ (0, 1)), or Kingman coalescent x t (α > 2), started at x 0 = N . This curious fact is unexplained so far.

Forward in time selection model and genealogies

In this Section, in the spirit of [START_REF] Brunet | Genealogies in simple models of evolution[END_REF], we indicate that the coalescent processes just discussed may be viewed as the genealogical processes of some forward in time evolving branching population models with selection. As it is often the case in population genetics, the process we are interested in is in the class of branching processes conditioned on having a fixed population size over each generation, in the spirit of [START_REF] Karlin | Direct product branching processes and related Markov chains[END_REF].

A Poisson-point process model with selection.

Start with N individuals at generation t = 0 and assume that each individual has an initial fitness x n (0) > 0, n = 1, ..., N.

To describe the state of the population at the next generation, assume first that, independently of one another, each individual potentially generates an infinite number of offspring along a Poisson point process (PPP) with intensity (or occupation) density [START_REF] Huillet | On the extended Moran model and its relation to coalescents with multiple collisions[END_REF] π xn(0

) (x) = -π ′ xn(0) (x) = αx n (0) α x -(α+1) , 7 7
The ′ symbol indicates derivative with respect to x.

where π xn(0) (x) := (x/x n (0)) -α , x > 0, α > 0, n = 1, .., N. We observe that, with

π (x) := x -α π xn(0) (x) = π (x) π (x n (0)) and if we let π (x) = -π ′ (x) = αx -(α+1) , then π xn(0) (x) = x n (0) α π (x) .
So the fitnesses of the offspring of each individual is generated according to a PPP depending on the fitness of its parent. From these simple assumptions, and as conventional wisdom suggests, we get: Proposition 10. In a PPP model for fitness-dependent offspring reproduction with occupation density [START_REF] Huillet | On the extended Moran model and its relation to coalescents with multiple collisions[END_REF], the fitnesses of the offspring of some parental individual with fitness x n (0) are x n (0) times the fitnesses of the offspring of some canonical individual with unit fitness: In this sense, the larger the fitness x n (0) of some individual is, the more he will, proportionally, produce offspring with large fitness.

Proof: Let (τ n ; n ≥ 1) be the points of a standard homogeneous Poisson point process (PPP) on the half-line with rate 1, and let π -1 (s) = s -1/α be the decreasing inverse of π. Then, with π -1 xn(0

) (s) = x n (0) • s -1/α , π -1 xn(0) (τ n )
; n ≥ 1 are the (ordered) points of the offspring PPP on the positive half-line (or here the fitness space) with occupation density π xn(0) . So, the fitter the individuals, the fitter their offspring, proportionally to the parental fitness. ⋄ Let us also briefly emphasize that, avoiding the fitness dependence on the parent of the PPP, π -1 (τ n ) ; n ≥ 1 are just the (ordered) points of a PPP on the half-line with occupation density π. Whenever, as in our case study here, the rate function π is not integrable up to x = 0, there are infinitely many such points, with 0 as an accumulation point whereas there is of course a finite Poisson (with mean π (ε)) number of them above some threshold ε > 0.

It is well-known that when α ∈ (0, 1), with π (x) := x -α and π -1 (s) = s -1/α , the positive cumulative rv [START_REF] Huillet | Asymptotics of symmetric compound Poisson population models[END_REF] [START_REF] Lamperti | Continuous state branching processes[END_REF]. When α = 1, the law of χ is the one of a positive Neveu rv with LST E e -λχ = e λ log λ (see [START_REF] Neveu | A continuous state branching process in relation with the GREM model of spin glass theory[END_REF]); the latter may be viewed as a Lamperti rv in the limit α → 1 + , [START_REF] Huillet | Energy cascades as branching processes with emphasis on Neveu's approach to Derrida's random energy model[END_REF].

χ d = n≥1 π -1 (τ n ) is a one-sided α-stable rv on (0, ∞) with LST E e -λχ = e -κλ α , κ = Γ (1 -α) > 0, λ ≥ 0. When α ∈ (1, 2), the law of χ is the one of a positive Lamperti rv with LST E e -λχ = e -cλ+κλ α , κ = -Γ (1 -α) > 0, c := E (χ) > 0,
Typically indeed, by the Lévy-Khintchine formula, π (x) dx stands for the Lévy measure for the non-negative jumps of χ which is the value at time t = 1 of an infinitely divisible subordinator (χ t ; t ≥ 0) with LST E e -λχ t , [START_REF] Bertoin | Lévy processes with no negative jumps and branching processes[END_REF]. In [START_REF] Huillet | Asymptotics of symmetric compound Poisson population models[END_REF], the terms π -1 (τ n ) are thus the ranked jumps of χ in its Lévy decomposition (see [START_REF] Perman | Order statistics for jumps of normalized subordinators[END_REF] for example).

Reproduction step.

Because we consider the offspring of all the N initial individuals with fitnesses x n (0) , n = 1, ..., N , the step-1 state of the whole population is thus obtained from a PPP with global equivalent occupation density

(21) π xN,α(0) (x) = αx -(α+1) N n=1 x n (0) α = αx N,α (0) α x -(α+1) ,
where

x N,α (0) := N n=1
x α n (0)

1/α
is the global equivalent initial fitness of the whole population at generation 0 to consider. We get Proposition 11. In a PPP with occupation density [START_REF] Huillet | On the extended Moran model and its relation to coalescents with multiple collisions[END_REF] for the descent of each individual with fitness x n (0) , n = 1, ..., N , the occupation density of the population as a whole is given by ( 21), where x N,α (0) :

= N n=1 x α n (0) 1/α
is the global equivalent fitness.

Proof: This follows from the superposition principle of Poisson point processes (see [START_REF] Kingman | Poisson processes[END_REF] p. 16). The fact that the intensity of the superposed PPP is in the same class as the one of a single PPP descending from x n (0) is a remarkable scaling property of π xn(0) (x) . ⋄ In this setup therefore, π xN,α(0) (x) := -π ′ xN,α(0) (x) stands for the occupation density that there is a point at x descending from any of the N initial individuals with fitnesses x n (0) , n = 1, ..., N.

Note that the cumulated fitness of all first-generation offspring is

n≥1 π -1 xN,α(0) (τ n ) d = x N,α (0) • χ,
where χ is given by [START_REF] Huillet | Asymptotics of symmetric compound Poisson population models[END_REF].

Selection step.

In order to model a population with fixed size over the generations, the final state of the population at time 1 is obtained while selecting the N individuals of the whole population whose fitnesses are the largest (the selection step), truncating therefore the latter sum to its N first terms.

The whole process (including reproduction and selection steps) is then iterated independently over the next generations.

From this definition of the process, if the x n (k)s are the fitnesses of the N fittest individuals at generation k, the ordered ones x (n) (k + 1) at generation k +1 (x (1) > ... > x (N ) ) descending from the whole population at step k are given by x

(n) = π -1 xN,α(k) (τ n ) = x N,α (k) τ -1/α n
. Here the τ n s are the ordered points of a standard Poisson process on the half-line with τ 1 < ... < τ n < ... < τ N . The law of τ n is thus the one of an Erlang gamma(n) rv with density f τ n (s) = s n-1 e -s /Γ (n) . Further, given τ N +1 = s, the probability density of τ 1 , ..., τ N is The joint law of the ordered x n (k + 1)s is thus the one of the images π -1 xN,α(k) (τ n )s, namely

f τ 1 ,...,τN (s 1 , .., s N | τ N +1 = s) = N ! s N 1 0<s1<...
f x (1) (k+1),...,x (N ) (k+1) (x 1 , .., x N ) = e -π x N,α (k) (xN ) N n=1 π xN,α(k) (x n ) 1 x1>...>xN .
Clearly also (with the two terms in the right-hand side term mutually independent),

x (N +1) (k + 1) d = x N,α (k) x * N +1 (k + 1)
, where x * N +1 (k + 1) is the (N + 1) -st largest point of a PPP with occupation density π (x) = αx -(α+1) , x, α > 0. By the image measure theorem, the density of

x * N +1 (k + 1) d = π -1 (τ N +1
) is obtained as a power-gamma density [START_REF] Janson | On the total external length of the Kingman coalescent[END_REF]

f x * N +1 (k+1) (x) = α N ! x -((N +1)α+1) e -x -α , x > 0.
Next, the conditional density of each x n (k + 1) given

x (N +1) (k + 1) = x is f xn(k+1) x n | x (N +1) (k + 1) = x = π (x n ) π (x) 1 xn>x ,
showing that (with the two terms in the right-hand side term mutually independent), x n (k + 1)

d = x (N +1) (k + 1) X n (k + 1)
where X n (k + 1) is a Pareto(α) distributed rv with density f (x) = αx -(α+1) on (1, ∞) 8 . Putting all this together, we obtained Proposition 12. Independently for each k ≥ 0, with x * N +1 (k + 1) having the power-gamma distribution [START_REF] Janson | On the total external length of the Kingman coalescent[END_REF] and X n (k + 1) being Pareto(α) distributed and with the three right-hand side terms being mutually independent, [START_REF] Karlin | Direct product branching processes and related Markov chains[END_REF] x n (k + 1)

d = x N,α (k) x * N +1 (k + 1) X n (k + 1)
, n = 1, ..., N, indicates how to update multiplicatively the fitness of the n-th individual at generation k+1, when the fitnesses of the previous generation are summarized in x N,α (k) .

We also clearly have an update of the global equivalent fitness x N,α from step k to step k + 1 as:

Corollary 13. With x N,α (k) := N n=1 x n (k) α 1/α
the global equivalent fitness of the whole population at generation k, the following recursion holds [START_REF] Kersting | The asymptotic distribution of the length of Beta-coalescent trees[END_REF] x N,α (k + 1)

d = x N,α (k) x * N +1 (k + 1) N n=1 X n (k + 1) α 1/α
. 8 We used the scaling property of Pareto(α) rvs X on (1, ∞) stating that X | X > a d = aX.

In the latter sum term, each X α n is thus a Pareto(1) distributed rv with density f (x) = x -2 on (1, ∞) and we need to sum N of them independently which is reminiscent of [START_REF] Freund | On the time back to the most recent common ancestor and the external branch length of the Bolthausen-Sznitman coalescent[END_REF].

Large N asymptotics of the α-mean fitness.

Defining x N,α (k) := 1 N N n=1 x α n (k) 1/α
, to be the generalized (Hölder) α-mean of the fitnesses x n , n = 1, ..., N, at generation k, it follows from ( 23) that [START_REF] Kingman | Poisson processes[END_REF] x N,α (k + 1)

d = x N,α (k) x * N +1 (k + 1) N n=1 X n (k + 1) α 1/α
.

By Jensen inequality, the Hölder α-means x N,α (k) are non-decreasing functions of α.

Corollary 14. The β-moments of x N,α at generation k are given by

E x N,α (k) β = x N,α (0) β   E x * β N +1 E   N n=1 X α n β/α     k .
Proof: This follows from the independence of the x * N +1 s and the X α n and their i.i.d. character within each generation k and from [START_REF] Kingman | Poisson processes[END_REF]. ⋄ Proposition 15. For large N , with v N := log log N

1 k log x N,α (k) a.s. → k→∞ 1 α v N .
With

F N (β) := - β α log log N - β α log N (ψ (1 -β/α) -log log N -1)
and f N (a) , a < 0, its Legendre transform, the large deviation regime is given by

1 k log P - 1 k log x N,α (k) → a → k→∞ f N (a) ≤ 0. Proof: Observing that E x * β N +1 = Γ (N + 1 -β/α) /Γ (N + 1) ∼ N →∞
N -β/α and applying [START_REF] Freund | On the time back to the most recent common ancestor and the external branch length of the Bolthausen-Sznitman coalescent[END_REF] giving the moments of a partial sum of N i.i.d. Pareto(1) distributed rvs, it follows that, with

β < α E   N n=1 X α n β/α   ∼ N large (N log N ) β/α 1 + β α log N (ψ (1 -β/α) -log log N -1) .
As a result, we get

E x N,α (k) β x N,α (0) β ∼ N large (log N ) (βk)/α 1 + β α log N (ψ (1 -β/α) -log log N -1) k .
Thus, for large N , with

F N (β) = - β α log log N - β α log N (ψ (1 -β/α) -log log N -1)
defining the concave thermodynamical 'pressure',

- 1 k log E x N,α (k) β → k→∞ F N (β) .
Thus, with a = F ′ N (β) < 0, by the large deviation principle

(26) 1 k log P - 1 k log x N,α (k) → a → k→∞ f N (a) ≤ 0,
where f N (a) = inf β<α (aβ -F N (β)) is the concave Legendre transform of F N , giving the large deviation rate function oflog x N,α (k) /k. In particular, for large N , with v N := log log N

(27) 1 k log x N,α (k) a.s. → k→∞ -F ′ N (0) ∼ 1 α v N ,
gives the limiting right shift of the Hölder α-mean fitness induced by selection effects. ⋄

Remarks:

(i) The limiting right-hand-side term in [START_REF] Lamperti | Continuous state branching processes[END_REF], although increasing very slowly with N , does not stabilize to a limit, in contrast to other similar models [START_REF] Bérard | Brunet-Derrida behavior of branching-selection particle systems on the line[END_REF] of branching with selection where, in each generation, each individual produces only two offspring with randomly shifted fitnesses.

(ii) With α > 0, let f (x) := x α > 0 define some (increasing) output map of the individuals fitnesses x, with x > 0. Defining

f (x) N (k) := 1 N N n=1 f (x n ) (k) = 1 N N n=1
x n (k) α to be the mean output fitness in generation k of the whole population, then, whatever α, (27) is also

1 k log f (x) N (k) a.s. → k→∞ v N ,
interpreting the speed v N itself. This suggests that it is of interest to work not only on the fitnesses x n themselves (and their α-mean x N,α ) but rather on some deformed version of the fitnesses x α n (and their standard mean f (x) N ). Clearly f (x) N itself obeys the recursion

f (x) N (k + 1) d = f (x) N (k) • x * N +1 (k + 1) α • N n=1 X n (k + 1) α ,
deriving again from [START_REF] Karlin | Direct product branching processes and related Markov chains[END_REF] and the definition of the equivalent global fitness x N,α (k).

(iii) Note finally that, given x N,α (k), the cumulative distribution function (cdf) of the fitness x (1) (k + 1) of the fittest individual among the N individuals at generation k + 1 is given by

P xN,α(k) x (1) (k + 1) ≤ x = P π -1 xN,α(k) (τ 1 ) ≤ x = e -(x/xN,α(k)) -α ,
where τ 1 is exp(1) distributed. Thus the fitness of the fittest individual obeys

x (1) (k + 1) d = x N,α (k) Y (k + 1) , where Y (k), k ≥ 0 is a sequence of i.i.d. Fréchet rvs with cdf P (Y ≤ x) = e -x -α . The conditional mean given x N,α (k) of x (1) (k + 1) is x N,α (k) Γ (1 -1/α
) and its median value x N,α (k) (log 2) -1/α . More generally, when N is large and unconditionally, due to the recursion (24) on the x N,α (k)s:

E x (1) (k + 1) β ∼ x N,α (0) β Γ (1 -β/α) e -kFN (β) .

5.2.

Genealogies. Now we turn to the genealogies of this branching process with selection.

The beta(1, 1β) and Bolthausen-Sznitman coalescents.

Recall that π xN,α(k) (x) := -π ′ xN,α(k) (x)
is the occupation density that there would be a point (an offspring) of the PPP at position (with fitness) x at generation k + 1, given a global population state x N,α (k) . Proposition 16. Looking backward in time, upon scaling time using c N ∼ 1/ log N , the genealogy of the branching model with selection is a beta(1, 1β) coalescent, reducing to the Bolthausen-Sznitman coalescent if β = 0.

Proof: Suppose first β = 0. Given there is an offspring at x at generation k + 1, the sampling probability that it would be an offspring of the individual with fitness x n (k) is thus

π xn(k) (x) π xN,α(k) (x) = αx n (k) α x -(α+1) αx N,α (k) α x -(α+1) = x n (k) α x N,α (k) α ,
which is independent of x. Observing from ( 23) and ( 24) that

x n (k) d = x N,α (k -1) x * N +1 (k) X n (k) , n = 1, ..., N, x N,α (k) d = x N,α (k -1) x * N +1 (k) N n=1 X n (k) α 1/α
, this random probability is also

X n (k) α N n=1 X n (k) α ,
where, for each k independently, the X n s are i.i.d. Pareto(α) distributed on (1, ∞) . Because in each generation, parents generate offspring independently and independently of one another, upon averaging, the probability that, at generation k + 1, i individuals share the same common ancestor is independent of k, with

P (N ) i,1 := N n=1 E   X n (k) α N n=1 X n (k) α i   = N E   X 1 (k) α N n=1 X n (k) α i   .
The X n (k)s being i.i.d. Pareto(α) distributed rvs, the X n (k) α s are i.i.d. Pareto(1) distributed rvs and we are thus back to the results of Proposition 6 with β = 0, stating that with c N ∼ 1/ log N,

c -1 N P (N ) i,1 → N →∞ 1 0 u i-2 Λ (du) = Γ (2) Γ (i -1) Γ (i) = 1 i -1 ,
where Λ d ∼beta(1, 1) , uniform. We can proceed similarly to derive the probabilities P (N ) i,j that i individuals have j < i parents, behaving consistently with ( 16) with β = 0. We conclude that, whatever α, in the large N limit, the time-scaled genealogy of the branching model with selection is a Bolthausen-Sznitman coalescent process, obtained while sampling from N Pareto(1) i.i.d. rvs. Would the sampling probabilities P (N ) i,j include a β-size biasing effect on total length Σ N := N n=1 X α n , the genealogy of this branching model with selection would be a full beta(1, 1β) coalescent, provided β < 1. While adopting this sampling point of view to compute the coalescence and merging probabilities, we therefore obtain a limiting genealogical coalescent process which is independent of α. ⋄ Genealogies from the output PPP.

What now if we set that the occupation density that, at generation k + 1, there is an offspring at x descending from some individual with fitness x n (k) at generation k, is instead given by [START_REF] Maruyama | Stochastic problems in population genetics[END_REF] π xn(k) (x) := x n (k) x -2 , while distorting the original occupation intensity π xn(k) (x) = αx n (k) α x -(α+1) ?

Then the occupation density that, at generation k + 1, there is an offspring at x descending from any individual of the whole population would take the form π xN,1(k) (x) := x N,1 (k) x -2 , where x N,1 (k) := N n=1 x n (k) is the cumulative fitness in generation k. If this were to be the case, given there is an offspring at x at time k + 1, the sampling probability that it is an offspring of the individual with fitness x n (k) would be, thanks to [START_REF] Karlin | Direct product branching processes and related Markov chains[END_REF] and [START_REF] Kersting | The asymptotic distribution of the length of Beta-coalescent trees[END_REF] with α = 1 :

(29) π xn(k) (x) π xN,1(k) (x) = x n (k) x -2 x N,1 (k) x -2 = x n (k) x N,1 (k) d = X n (k) N n=1 X n (k)
, again independently of x. This probability now involves a normalized sum of the X n s, which are i.i.d. Pareto(α) distributed and the strategy to compute the merging probabilities of the ancestral process will be modified.

Under this hypothesis indeed, the probability that, at generation k +1, i individuals share the same common ancestor reads

P (N ) i,1 := N n=1 E   X n (k) N n=1 X n (k) i   = N E   X 1 (k) N n=1 X n (k) i   =: N E S 1 (k) i ,
where S 1 is the normalized segment size now obtained from N i.i.d. Pareto(α) distributed rvs, normalized by their sum Σ N (k) := N n=1 X n (k). We can proceed similarly to derive the probabilities P (N ) i,j that i individuals have j < i parents and we are back to the studies of Sections 2 -4. And we can as well β-size-bias these sampling probabilities on the total lengths Σ N . We call this sampling procedure the distorted sampling procedure. -a continuous-time beta(2α, αβ) coalescent if α ∈ (1, 2) and β < α (upon scaling time with c N ∝ N -(α-1) ).

-a continuous-time beta(1, 1β) coalescent if α = 1 and β < 1 (upon scaling time with c N ∝ 1/ log N ).

-a discrete-time Poisson-Dirichlet(α, -β) coalescent if α ∈ (0, 1) and β < α.

Proof: It remains to interpret the distorted size-biased sampling procedure which is proposed to compute the merging probabilities of the ancestral process: Assume that the fitness dependent PPP describing the descent of an individual with fitness x n (k) is now the output image of the original one, given by the canonical application f xn(k) (x) = x n (k) (x/x n (k)) α , f xn(k) : R + → R + . Its gives rise to a new PPP with distorted intensity π f xn (k) (dx) = x n (k) x -2 dx, the image measure of π xn(k) (dx) = π xn(k) (x) dx = αx n (k) α x -(α+1) dx by f xn(k) (as a result of the classical Campbell formula for PPPs, see [START_REF] Kingman | Poisson processes[END_REF] p. 28). This was the starting point in [START_REF] Maruyama | Stochastic problems in population genetics[END_REF]. The skewed computation in [START_REF] Möhle | On the number of segregating sites for populations with large family sizes[END_REF] of the probability that some offspring is descending from one individual with fitness x n (k) is thus based not on the original PPP attached to x n (k) but rather on a deformed version of it through f xn(k) . We note that f xn(k) (x), as a function of the two arguments (x n (k) , x) is homogeneous with If α ∈ [0, 1) , β < α, the obtained large N genealogical coalescent will coincide with the discrete-time-k Poisson-Dirichlet coalescent with parameters α and -β.

Only when α = 1 do we get as in [START_REF] Brunet | Genealogies in simple models of evolution[END_REF] (upon scaling time logarithmically with N ) the Bolthausen-Sznitman coalescent (β = 0) or more generally the beta(1, 1β) coalescent, provided β < 1. ⋄

4

  Here, because two processes are involved, the symbol d = means convergence of all the finitedimensional distributions of x (N) [t/c N ] , t ≥ 0 to the ones of xt, t ≥ 0.

  with a, b > 0 and with B (a, b) the beta function: we get the beta(a, b) coalescents.

Proposition 7 .

 7 When α > 2, for any value of β, with µ := α α-1 > 0 and ρ := α α-2 > 0, upon scaling time using a linear effective population size N e = c -1 N = N µ 2 /ρ, the limiting process x t d = lim N →∞ x (N ) [t/cN ] exists and is the continuous-timet Kingman coalescent.

Proposition 9 .

 9 similar scaling behavior for c N was recently obtained in Theorem 2.4 of [20], dealing with coalescents arising from compound Poisson discrete reproduction models, in the critical case. When α = 2, upon scaling time using an effective population size N e = c -1 N = (2N ) / log N, the limiting process x t d = lim N →∞ x (N ) [t/cN ] is the continuous-time-t Kingman coalescent.

<sN <s and so f τ 1

 1 ,...,τN (s 1 , .., s N ) = 1 0<s1<...<sN ∞ sN dse -s = e -sN 1 0<s1<...<sN .

Proposition 17 .

 17 Looking backward in time, using a distorted size-biased sampling procedure, the genealogy of the branching model with selection is-a continuous-time Kingman coalescent if α ≥ 2 (upon scaling time with c N ∝ 1/N if α > 2 or c N ∝ log N/N if α = 2).

f

  λ a xn(k) λ b x = λ a(1-α)+bα f xn(k) (x), λ > 0, leading obviously, if λ = x n (k) -1/a , b/a = 1 and f (x) := f 1 (x) = x α to: f xn(k) (x) = x n (k) f (x/x n (k)). The function f (x) = x α is the output fitness function introduced in Subsection 5.1, Remark(ii).Using this distorted size-biased sampling procedure therefore, following the introductory arguments, the full class of the Pareto-coalescents (described in Sections 2 -4) are obtained.Suppose for instance α ∈ (1, 2), β = 0. Based on the previous computations of Sections 2 -4, we conclude that the large N distorted genealogy of the branching model with selection coincides (upon scaling time correspondingly: k → [t/c N ]) with a beta(2α, α) coalescent process, obtained while sampling from N Pareto(α) i.i.d. rvs, normalized by their sum. Would this probability involve a β-size biasing effect, the genealogy of this branching model with selection is identified to a beta(2α, αβ) coalescent, β < α.

E Σ β N = N E S i 1 = N Γ (N θ) Γ (θ)

c N = N α E Σ β-α N -1 E Σ β N Γ (2α) Γ (αβ) Γ (2β) ∼ α Γ (1β) 1 -β α Γ (1 -(βα)) Γ (1α) Γ (2α) Γ (αβ) Γ (2β)

(N ) i,j → P i,j where (1 ≤ j ≤ i)(10) P i,j = i! j! α j-1 Γ (1β) Γ (1β/α) Γ (jβ/α) Γ (iβ) * i1+...+ij=i j l=1

* i1+...+ij =i i i 1 ...i j φ j (i 1 , .., i j ) ,where the φ j (i 1 , .., i j )s define the probabilities of a (i 1 , .., i j ) -merger (i 1 + .. + i j = i). These φ j (i 1 , .., i j ) , which can be read from[START_REF] Dhersin | On the length of an external branch in the Beta-coalescent[END_REF], may be written under the alternative formφ j (i 1 , .., i j ) : = α j-1 Γ (1β) Γ (1β/α) Γ (jβ/α) Γ (iβ) j l=1 Γ (i lα) Γ (1α) = c j,α,β Γ (αjβ) Γ (iβ) j l=1 Γ (i lα) ,where c j,α,β :=

Acknowledgments:

The author acknowledges partial support from the ANR Modélisation Aléatoire en Écologie, Génétique et Évolution (ANR-Manège-09-BLAN-0215 project) and from the labex MME-DII (Modèles Mathématiques et Économiques de la Dynamique, de l' Incertitude et des Interactions). The author is also indebted to his referees for pointing out some errors in an earlier version of the draft and for encouraging him to write down a more concise and complete version.

Laboratoire de Physique Théorique et Modélisation, CNRS-UMR 8089 et Université de Cergy-Pontoise, 2 Avenue Adolphe Chauvin, F-95302, Cergy-Pontoise, France, E-mail: Thierry.Huillet@u-cergy.fr


Valeria Loscrì
email: vloscri@deis.unical.it

Enrico Natalizio

Valentina Mannara

Gianluca Aloi

Valeria Loscrí

A novel Communication Technique for Nanobots based on acoustic signals

Keywords: Acoustic Communication, Nanobots, In-vivo Applications

In this work we present the simulation of a swarm of nanobots that behave in a distributed fashion and communicate through vibrations, permitting a decentralized control to treat endogenous diseases of the brain. Each nanobot is able to recognize a cancer cell, eliminate it and announces through a communication based on acoustic signals the presence of the cancer to the other nanobots. We assume that our nanodevices vibrate and these vibrations cause acoustic waves that propagate into the brain with some intensity that we evaluated by taking into account the specific physical factors of the context, the nano-metric nature of the vibrant devices and the characteristic of the fluid where the devices are immersed. An important aspect of our approach is related to the communication based on vibrations. This choice is related to the application context where is not advisable either to use indiscriminate chemical substances or electromagnetic waves. Whereas, ultrasonic waves are used in the most frequent diagnostic techniques and the use of this kind of techniques should not have negative collateral effects. Specifically, we propose an approach based on bees' behavior in order to allow our devices to communicate, coordinate and reach the common objective to destroy the cancerous tissues. In order to evaluate the effectiveness of our technique, we compared it with other techniques known in literature and simulation results showed the effectiveness of our technique both in terms of achievement of the objective, that is the destruction of the cancerous cells, and velocity of destruction.

Introduction

Nanotechnologies are a new approach based on comprehension and deep knowledge of the properties of the matter at the nanoscale level: a nanometer corresponds to the length of a small molecule. At nanoscale, behaviors and characteristics of the matter drastically change, making it necessary a synergy among many different disciplines. A great percentage of the nanosytems world is represented by technologies derived from micro-electronics, pushed up to nanoscale, in order to obtain electronic, optical, fluidic, integrated mechanic functions to be applied to different fields ranging from microelectronic, telecommunications and sensor, to environmental and bio-medical. The latter field could be literally revolutionized from the potential applications of nano-devices, both in the diagnostic and pharmaceutic fields [START_REF] Freitas | Nanomedicine, Volume I: Basic Capabilities[END_REF]. One of the most interesting applications is the controlled release of drugs over time and exactly localized in cells or organs that need it, drastically reducing the side effects. Nanotechnology is already used in the field of diagnostics through the use of synthetic tracer molecules for investigating biological processes in a non-invasive fashion. In the cancer fight, as for the treatment of some diseases related to the cardiovascular system more and more often you hear speak about nano-therapy. Among the first applications to be postulated in the early '90s and also among the most fascinating of nanomedicine is the idea of using nanorobots [START_REF] Atakan | Body Area NanoNetworks with Molecular Communications in Nanomedicine[END_REF]. A device of a few nanometers (nanorobots will typically be 0.5 to 3 microns large with 1-100 nm parts), could be introduced into the body without causing injury and, if equipped with sensors that transmit precise images could facilitate the early diagnosis of cancer and carry drug to the target or to perform other tasks that would otherwise require invasive surgery [START_REF] Cavalcanti | Nanorobot Hardware Architecture for Medical Defense[END_REF], [START_REF] Amato | Early-Stage Diagnosis of Endogenous Diseases by Swarms of Nanobots: an Applicative Scenario[END_REF]. In order to prevent attacks from immune system a nanorobot in vivo should be characterized with a smooth and flawless diamond exterior, because this prevents Ieukocytes activities since the exterior is chemically inert and have low bioactivity [START_REF] Davis | Nanorobots in Brain Tumor[END_REF].

A nanorobot is a system able to modify the surrounding in a controlled and predictable fashion, with size at the molecular or even atomic scale. In practice, a nanorobot is either a passive or active structure able to detect, signal and elaborate information. The very limited size of the devices implies limited capabilities and reduced computation resources, therefore it is necessary to make devices collaborate by applying design techniques such as Swarm Intelligence in order to realize complex systems through the interaction and the cooperation of very simple agents [START_REF] Dorigo | Optimization, Learning and Natural Algorithms[END_REF].

This work is inspired by the paper of Lewis Anthony and Bekey George developed in 1992 at Southern California University and entitled "The Behavioral Self-Organization of Nanorobots Using Local Rules" [START_REF] Lewis | The Behavioral Self-organization of Nanorobots using Local Rules[END_REF], where the objective was the evaluation of a swarm of nanobots, organized without any centralized unit, only trough simple local rules, to destroy cancer cells inside the brain. Thanks to the technological development, the assumptions made in [START_REF] Lewis | The Behavioral Self-organization of Nanorobots using Local Rules[END_REF] are plausible and it is possible to design molecular nanobots constituted from two polymers able to deactivate the production of a protein that causes the death of the cells attached [START_REF] Davis | MoEvidence of RNAi in humans from systematically administered siRNA via targeted nanoparticle[END_REF]. Lewis's work has an early work for the assumption of using robots for in vivo applications through a swarm technique. In [START_REF] Lewis | The Behavioral Self-organization of Nanorobots using Local Rules[END_REF] and in other next works [5] [6], the communications are performed by chemical mechanisms, that is, the report of either some events or object is realized through the release of chemical substances that diffuse by attracting other nanorobots through the gradient associated to the signal intensity. Specifically, Lewis's strategy is based on three chemical substances that attract the swarm of nanobots in a gradual fashion to achieve the total elimination of the diseased cells. An early part of the work is devoted to the tuning of the characteristic parameters of the algorithm where objective difficulties emerged, first of all the extremely reduced dimensions of the devices, their limited computational capabilities and the structural space needed to install sensors and interfaces to capture all the chemical substances involved. For these reasons, we considered a different approach inspired from bees' behaviors that use vibration to communicate the distance and the position of food sources to the rest of the colony. Our choice for this specific kind of communication has been lead by the simplicity of the components required from the devices, and the fact that many diagnostic instruments are already based on microwaves, which are considered not dangerous for the human. A useful summary of communication paradigms to interconnect nanobots in a body area network is given in [START_REF] Akyildiz | Nanoneworks: A New Communication Paradigm[END_REF]. Preliminary results of Bee's algorithm are presented in [START_REF] Loscrí | Efficient Acoustic Communication Techniques for Nanobots[END_REF]. The rest of the paper is structured as follows. In Section II we give some details about Lewis's approach. In Section III we describe our bees's approach. In Section IV we give the details about the simulation environment and the results. Finally, we conclude the paper in Section V and we give some directions for future works in Section VI.

Lewis Algorithm

When Lewis and Bekey presented their work in 1992, the Swarm Intelligence concept had not yet been fully defined. In fact, authors refer to a μ -colony as a set of robots that coordinate activities autonomously to perform some specific tasks. The colony could be constituted by hundreds or thousands of μ -robots, and each of them has simple computational capabilities like a colony in nature. As we already outlined, the technique they proposed was based on chemical communication. Each μ -robot is able to mark its surrounding through chemical substances, recognizes the different chemical signals and follows the different gradients until it reaches the cancer cells. The colony is injected close to the cancer. Once injected into the body, the μ -robots move randomly until they reach cancer cells. After the first contact the μ -robot emits a substance in its surrounding called CHEM-1. This substance is absorbed by the body after a certain time. A certain percentage of μ -robots differentiate in guidepost, stop and start to secrete substances that permit the transmission of the signal over long distances. These substances are CHEM-2 and CHEM-3 and are used as repeaters. The number of μ-robots that differentiate in guideposts determines the efficiency, which is the convergence of μ -robots to defeat the tumor. The number of guideposts is very important. In fact, if too many nanobots differentiate it is difficult to reach and destroy the tumor. On the other hand, if the number of guideposts is too low, the colony is not able to complete the task. For this reason, authors considered a differentiation probability p = 0.01 and related the total number of guideposts to the total number of μ -robots n and the current time t, pnt. In what follows we show the pseudo-code of the Lewis-Bekey approach: The algorithm's efficiency is strictly related to some additional parameters, like the threshold value of CHEM-1, Θ. The environment considered is the cerebral cortex that is represented in the simulator NSl [START_REF] Weitzenfeld | The Neural Simulation Language[END_REF] as a cells array. The modeled chemical communication takes into account the generation, the diffusion of the substances immersed in a fluid and the absorption factor. A single cell is the measure unit of the space and only one nanobot at a time can occupy a single cell. The cancer and the colony are positioned inside the grid. Each robot is able to move step by step in one of the 8 adjacent cells. Once a μ -robot reaches a cancer cell it is able to perceive the tumor and to signal its presence to other robots. In our simulation, we adopted two different convergence criteria. The first one is represented by the total elimination of the cancer cells and the other is a maximum number of iterations (1000), that is not an actual convergence criterion but we need to stop the simulation when the algorithm fails to eliminate all the cells. Even if the algorithm is strictly related to the setting of different parameters, the aim of the authors was mainly to show how behavioral algorithm is able to make the colony reach a common objective. It is worth to note that, the use of three chemical signals implies at least three different interfaces to capture the different substances. Moreover, a single nanobot needs to physically transport three different substances. By considering the reduced dimensions it seems very difficult to imagine how this algorithm could be realized.

Bees' Approach

As we already outlined above, the acoustic communication we propose in this work is inspired from bees exploration technique for food searching. In practice, our nanobots borrow from the bees the capabilities to communicate through vibration (waggle dance for the bees). The self-organization of the bees is based on very simple rules related to the behavior of each individual. Moreover, the concept of swarm applied in the context of in vivo application has several advantages compared to isolated nanobot. For example, acoustical nanobots could form in vivo communication networks that could transfer data across much larger distance than possible with direct transmission by considering the attenuation at high frequencies [START_REF] Stojanovic | On the relationship between capacity and distance in an underwater acoustic communication channel[END_REF]. Generally, nanorobots can improve their performance and they are able to accomplish complex tasks, by coordinating their actions in a decentrilezed fashion. A similar approach is considered in [START_REF] Hogg | Controlling tiny multi-scale robots for nerve repair[END_REF]. We refer to this bees' inspired technique as NanoBee and is supported by the possibility of using acoustic waves as transmission means in communications in vivo without specific risks associated to. We exploited the analysis made in [START_REF] Freitas | Nanomedicine, Volume I: Basic Capabilities[END_REF] to characterize our simulation model and tune the parameters. The vibrations associated to the devices generate acoustic waves that propagate in an elastic medium and cause pressure variations and movement of the particles that compose the medium and that can be perceived and detected from an acoustic detector. This perturbation, that carries both the information and the energy, propagates while every particle, also in the case of a fluid, remains nearby its original position. In practice, there are local vibrations (compression and rarefaction) of the particles and in the case of fluid, where there can not be cutting efforts, the vibrations are parallel to the propagation direction of the wave, i.e. longitudinal waves. The dimension or magnitude of sound can be indifferently expressed as sound power, sound intensity or sound pressure. Sound power is the total amount of acoustic energy emitted from a sound source and is measured in watt. Sound intensity is the ratio of the power of a sound wave and the crossed superficial area, it is usually measured in watt/meter 2 . Sound pressure is the value of the pressure variation of a corps in a generic point inside the sound field and is measured in newton/meter 2 . It is worth to note that for a given sound, sound power is constant while both intensity and pressure depend on the specific measure conditions. The simplest case is sound or noise that propagate freely without any obstacle. With this assumption and when the medium considered is non-dissipative, intensity (I), power (P) and pressure (p) are correlated as it follows:

I = W 2πr 2 = p 2 ρc ( 1
)
where ρ is the density of the fluid and c is the light speed. It is worth to note that both intensity and pressure decrease with the square of the distance (r). Based on these considerations we assume the possibility to design a technique based on acoustic communication by modeling the propagation of signals with sound intensity. NanoBee technique considers nanobots as point sources, since their size is small compared to the distance from the receiver. Each nanobot is able to signal the presence of a cancer cell by stopping its movement and starting vibrating, mimicking the bees behavior, in order to transmit a sound signal to alert other agents (nanobots). Every nanobot that receives the signal moves towards the gradient of the intensity received. Vibrating time is limited since it depends on the force associated to each device to accomplish the movement. In order to simulate the spatial variation of intensity sound, we assume three different probabilistic intervals, that simulate both the spatial propagation and the attenuation of the sound. Specifically, we choose higher probabilities values close to the nanobot that is dancing and probabilities values smaller when distances are greater by simulating in this way the attenuation with the square of the distance. A more realistic version of the algorithm must take into account the temporal attenuation of a sound signal which introduces a decrease of both the intensity and the distance reached by the signal. The first version, where we do not take into account the temporal attenuation of the signal is referred as NanoBee ON-OFF, where ON and OFF indicate whether the device is "ON" (it vibrates with the maximum power) otherwise it is "OFF" (and does not vibrate at all). The second version of our algorithm is referred as NanoBeeEvan and is related with a temporal attenuation of the sound signal. From a computational point of view, our devices are very simple since the algorithm only requires that each device has capabilities of:

-recognition of a cancer cell; -destruction of a cancer cell; -emission of vibrations to signal the position of the cancer; -detection of acoustic waves.

In what follows we give the pseudo-code of the NanoBeeAlgorithm: Algorithm 2 NanoBee Repeat for each "active" nanobot i: pick up any sound signal in its surrounding; IF there are not any signal/cancer cells search randomly; IF discovered a cancer cell THEN eliminate the cell and starts to dancing; IF a signal has been received THEN moves towards the higher intensity of the signal; Until there is an "active" nanobot A nanobot is active when it has enough energy, above a certain predetermined threshold and it has not discovered any cancer cell yet.

Simulation Results

We simulate the systems previously described with Neural Simulation Language (NSL), [START_REF] Weitzenfeld | The Neural Simulation Language[END_REF]. NSL possesses many features that facilitated the simulation development, including graphics capabilities, and language constructs for handling layers of grid object. In order to evaluate our techniques in different situations we consider the possibility to have different initial configurations of the cancer. In fact, based on the specific geometry of the processes, different cancer forms characterize different cancer types. Specifically, we consider the possibility to have metastases detached from the original tissue. Usually, metastases are a group of cells and only rarely they are isolated, but for completeness we considered the three different scenarios: "aggregated", "scattered groups" and "isolated".

Different Initial Configuration

The space where nanobots move and act is a grid of cells modeling the vertebral cortex. In Figure 1 we show the initial configuration when the "aggregated" case is considered, in Figure 2 the initial configuration of the cancer is with "groups" of cells and in Figure 3 we show the initial configuration of the cancer with "isolated" cells. Specifically, the number of cancer cells is 42 in every configuration and the cancer cells are represented in blue, while the nanobots, that are 289 in every configuration, (as in [START_REF] Lewis | The Behavioral Self-organization of Nanorobots using Local Rules[END_REF]), are represented in black.

Convergence Time

The input parameters we consider in order to evaluate the performance of the various techniques are summarized in table 1. The algorithms considered in this work are based on the concept of nanobots swarm, which allows us to consider a reduced dimension of devices with available power smaller than that necessary if only a device is considered as observed in [START_REF] Hog | Acoustic Communication for Medical Nanorobots[END_REF]. Frequencies ranging from 1 to 5 MHz are typically used for the head [START_REF] Hog | Acoustic Communication for Medical Nanorobots[END_REF], thus, we choose the value of frequency equal to 1 MHz. Moreover, in our bees' approach, communication activity is only related to the achievement of the target. In practice, nanobots start to vibrate when they discover a tumor and not when they move. We apply the method of Independent Replications with a confidence interval of 95%.

In Fig. 4 we show the convergence time of each considered algorithm for the "aggregated" case. The first important observation is related to Lewis approach (referred in the diagram as NanoLewis) that is not able to totally destruct the cancer cells even if we consider the maximum number of nanobots. Concerning Lewis'approach we had to set many parameters and after conducting a sensitivity analysis we chose the set of parameters that guarantee the best performance in terms of convergence. During the simulation we observed as the combination of the three chemical substances creates a kind of "barrier effect". In practice, after the first nanobots enters in contact with cancer cells and a part of the devices differentiate and start working as a kind of relay by sending CHEM-1. At this point, it is very difficult for the other nanobots to attack the most internal cancer cells. When the stations differentiate like "guidepost" and start send CHEM-2 if they already entered in contact with CHEM-1 or CHEM-3 if they were touched from CHEM-2, the situation worsens. Bees' approaches behave similarly, with the difference that NanoBee ON-OFF is able to eliminate all the cancer cells with less nanobots than NanoBeeEvan. NanoBee ON-OFF is more effective in this case cause of the signal propagation that reaches farther during the time and is able to attract more nanobots. NanoBeeEvan consideres the attenuation of the signal power during the time. Hence, after a certain period the signal only reaches a shorter distance.

In Fig. 5 we can observe how the NanoLewis is not yet able to destruct all the cancer cells, in the "scattered group" algorithm. In this case we observe a little improvement of the NanoBeeEvan. This latter, through the "evanescence" of the signal in the time, implies more randomness of movement of the nanobots. Hence, when tumor cells are scattered in the region of interest, this randomness is useful to reach the various points of interest more quickly.

This reasoning is confirmed in Fig. 6 where NanoBeeEvan overcomes the NanoBee ON-OFF, for the case of "isolated" cells.

It is worth to note how both the nano-bee approaches require less nanobots to reach the goal to eliminate all the cancer cells in respect of those needed by Lewis algorithm. In our opinion the analysis of the performance in terms of both effectiveness, that is the capability to destruct all the cancer cells, and the rapidity to do it by varying the number of nanobots, is very important because it is very important to use the minimum number of nanobots in a similar context. Furthermore, it is very important to try to reduce both the chemical and the acoustic messages as much as possible.

Amount of Eliminated Cancer Cells

In this section we present another type of analysis related to the variation of the number of nanobots but also related to the capability to eliminate the cancer cells. The simulation parameters are always those summarized in the Table 1.

In Fig. 7 we can observe how the number of cancer cells eliminated by the NanoLewis technique is always smaller than the other two techniques. As already outlined before, we observe that the algorithm is effective when the first nanobot attaches the tumor, but the barrier created through the chemical substances does not allow to penetrate inside the body of the cancer. In fact, even if the number of nanobots increases, above all, in the first configuration (where we have only a single tumor), the effect does not change. This shows that the scarce efficiency is not related to the small number of nanobots. On the contrary, with our bees techniques we have that an increasing number of nanobots corresponds to an increased capability to attack the cancer. In fact, in the "aggregated" case, we pass from around 18 cells to 42 cells. In Fig. 8 we can observe how the barrier effect of the NanoLewis approach is mitigated in the "scattered groups" configuration. Moreover, for this specific configuration the bees approaches behave very similarly. In Fig. 9 we can notice how the effects of the randomness improves the behavior of the NanoBeeEvan that is very close in terms of performance to the NanoBee ON-OFF. From the analysis conducted is clear that the algorithms considered with the associated parameters in terms of size of devices, available power, etc. are feasible for nanobots operating in the specific environment as we defined (i.e. the brain). When a different environment is considered, the characteristic parameters have to be re-defined and the effectiveness of the technique needs to be proven. Invivo applications need a very accurate choice of available parameters, that can change based on the specific application they are considered for.

Conclusions

In this work we considered and evaluated the possibility of using and exploit a swarm of nanobots in medical applications for the treatment of endogenous diseases and through this preliminary study we showed how this is possible. Specifically, we focused on the communication and organization aspects, necessary for the control of a single nanobot through very simple rules and based only on local information exchange. We took inspiration from a previous work where the behavior of the swarms of nanobots was defined through chemical signals able to attract the devices. After an accurate evaluation of the types of communication physically feasible for in vivo communications and based on the bees behavior in the supply phase, we proposed an algorithm based on acoustic signals, NanoBee. This choice is supported from the common use of sound waves for medical applications. By assuming the possibility for a nanobot to generate a vibration when it picks up a cancer cellular tissue, we evaluated this approach and compared the algorithms via simulation by a neural simulation tool, NSL. We showed that nano-bee techniques are more effective to face cancer with different initial shape than Lewis'approach, by eliminating all the malignant cells. Communications are more efficient and nodes only need to perceive the "sound" by moving a step at each time. We showed how the proposed techniques based on the bees behaviors are effective also when the number of nanobots is smaller than that considered in the Lewis and Bekey technique.

Future Works

Concerning the future works there are many aspects that are worth to be considered in a deep way. The first one is concerning the chemical communication approach considered from Lewis and Bekey. In fact, during the simulation we were able to observe that the "barrier effect" worsens when CHEM-2 and CHEM-3 are activated. Maybe, an approach with only a single chemical substance would be more effective, and in terms of nanobot design it would be more feasible. An other important aspect is the validation of the probabilistic approach we adopted to simulate the attenuation of the signal with the square of the distance. It would be interesting to evaluate a real model of the brain and compare it with our probabilistic approach. Another important aspect to be considered is related to a more realistic simulation environment, where "obstacles" (i.e. others cells) are explicitly taken into account and the interaction of the signals is analyzed.

Algorithm 1 1 IF

 11 Lewis and Bekey AlgorithmIF There is no chemical Markers and No Tumor THEN Do a random Walk IF A Tumor is detected THEN Destroy the cell, Differentiate and Broadcast CHEM-1; w/prob The Mag. of CHEM-1 is greater than Θ THEN Do a random Walk IF CHEM-1 is detected THEN Move up the gradient of CHEM-1; w/prob p or differentiate and Broadcast CHEM-2; w/prob(1-p) IF CHEM-2 is detected THEN Move up the gradient of CHEM-2; w/prob p or Differentiate and broadcast CHEM-3; w/prob(1-p) IF CHEM-3 is detected THEN Move up the gradient of CHEM-2; w/prob 1

Fig. 1 .

 1 Fig. 1. Initial Configuration of the "aggregated" case.

Fig. 2 .

 2 Fig. 2. Initial Configuration of the "scattered groups" case.

Fig. 3 .

 3 Fig. 3. Initial Configuration of the "isolated" case.

Fig. 4 .

 4 Fig. 4. Convergence Time in the case of "aggregated" configuration.

Fig. 5 .

 5 Fig. 5. Convergence Time in the case of "scattered groups" configuration.

Fig. 6 .

 6 Fig. 6. Convergence Time in the case of "isolated" configuration.

Fig. 7 .

 7 Fig. 7. Number of eliminated cancer cells in the case of "aggregated" configuration.

Fig. 8 .

 8 Fig. 8. Number of eliminated cancer cells in the case of "scattered groups" configuration.

Fig. 9 .

 9 Fig. 9. Number of eliminated cancer cells in the case of "isolated" configuration.

Table 1 .

 1 Simulation Parameters

	Searching Space	50x50 cells
	N Nanobots	289
	Cancer size	42 cells
	Size of nanobots	0.5 μm
	Frequency	10 MHz
	Power	0.5 pW
	Simulation Time	1000 runs
	Prob Values	[1-0.7],[0.69-0.3],[0.29,0]
	Confidential Interval	95%