Asymptotics of the allele frequency spectrum associated with the Bolthausen-Sznitman coalescent

We work in the context of the infinitely many alleles model. The allelic partition associated with a coalescent process started from n individuals is obtained by placing mutations along the skeleton of the coalescent tree; for each individual, we trace back to the most recent mutation affecting it and group together individuals whose most recent mutations are the same. The number of blocks of each of the different possible sizes in this partition is the allele frequency spectrum. The celebrated Ewens sampling formula gives precise probabilities for the allele frequency spectrum associated with Kingman's coalescent. This (and the degenerate star-shaped coalescent) are the only Lambda coalescents for which explicit probabilities are known, although they are known to satisfy a recursion due to Moehle. Recently, Berestycki, Berestycki and Schweinsberg have proved asymptotic results for the allele frequency spectra of the Beta(2-alpha,alpha) coalescents with alpha in (1,2). In this paper, we prove full asymptotics for the case of the Bolthausen-Sznitman coalescent.


Exchangeable random partitions
In recent years, the topic of exchangeable random partitions has received a lot of attention (see Pitman [35] for a lucid introduction). A random partition of N is said to be exchangeable if, for any permutation σ : N → N such that σ(i) = i for all i sufficiently large, we have that the distribution of the partition is unaffected by the application of σ. It was proved by Kingman [28; 29] that if the partition has blocks (B i , i ≥ 1) listed in increasing order of least elements then the asymptotic frequencies, exist almost surely. Let (f ↓ i ) i≥1 be the collection of asymptotic frequencies ranked in decreasing order. Then we can view (f ↓ i ) i≥1 as a partition of [0, 1] into intervals of decreasing length. In general, since it is possible that i≥1 f ↓ i < 1, there will also be a distinguished interval of length 1 − i≥1 f ↓ i . Consider now the following paintbox process, which creates a random partition of N starting from the frequencies. Take independent uniform random variables U 1 , U 2 , . . . on [0, 1]. If U i and U j land in the same non-distinguished interval of the partition then assign i and j to be in the same block. If U i lands in the distinguished interval, assign i to a singleton block. The partition we create in this way is exchangeable and has the same distribution as the partition with which we began. This procedure can also be thought of in terms of a classical balls-in-boxes problem with infinitely many unlabelled boxes, see in particular Karlin [27] and Gnedin, Hansen and Pitman [18].
There are several natural questions that we may ask about an exchangeable random partition restricted to the first n integers (or, equivalently, about the partition formed by the first n uniform random variables in the paintbox process). How many blocks does this partition have? How many blocks does it have of size exactly k, for 1 ≤ k ≤ n? Even in the absence of precise distributional information for finite n, can we obtain n → ∞ limits for these quantities, in an appropriate sense? These questions have been studied for various classes of exchangeable random partitions and random compositions, see in particular [1; 17; 20; 21; 23; 24].

Coalescent process and allelic partitions
In this paper, we study a particular exchangeable random partition which derives from a coalescent process. The origins of this partition lie in population genetics and we will now describe how it arises and give a brief review of the relevant literature. For large populations, genealogies are often modelled using Kingman's coalescent [30]. This is a Markov process taking values in the space of partitions of N (or [n] def = {1, 2, . . . , n}), such that the partition becomes coarser and coarser with time. Whenever the current state has b blocks, any pair of them coalesces at rate 1, independently of the other blocks and irrespective of the block sizes. We start with a sample of genetic material from n individuals. Here, n is taken to be small compared to the total underlying population size. We imagine tracing the genealogy of the sample backwards in time from the present. Then the blocks of the coalescent process at time t correspond to the groups of individuals having the same ancestor time t ago (where time is measured in units of the total underlying effective population size). See Ewens [16] or Durrett [14] for full introductions to this subject. In the population genetics setting, it is natural to introduce the concept of mutation into this model. One of the most celebrated results in this area is the Ewens Sampling Formula, which was proved by Ewens [15] in 1972. It concerns the infinitely many alleles model, in which every mutation gives rise to a completely new type. It says that if we take a sample of n genes subject to neutral mutation (that is, mutation which does not confer a selective advantage) which occurs at rate ρ for each individual, then the probability q(m 1 , m 2 , . . .) that there are m j types which occur exactly j times is given by where θ = 2ρ, (θ) n↑ = θ(θ + 1) · · · (θ + n − 1) and we must have j≥1 jm j = n. Another way of expressing this (due to Kingman [28]) is to picture the coalescent tree associated with Kingman's coalescent and place mutations along the length of the skeleton as a Poisson process of intensity ρ = θ/2. For each individual, trace backwards in time (i.e. forwards in coalescent time) to the most recent mutation. Group together those individuals whose most recent mutations are the same; this gives the allelic partition. Then m j is the number of blocks in the allelic partition containing exactly j individuals.
It is natural to extend these ideas to more general coalescent processes. See Figure 1 for an example of a general coalescent tree and its allelic partition. The Λ-coalescents are a class of Markovian coalescent processes which were introduced by Pitman [34] and Sagitov [37]. Like Kingman's coalescent, they take as their state-space the set of partitions of [n] (or, indeed, of the whole set of natural numbers). Their evolution is such that only one block is formed in any coalescence event and rates of coalescence depend only on the number of blocks present and not on their sizes. Take Λ to be a finite measure on [0, 1]. In order to give a formal description of the Λ-coalescent, it is sufficient to give its jump rates. Whenever there are b blocks present, any particular k of them coalesce at rate Note that, in contrast to Kingman's coalescent, here we allow multiple collisions; that is, we allow more than two blocks to join together. However, we do not allow more than one group of blocks to coalesce at once. Kingman's coalescent is the case Λ(dx) = δ 0 (dx), where unit mass is placed at 0. The case Λ(dx) = dx, called the Bolthausen-Sznitman coalescent, was introduced by Bolthausen and Sznitman [7] in the context of spin glasses. It has many nice properties and appears to be more tractable than most Λ-coalescents. For example, its marginal distributions are known explicitly [34]. It has been studied in some detail: see, for example, Pitman [34], Bertoin and Le Gall [5], Basdevant [2] and Goldschmidt and Martin [25].
Another sub-class of the Λ-coalescents which has recently been particularly studied is the Beta coalescents, so-called because Λ here is a beta density: for some α ∈ (0, 2). (The α = 1 case is the Bolthausen-Sznitman coalescent and, in some sense, α = 2 corresponds to Kingman's coalescent.) See Birkner et al [6] for a representation in terms of continuous-state branching processes when α ∈ (0, 2).
If we suppose that, instead of Kingman's coalescent, the genealogy of the population evolves according to a general Λ-coalescent then, except in the special case of the degenerate star-shaped coalescent (where Λ(dx) = δ 1 (dx)), there is no known explicit expression for the probability q(m 1 , m 2 , . . .) of having m j blocks in the allelic partition of size j. However, Möhle [31] has shown that the probabilities q must satisfy the following recursion: where λ n = n k=2 ( n k ) λ n,k , ρ = θ/2, m = (m 1 , m 2 , . . .) and e i is the vector with a 1 in the ith co-ordinate and 0 in all the rest. He has also shown [33] that, except in the cases of the starshaped coalescent and Kingman's coalescent, the allelic partition is not regenerative in the sense of Gnedin and Pitman [19]. Dong, Gnedin and Pitman [12] have studied various properties of the allelic partition of a general Λ-coalescent. In particular, they view the allelic partition as the final partition of a coalescent process with freeze (see Section 2 where we use this formalism) and also give an alternative description of q as the stationary distribution of a certain discrete-time Markov chain.
Consider again the Beta coalescents. Suppose that we start the coalescent process from the partition of [n] into singletons. Let N k (n) be the number of blocks of size k, for k ≥ 1, and let N (n) = n k=1 N k (n). From the biological point of view, N k (n) is the number of types which appear exactly k times in a sample of size n and N (n) is the total number of types in the sample. The complete allele frequency spectrum is the vector In the case of α ∈ (1, 2), Berestycki, Berestycki and Schweinsberg [3; 4] have proved that and, for k ≥ 1, that as n → ∞.
The corresponding convergence results for Kingman's coalescent can be derived from the Ewens sampling formula: without rescaling, we have where Z 1 , Z 2 , . . . are independent Poisson random variables such that Z i has mean θ/i. It follows that N (n) log n a.s.
as n → ∞ and, moreover, that It is clear that the Beta coalescents belong to a completely different asymptotic regime.
A related problem concerns the infinitely many sites model. Here, as before, we put mutations on the coalescent tree, but this time we imagine that we trace the genealogy of long stretches of chromosome from each of our n individuals. Each time a mutation arrives, it affects a different site on the chromosome. The number of segregating sites is the number of sites at which there exists more than one allele in our sample of chromosomes. This is simply the number of mutations on the skeleton of the coalescent tree. Let S(n) be the number of segregating sites when we start with a sample of n individuals. Clearly the distributions of S(n) and N (n) are related, in that in both cases we count mutations along the skeleton of the coalescent tree; for N (n), we discard any mutation which arises on a lineage all of whose members have already mutated. For the Beta coalescents with α ∈ (1, 2), the asymptotics of S(n) are given in [3] and are the same as those of N (n) given at (1). In [32], Möhle has studied the limiting distribution of S(n) in the case where the measure x −1 Λ(dx) is finite (which includes the Beta coalescents with α ∈ (0, 1)). He proves that where (σ t ) t≥0 is a drift-free subordinator with Lévy measure given by the image under the transformation The number of segregating sites is, in turn, closely related to the length of the coalescent tree (i.e. the sum of the lengths of all of the branches) and to the total number of collisions before absorption. This has been studied for various Λ-coalescents in [3; 4; 11; 13; 22; 26].

The Bolthausen-Sznitman allelic partition
Turning now to the Bolthausen-Sznitman coalescent, Drmota, Iksanov, Möhle and Rösler [13] have proved that log n n S(n) where S(n) is the number of segregating sites. They also give the fluctuations of S n : where a n = n log n + n log log n (log n) 2 , b n = n (log n) 2 and S is a stable random variable with E e itS = exp − 1 2 π|t| + it log t , t ∈ R. The purpose of this paper is to prove the following theorem concerning the complete allele frequency spectrum of the Bolthausen-Sznitman coalescent.
be the number of blocks of the allelic partition of size k when we start with n singleton blocks. Then .
We note that the same asymptotics hold for S(n) and N 1 (n), which bound N (n) above and below respectively. Thus, as a corollary, we also get the asymptotics for N (n): Suppose that we start a general Λ-coalescent (Π(t)) t≥0 from the partition of N into singletons. Then it has been proved by Pitman [34] that either Π(t) has only finitely many blocks for all t > 0 ((Π(t)) t≥0 comes down from infinity) or Π(t) has infinitely many blocks for all time ((Π(t)) t≥0 stays infinite). See Schweinsberg [38] for an explicit criterion for when a Λ-coalcescent comes down from infinity, in terms of the λ b,k 's. The fundamental difference between the Beta coalescents for α ∈ (1, 2) and α ∈ (0, 1] (including the Bolthausen-Sznitman coalescent) is that the former coalescents come down from infinity and the latter do not. This accounts for the fact that in Berestycki, Berestycki and Schweinsberg's result, the scalings are the same for all different sizes of block as n becomes large, whereas in our theorem, the singletons must be scaled differently. Essentially, coalescence occurs rather slowly and the overwhelming first-order effect is mutation, which causes the allelic partition to consist mostly of singletons. However, at the second order (i.e. considering (N 2 (n), N 3 (n), . . .)), we can feel the effect of the coalescence.
We do not claim that our results are of any application in population genetics: to the best of our knowledge, the Bolthausen-Sznitman coalescent has not been used to model the genealogy of any biological population. Nonetheless, our method may extend to the case of coalescents which are more biologically realistic.
Our method of proof is of some interest in itself. We track the formation of the allelic partition using a certain Markov process, for which we then prove a fluid limit (functional law of large numbers). The terminal value of our process gives the allele frequency spectrum and the fluid limit result, after a little extra work, allows us to read off the asymptotics.
Fluid limits have been widely used in the analysis of stochastic networks (see, for example, [8], [39]) and in the study of random graphs ( [9], [36], [40]). In some sense, the prototypical result of the type in which we are interested is the following: suppose we take a Poisson process, (X(t)) t≥0 of rate 1, started from 0. Then the re-scaled process (N −1 X(N t)) t≥0 stays close (in a rather strong sense) to the deterministic function x(t) = t, at least on compact time-intervals. For a general pure jump Markov process, the fluid limit is determined as the solution to a differential equation. In this article we have relied on the neat formulation in Darling and Norris [10]. However, our fluid limit is somewhat unusual. Firstly, instead of scaling time up, we actually scale it down, by a factor of log n. Moreover, we have three different "space" scalings for different co-ordinates of our (multidimensional) process.

Fluid limit
Consider the formation of the allelic partition, starting from the partition into singletons and run until every individual has received a mutation. The easiest way to think of this is to use the terminology of Dong, Gnedin and Pitman [12] in which blocks have two possible states: active and frozen. We start with all blocks active and equal to singletons. Active blocks coalesce according to the rules of the Bolthausen-Sznitman coalescent: if there are b active blocks present then any particular k of them coalesce at rate Moreover, every active block becomes frozen at rate ρ and stays frozen forever (this act of freezing creates a block in the allelic partition).
The data we will track are as follows. Let X n k (t) be the number of active blocks of the coalescent partition at time t containing k individuals, k ≥ 1, where we start at time 0 with n active individuals in singleton blocks. For k ≥ 1, let Z n k (t) be the number of blocks of the allelic partition of size k which have already been formed by time t (this is the number of times so far that an active block containing precisely k individuals has become frozen). For d ≥ 1, let Y n d+1 (t) = ∞ k=d+1 X n k (t), the number of active blocks containing at least d + 1 individuals. It is straightforward to see that, for any d ≥ 1, is a (time-homogeneous) Markov jump process taking values in {0, 1, 2, . . . , n} d+2 , with Fix d ≥ 1 and writeX n,d (t) = (X n 1 (t),X n 2 (t), . . . ,X n d (t),Ȳ n d+1 (t),Z n d (t)) and define a stopping time (Note that T n is the same regardless of the value of d.) For t ≥ 0, let Finally, let We write · for the Euclidean norm on R d+2 .
This is the key to the following result.
Theorem 1.1 now follows directly, since N k (n) = Z n k (T n ) for k ≥ 1. Note that Proposition 2.1 tells us how the allele frequency spectrum is formed.
Remark. Delmas, Dhersin and Siri-Jegousse [11] have recently considered the lengths of coalescent trees associated with Beta coalescents for α ∈ (1, 2). Part (1) of their Theorem 5.1 appears to be a result analogous to our Proposition 2.1.

Asymptotic frequencies
It would be very interesting to have a better understanding of the distribution of the asymptotic frequency sequence of the allelic partition associated with the Bolthausen-Sznitman coalescent. In [18], Gnedin, Hansen and Pitman obtain relations between the total number of blocks N (n) of an exchangeable random partition restricted to the set {1, . . . , n} and the asymptotic form of the sequence (f ↓ i ) i≥1 . More precisely, they prove that, for any α ∈ (0, 1) and any function ℓ : R + → R + , slowly varying at infinity, we have where ℓ * is also a slowly varying function which can be expressed in term of α and ℓ.
It would be nice to have a similar result for the allelic partition associated with the Bolthausen-Sznitman coalescent. There are, however, two main difficulties: first, we would need almost sure convergence of the rescaled process N (·), whereas here we have only established convergence in probability. Second, the Bolthausen-Sznitman coalescent corresponds to the critical case α = 1 for which the first of the above equivalences no longer holds. In this setting, according to Proposition 18 of Gnedin, Hansen and Pitman [18], we have only the implication: and, in addition, that log n n N 1 (n) a.s.
The form of the limits is, of course, basically the same as in our Theorem 1.1 and so we might expect the following result.
be the asymptotic frequency sequence of the allelic partition associated with the Bolthausen-Sznitman coalescent. Then

Fluctuations
Another interesting and natural question concerns the form of the fluctuations around the deterministic limits in Theorem 1.1. The methods used in this paper do not give us access to this information. In view of the fact that S(n) and N (n) have the same first-order asymptotics, one might be tempted to conjecture that they should have the same second-order asymptotics as well. However, we do not have a cogent argument for why this should be the case.

Beta coalescents
The fluid limit methods used in this paper can, in principle, be extended to deal with other classes of coalescent process. For instance, the method seems to work for the Beta coalescents with parameter α ∈ (1, 2). However, the calculations are more complicated than in the Bolthausen-Sznitman case. Indeed, for the Bolthausen-Sznitman coalescent, the active partition is mostly composed of singletons at any time, which essentially enables us to neglect collisions between non-singleton blocks. This approximation does not hold for the Beta coalescents with α ∈ (1, 2). Since the relevant result has already been proved by Berestycki, Berestycki and Schweinsberg [3; 4] by other methods, we will not give the details.
We may also consider the Beta coalescents with parameter α ∈ (0, 1). Möhle's result (2) that the total number of mutations along the coalescent tree, re-scaled by n, converges in distribution to some non-degenerate random variable suggests that here we may expect to have convergence in distribution of the allelic partition to a random vector. Clearly, the fluid limit methods used in the present paper do not adapt to this situation, but we can still use them to investigate the expected value of the number of blocks of different sizes. Indeed, the drift of the re-scaled process converges to an explicit function b (d) (but the varianceᾱ n,d does not tend to 0). This enables us to conjecture that where C 1 , C 2 , . . . are strictly positive random variables. We intend to address this problem in a future paper.

Proofs
In this section, we prove Proposition 2.1 and deduce Proposition 2.2. In order to do so, we use the fluid limit methodology described in Darling and Norris [10]. Firstly, we need to set up some notation. Let β n,d (m) be the drift of the process X n,d when it is in state m = (m 1 , m 2 , . . . , m d+2 ) ∈ {0, ..., n} d+2 , so that where q n,d (m, m ′ ) is the jump rate from m to m ′ . Let α n,d (m) be the corresponding variance of a jump, in the sense that Let us also introduce the notation for 1 ≤ k ≤ d + 2, so that we may decompose α n,d (m) as Finally, let M def = d+1 k=1 m k denote the total number of active blocks in the partition. We will need to compute the drift and infinitesimal variance of the re-scaled processX n,d , which takes values in the set S n,d def = 0, 1 n , . . . , 1 × 0, log n n , 2 log n n , . . . , log n d × 0, (log n) r n , 2 (log n) r n , . . . , (log n) r , where r = 1 if d = 1 and r = 2 if d ≥ 2. Denote byβ n,d (ξ) andᾱ n,d (ξ) the drift and infinitesimal variance ofX n,d when it is in the state ξ = (ξ 1 , ξ 2 , . . . , ξ d+2 ) ∈ S n,d . Then, letting m = (nξ 1 , n log n ξ 2 , . . . , n log n ξ d+1 , n (log n) r ξ d+2 ), we havē Then the vector field b (d) is Lipschitz in the Euclidean norm with constant K def = ρ 2 + π 2 /3. The function x (d) (t) of the previous section is the unique solution of the differential equation In order to prove Proposition 2.1, we need a few lemmas. Firstly, we prove two analytic results. For n ∈ N, let h(n) = n−1 i=1 1 i , the (n − 1)th harmonic number.
Proof. It is an elementary fact that, for k ≥ 2, This entails that in the specified range of x.
Lemma 4.2. For 0 ≤ j ≤ n and k ≥ 0, By the mean value theorem and monotonicity of the logarithm, Hence, the result follows.
We now have the necessary tools to begin proving the fluid limit result.
Fix R ≥ e and let l(n, R, d) = R −1 + d/n and S n,d = ξ ∈ S n,d : It follows that for t 0 < ∞, Proof. We must perform some elementary (but rather involved) calculations. From the rates of the process we will calculate the co-ordinates of β n,d (m) in turn. Recall first that the Λcoalescents are consistent in the sense that, regardless of how many blocks are present in total, if we look just at a subcollection of size b of them, any k of those b blocks coalesce at rate λ b,k . Furthermore, freezing occurs in a Markovian way. Hence, if M active blocks are present in the partition, the next event involves the coalescence of precisely j of them at rate M j λ M,j = M j(j−1) (see Theorem 3.1 of [12]). Thus, we have Then, using the fact that b 1 Thus, For the (d + 1)th co-ordinate we have Note that we have (split the sum on the left-hand side according as b d+1 = 0 or b d+1 = 0). Thus, we get Finally, β n,d d+2 (m) = ρm d .
Using (3) and the notation m = (m 1 , . . . , m d+2 ), we obtain the following expressions: Bearing in mind that M = n ξ 1 + 1 log n d+1 i=2 ξ i and that ξ ∈S n,d (defined at (5)), we can apply Lemma 4.1, to get from (6) that Consider now the sum in the expression (7) forβ n,d k (ξ) when 2 ≤ k ≤ d. We split it into two parts, j = k and 2 ≤ j ≤ k − 1. The j = k term is

By Lemma 4.2 we have
Turning now to the other term, if 2 ≤ j ≤ k − 1, we have With another application of Lemma 4.1, it follows that, for ξ ∈S n,d We turn finally to the expression (8) forβ n d+1 (ξ). Consider the sum which constitutes the third term. We have But then and so, arguing as before, we obtain, for ξ ∈S n,d d+2 (ξ)| = 0. Putting everything together, we obtain that for some constant C(R), whenever ξ ∈S n,d . The final deduction follows easily. It follows that for t 0 < ∞, Proof. Recall that for 1 ≤ k ≤ d + 2 we have We will deal with the co-ordinates in turn.
Proof of Proposition 2.1 We will, in fact, prove the stronger result that, for any d ≥ 1 and any 0 < δ < 1, as n → ∞. Fix d ≥ 1. We follow the method used in Theorem 3.1 of Darling and Norris [10] and start by noting thatX n,d (t) has the following standard decomposition X n,d (t) =X n,d (0) + M n,d (t) + t 0β n,d (X n,d (s))ds, where (M n,d (t)) t≥0 is a zero-mean martingale in the natural filtration ofX n,d . Since Recall that K is the Lipschitz constant of b (d) . Fix R ≥ e and let Ω n,d,1 = and Ω n,d,3 = From (12) we obtain that for t < T R,d n ∧ t 0 and on the event Ω n,d,1 ∩ Ω n,d,2 , Hence, by Gronwall's lemma, n,d (X n,d (s))ds .
In fact, we wish to prove this result for t 0 rather than T R,d n ∧ t 0 . Set Ω n,R,d = sup Since T R,d n = T R,d,1 n ∧ T R,d,2 n , it will suffice for us to show that T R,d,1 n > t 0 and T R,d,2 n > t 0 on Ω n,R,d for all large enough n and R.
Note now that 0 ≤ y 2 (t) ≤ e −1 < R for all t ≥ 0. Take n to be sufficiently big that (log n) The desired result (10) follows.

Proof of Proposition 2.2
For convenience we will write .
Since we can take t 0 arbitrarily large, we can make z d (t 0 ) arbitrarily close to z d (∞). With high probability, on the time interval [0, t 0 ], log n n Z n 1 ( t log n ) stays close to z 1 (t) and, likewise, for d ≥ 2, (log n) 2 n Z n d ( t log n ) stays close to z d (t). So the work of this proof will be to demonstrate that Z n d (t) does not do anything "nasty" between times t 0 /(log n) and T n . (Note that this interval is potentially quite long: absorption for the coalescent takes place at about time log log n; see Proposition 3.4 of [25].) We will have to split the interval [t 0 /(log n), T n ] into two parts and deal with the process separately on each.
(We can do this uniformly in k because of the special form of the functions x 1 (t), y 2 (t) and z k (t), k ≥ 1.) Take 0 < ǫ < δ 3 ∧ δη 72ρ ∧ x 1 (t 0 ). Let Since (log n)T R,d,1 n ≤ T n and ǫ < x 1 (t 0 ), we know by the argument in the proof of Proposition 2.1 that T n > t 0 /(log n) on Ω n,1 for all sufficiently large n and R. Let τ n = inf t ≥ t 0 log n : X n 1 (t) < n (log n) 3 and Y n 2 (t) < n (log n) 3 .
Fix t 1 > 0. For any particular n, A n (t) and Now consider the case d ≥ 2. On Ω n,d we have and, by taking n sufficiently large, we have by Proposition 2.1 that P Ω c n,1 + P Ω c n,d < η/2. Hence, But η was arbitrary and so this completes the proof of Proposition 2.2.