FUNCTIONAL LIMIT THEOREMS FOR MULTITYPE BRANCHING PROCESSES AND GENERALIZED PÓLYA

URNS

SVANTE JANSON

Abstract. A functional limit theorem is proved for multitype continuous time Markov branching processes. As consequences, we obtain limit theorems for the branching process stopped by some stopping rule, for example when the total number of particles reaches a given level.

Using the Athreya–Karlin embedding, these results yield asymptotic results for generalized Pólya urns. We investigate such results in detail and obtain explicit formulas for the asymptotic variances and covariances. The general formulas involve integrals of matrix functions; we show how they can be evaluated and simplified in important special cases. We also consider the numbers of drawn balls of different types and functional limit theorems for the urns.

We illustrate our results by some examples, including several applications to random trees where our theorems and variance formulas give simple proofs of some known results; we also give some new results.

1. Introduction

Consider a generalized Pólya urn process (X_{n})_{n=0}^{∞} defined as follows. (This process is also known as a generalized Pólya–Eggenberger urn or a generalized Friedman urn, cf. [17], [47], [19].) There are balls of q types (or colours) 1, …, q, and each X_{n} is a vector (X_{n1}, …, X_{nq}), where X_{ni} ≥ 0 is the number of balls of type i in the urn at time n. The urn starts with a given vector X_{0}, random or not. We are further given, for each type i, an activity (or weight) a_{i} ≥ 0 (typically a_{i} = 1, but sometimes different a_{i} are useful [2]; we will even find use for a_{i} = 0), and a random q-dimensional vector ξ_{i} = (ξ_{i1}, …, ξ_{iq}) with integer coordinates. (Actually, only the distribution of ξ_{i} matters.) We usually further assume that, almost surely,

ξ_{ij} ≥ 0,  j ≠ i,  (1.1)
ξ_{ii} ≥ −1.  (1.2)

(Relaxation of these requirements, allowing further negative values, will be discussed in Remark 4.2.) The urn evolves according to a Markov process. At each time n ≥ 1, one of the balls in the urn is drawn at random such that the probability of drawing a particular ball of type i is proportional to the activity a_{i}, i.e. the probability of drawing a ball of type i is a_{i}X_{n−1,i}/ ∑_{j} a_{j}X_{n−1,j}. (In particular, if every a_{i} = 1, a ball is drawn uniformly at random.) The drawn ball is returned to the urn together with ΔX_{nj} balls of type j, for each j = 1, …, q, where ΔX_{n} = (ΔX_{n1}, …, ΔX_{nq}) is a random vector such that if the drawn ball has type i, then ΔX_{n} has the same distribution as ξ_{i} and is independent of everything else that has happened so far. (In many applications, the replacement vectors ξ_{i} are deterministic and the randomness enters solely through the the random draws.)

Date: February 20, 2003.

Note that (1.1) means that we may add but never remove balls of other types than the drawn one, while (1.2) means that ΔX_{ni} = −1 is allowed when i is drawn, meaning that the drawn ball is removed (with or without addition of balls of other types). Indeed, the rule above may also be stated as: The drawn ball is removed and, if it had type i, we add a number of new balls with the distribution of (ξ_{ij} + δ_{ij})_{j=1}^{q}. This is often a more natural formulation, and explains (1.2) better. Note that, in both versions, ξ_{i} describes the change in the composition of the urn when a ball of type i is drawn.

If the urn ever becomes empty, or, more generally, there are no balls with non-zero activity left, the process stops (extinction). We are only interested in urns where there is a positive probability of non-extinction, and our main goal is to describe the asymptotics of the urn conditioned on non-extinction. Indeed, in typical applications extinction cannot occur at all.

For some specific examples and applications, see Section 7.

Urn models of this type have been studied by many authors, including [13], [14], [49], [19], [18], [8], [9], [30], [10], [20], [51], [11], [35], [12].

We will use the method of Athreya and Karlin [8], see also [9, §V.9], and study the urn process by embedding it into a multitype continuous time Markov branching process X (t) = (X_{1}(t), …, X_{q}(t)). This process is defined using the same data a_{i} and ξ_{i}, i = 1, …, q, as above and an initial vector X (0) = X_{0}. In this process we assume that a ball (particle) of type i lives an exponentially distributed time with mean a_{i}^{−1}, i.e. it dies with intensity a_{i}, and when it dies, it is replaced by a set of balls with the distribution given by (ξ_{ij} +δ_{ij})_{j=1}^{q}, all life times and offspring compositions being independent. Alternatively, when ξ_{ii} ≥ 0 a.s., the ball lives for ever and at random times according to a Poisson process with intensity a_{i}, it gives birth to a new litter of balls with the distribution given by ξ_{i}. (Unfortunately, the embedding is exemplified in [8] and [9] only with urn processes where the drawn ball is returned, which in our notation means ξ_{ii} ≥ 0. It is clear from [8], [9] that the results hold also when the drawn ball is removed, as observed in e.g. [10], [2], but some authors have overlooked this.) We assume throughout that all variables ξ_{i} have finite mean (and variance, see (A2) below); this is sufficient to prevent explosion and guarantees that the process X (t) exists for all t ≥ 0 [9, §V.7]. We define X (t) to be right-continuous.

Let τ_{0} = 0, and let τ_{n}, n ≥ 1, be the n:th time a ball dies (splits). It is easily shown [8], [9] that the process (X (τ_{n}))_{n=0}^{∞} equals (in distribution)

(X_{n})_{n=0}^{∞}; hence limit theorems for X_{n} can be derived from limit theorems for X (t). The processes X (t), t ≥ 0, and X_{n}, n ≥ 0, are the same up to a random time change (extending the parameter n to real values). However, since X (t) grows exponentially (Lemma 9.8), the time scales are different.

In several applications, see e.g. Examples 7.5 and 7.8, urn processes not satisfying (1.1) and (1.2) appear. Although the embedding method by Athreya and Karlin does not apply directly, it can be modified to handle this case too, in at least two ways. First, in many cases, it is possible to transform the urn process into a different urn process satisfying the conditions above, using a “superball” argument, see Remark 4.2. Secondly, even if ξ_{ij} are arbitrary integers, but we for simplicity assume that they are such that the process never gets stuck, it is possible to define the corresponding continuous time process; this is a generalized branching process where the death of one ball may force the removal of others (like the ancient practices of sacrificing slaves or burning widows), and under appropriate conditions, the results extend to this case too. For simplicity, we nevertheless assume (1.1) and (1.2) in the main parts of the paper, and discuss these extensions in Remark 4.2 and some examples.

Our main results are stated in Section 3; we first introduce some notation and basic assumptions in Section 2. Some extensions, and problems for future research, are discussed in Section 4. We aim at directly applicable results where, in the case of normal limits, the asymptotic variances and covariances are given explicitly, by formulas computable using linear algebra. The formulas in Section 3 are given by integrals. Some simplifications and evaluations of the formulas in important special cases are given in Section 5; see also the examples in Section 7. We discuss how our results and methods relate to some previous papers in Section 6.

In Section 7, we give some examples and applications of our results. In particular, urn processes have been used by several authors to study various classes of random trees, see e.g. the “fringe analysis” in [2]. We review several such applications and show how some old (and a few new) results follow easily from our theorems. See also [29], which may serve as an easier introduction to such applications. Finally, the proofs of the main results are given in Sections 8–10.

Acknowledgements. This research was begun more than ten years ago, inspired by a talk by Hosam Mahmoud at a conference on Random Graphs in Poznań 1991; I thank Hosam Mahmoud and others for interesting discussions over the years. The research has partly been carried out during visits to the Mittag-Leffler Institute in Djursholm, Sweden and the Isaac Newton Institute for Mathematical Sciences in Cambridge, U.K.

2. Preliminaries

We let A denote the q × q matrix

A := (a_{j} 𝔼 ξ_{ji})_{i,j=1}^{q}.  (2.1)

The matrix A and its eigenvalues will play a central role.

Note our choice of notation; in the main case when a_{j} = 1, A_{ij} = 𝔼 ξ_{ji} and the j:th column of A is the expected change when a ball of type j is drawn (splits). It may seem more natural to consider the transpose A′, and this is done by other authors (which should be remembered when comparing the results). The reason for our choice is that we will use the standard notation where a matrix is regarded as an operator acting on column vectors to the right. (In contrast to the standard notation for Markov chains, where the transition matrix acts on row vectors to the left.)

By (1.1), A + αI is a non-negative matrix if α is large enough, so by the standard Perron–Frobenius theory, A has a largest real eigenvalue λ_{1} such that every other eigenvalue λ satisfies Re λ < λ_{1} (see e.g. [50, Chapter 1 and Theorem 2.6] or [32, Appendix 2]). We order the eigenvalues with decreasing real parts: λ_{1} > Re λ_{2} ≥ Re λ_{3} ≥ . . ..

We write i ≻ j if it is possible to find a ball of type j in an urn beginning with a single ball of type i, i.e. if (A^{n})_{ji} > 0 for some n ≥ 0. The relation ≻ is transitive and reflexive, so it partitions the set of all types into equivalence classes C_{1}, …, C_{ν} such that i and j belong to the same class if and only if i ≻ j and j ≻ i; moreover, ≻ induces a partial order among the equivalence classes. We say that a type i is dominating if i ≻ j for every type j; similarly a class C_{k} is dominating if some (and then every) i ∈ C_{k} is dominating.

Note that if we order the classes suitably, and take the types in this order, A becomes a block triangular matrix, see [34] for a detailed treatment. Hence the set of eigenvalues of A (with multiplicities) is the union of the sets of eigenvalues of the restrictions of A to the classes C_{k}; we say that an eigenvalue belongs to a class if it is an eigenvalue of the restriction of A to this class.

The urn or branching process (or A) is irreducible (or positive regular, which is equivalent in continuous time) if there is only one equivalence class, i.e. if i ≻ j for any types i and j; equivalently, every type is dominating.

We are mainly interested in the irreducible case, but for an important technical reason, see for example the proof of Theorem 3.16, we will state our results somewhat more generally.

Our basic assumptions are as follows (but see Remark 4.2):

(A1) (1.1) and (1.2) hold, i.e. ξ_{ij} + δ_{ij} ≥ 0 a.s. for all i and j.

(A2) 𝔼 ξ_{ij}^{2} < ∞ for all i, j = 1, …, q.

(A3) The largest real eigenvalue λ_{1} of A is positive, λ_{1} > 0.

(A4) The largest real eigenvalue λ_{1} is simple.

(A5) There exists a dominating type i with X_{0i} > 0 (X (0)_{i} > 0), i.e. we

start with at least one ball of a dominating type. (A6) λ_{1} belongs to the dominating class.

We assume that the classes are ordered so that C_{1} is the dominating class.

(A1) is already discussed. (A2) is essential since we use L^{2} theory. (A3) says that the branching process X (t) is supercritical, and implies that X (t), and thus X_{n}, has a positive probability of non-extinction. (Non-extinction is also possible in some exceptional cases with λ_{1} = 0, for example when the balls always have exactly one child, and thus change type according to a Markov chain. We do not treat these cases.)

Note that (A4), (A5) and (A6) hold when A is irreducible [32], [50]. In the reducible case, (A5) is only a weak restriction; if we consider a case with a single ball initially, we may ignore all types that cannot occur and then (A5) holds. (A4) and (A6) are more significant restrictions; see [34] (which treats the related case of multitype Galton–Watson processes in discrete time) or [35] for some complications that otherwise can occur.

We say that the process becomes essentially extinct if at some time there are no balls of any dominating type left. Note that if we restrict attention to the balls of the dominating types, we have an irreducible multitype Galton– Watson process, and essential extinction means that this restricted process becomes extinct. For irreducible processes (our main concern), essential extinction is thus the same as extinction.

For most applications, the following lemma yields convenient criteria (possibly combined with Remark 4.2).

Lemma 2.1. If A is irreducible, (A1) and (A2) hold, ∑_{j} 𝔼 ξ_{ij} ≥ 0 for every i and ∑_{j} 𝔼 ξ_{ij} > 0 for some i, then (A1)–(A6) hold and (essential) extinction is impossible.

Proof. The conditions imply that the total number of balls never decreases, which guarantees non-extinction. Since the process is irreducible, this is the same as essential non-extinction. We have already remarked that (A4)–(A6) hold when A is irreducible. Finally, it is easy to see that the conditions imply (A3), cf. [50, Theorem 1.1, Corollary 1]. □

We collect various facts and notations that will be used throughout the paper, usually without further comment.

We will often regard A and other matrices as linear operators in ℝ^{q} or ℂ^{q}, or in some invariant subspace thereof. In this context, vectors in ℂ^{q}, in particular X_{n} and X (t), are always regarded as column vectors. Consequently, by an eigenvector of A we mean a right eigenvector; a left eigenvector is the same as an eigenvector of the transpose matrix A′.

Note that if u and v are vectors in ℂ^{q}, then u′v is a scalar while uv′ is a q × q matrix. We also use the notation u · v for u′v.

We use | · | for the norm of both vectors and matrices. (The choice of matrix norm is irrelevant.)

We let a denote the (column) vector (a_{1}, …, a_{q}) of activities, and let u_{1} and v_{1} denote left and right eigenvectors of A corresponding to the largest eigenvalue λ_{1}, i.e. vectors satisfying

u′_{1}A = λ_{1}u′_{1}, Av_{1} = λ_{1}v_{1}.

By (A4) u_{1} and v_{1} are unique up to scalar factors, and by the Perron– Frobenius theory [32], [50], (applied to A + αI for suitable α), they may be chosen non-negative.

If the process is irreducible, all entries of u_{1} and v_{1} are strictly positive [32], [50]. In general, it follows easily from this result applied to the restriction to the dominating class C_{1} together with (A6) that v_{1i} > 0 for every i while u_{1i} > 0 if i ∈ C_{1} (i.e., i is dominating) and u_{1i} = 0 otherwise.

The scalar products u_{1} · v_{1} and a · v_{1} thus are both positive, and we may assume that v_{1} and u_{1} are normalized such that

a · v_{1} = a′v_{1} = v′_{1}a = 1,  (2.2)
u_{1} · v_{1} = u′_{1}v_{1} = v′_{1}u_{1} = 1.  (2.3)

This determines u_{1} and v_{1}, and we fix this choice of u_{1} and v_{1} throughout the paper. (Otherwise, obvious normalization factors would enter into some formulas.)

We will use the Jordan decomposition of the matrix A in the following form, see e.g. [46, Theorem 7.6]: There exists a decomposition of the complex space ℂ^{q} as a direct sum ⨁ E_{λ} of generalized eigenspaces E_{λ}, such that A−λ is a nilpotent operator on E_{λ}; here λ ranges over the set Λ of eigenvalues of A. (A − λ denotes A − λI, where I is the identity matrix of appropriate size.) In other words, there exist projections P_{λ}, λ ∈ Λ, that commute with A and satisfy

∑_{λ∈Λ} P_{λ} = I,  (2.4)
AP_{λ} = P_{λ}A = λP_{λ} + N_{λ},  (2.5)

where N_{λ} = P_{λ}N_{λ} = N_{λ}P_{λ} is nilpotent. Moreover, P_{λ}P_{μ} = 0 when λ ≠ μ. We let d_{λ} ≥ 0 be the integer such that N_{λ}^{d_{λ}} ≠ 0 but N_{λ}^{d_{λ}+1} = 0. (Equivalently, in the Jordan normal form of A, the largest Jordan block with λ on the diagonal has size d_{λ} + 1.) Hence d_{λ} = 0 if and only if N_{λ} = 0, and this happens for all λ if and only if A is diagonalizable, i.e. if and only if A has a complete set of q linearly independent eigenvectors.

Note, by taking transposes in (2.4) and (2.5), that P′_{λ} and N′_{λ} are the corresponding projections and nilpotent operators for A′.

We define, for k = 0, 1, …, the quotient space E_{λ,k} := E_{λ}/N_{λ}^{k+1}E_{λ} and the projection Q_{λ,k} : E_{λ} → E_{λ,k}, noting that E_{λ,d_{λ}} = E_{λ} and Q_{λ,d_{λ}} = I. Then N_{λ} : E_{λ} → E_{λ} induces a map N_{λ} : E_{λ,k} → E_{λ,k+1}, and if 0 ≤ j ≤ k,

N_{λ}^{j}Q_{λ,k−j} = Q_{λ,k}N_{λ}^{j} : E_{λ} → E_{λ,k}.  (2.6)

Since we assume that λ_{1} is a simple eigenvalue, N_{λ_{1}} = 0 and d_{λ_{1}} = 0, and P_{λ_{1}} is the one-dimensional projection

P_{λ_{1}} = v_{1}u′_{1}.  (2.7)

In the sequel, λ will always denote an eigenvalue of A. (Formally, the results hold for other λ too if we then set P_{λ} = N_{λ} = 0.)

We recall that matrix exponentials may be defined by power series; for example e^{tA} = ∑_{j=0}^{∞} t^{j}A^{j}/j!. We have, using (2.5) and commutativity,

\begin{aligned}P_{λ}e^{tA} & = e^{tA}P_{λ} = P_{λ} ∑_{j=0}^{∞} \frac{t^{j}}{j!}(P_{λ}A)^{j} = P_{λ}e^{tP_{λ}A} = P_{λ}e^{λtP_{λ}+tN_{λ}} \\ & = P_{λ}e^{λtP_{λ}}e^{tN_{λ}} = P_{λ}e^{λt} ∑_{j=0}^{d_{λ}} \frac{t^{j}}{j!}N_{λ}^{j} \end{aligned}  (2.8)

and thus, by (2.4),

e^{tA} = ∑_{λ} ∑_{j=0}^{d_{λ}} \frac{t^{j}}{j!}e^{λt}P_{λ}N_{λ}^{j}.  (2.9)

Some immediate consequences are

|P_{λ}e^{tA}| ≤ C(1 + |t|)^{d_{λ}}e^{Re λt},  −∞ < t < ∞,  (2.10)

and more generally, for 0 ≤ k ≤ d_{λ},

|Q_{λ,k}P_{λ}e^{tA}| ≤ C(1 + |t|)^{k}e^{Re λt},  −∞ < t < ∞,  (2.11)

and, using (A4),

|e^{tA}| ≤ Ce^{λ_{1}t},  0 ≤ t < ∞,  (2.12)

where, as sometimes later, C denotes unspecified constants that may depend on the data q, a_{i}, ξ_{i}, X_{0}.

As is well-known since decades [33], [5], [6], the asymptotic behavior depends on whether there is any eigenvalue beside λ_{1} with a real part > λ_{1}/2. We thus define Λ_{I} := {λ ∈ Λ : Re λ < λ_{1}/2}, Λ_{II} := {λ ∈ Λ : Re λ = λ_{1}/2}, Λ_{III} := {λ ∈ Λ : Re λ > λ_{1}/2}; hence Λ is the disjoint union Λ_{I} ∪Λ_{II} ∪Λ_{III}. We further define P_{I} := ∑_{λ∈Λ_{I}} P_{λ}, the projection onto the sum of the generalized eigenspaces with Re λ < λ_{1}/2.

For later use, we define the following matrices.

B_{i} := 𝔼(ξ_{i}ξ′_{i}),  (2.13)
B := ∑_{i=1}^{q} v_{1i}a_{i}B_{i},  (2.14)

Σ_{I} := ∫_{0}^{∞} P_{I}e^{sA}Be^{sA′}P′_{I}e^{−λ_{1}s} ds,  (2.15)
Σ_{II} := ∑_{λ∈Λ_{II}} P_{λ}BP ′ ; λ̄ = ∑_{λ∈Λ_{II}} P_{λ}BP_{λ}^{∗},  (2.16)

Σ_{II,d} := \frac{1}{(2d + 1)d!^{2}} ∑_{λ∈Λ_{II}} N_{λ}^{d}P_{λ}BP_{λ}^{∗}(N_{λ}^{∗})^{d},  (2.17)

where ^{∗} denotes Hermite conjugation and d = 0, 1, …; thus Σ_{II,0} = Σ_{II}.

Remark 2.2. Let ξ_{∗} be the random vector obtained by choosing ξ_{i} for a random type i, with the probability a_{i}v_{1i} for type i. It follows from Theorem 3.21 below that this is the asymptotic distribution of the drawn types; thus ξ_{∗} is the asymptotic distribution of the added balls.

Then (2.14) says B = 𝔼 ξ_{∗}ξ′_{∗}. Since, 𝔼 ξ_{∗} = ∑_{i=1}^{q} v_{1i}a_{i} 𝔼 ξ_{i} and, by (2.1),

∑_{i=1}^{q} v_{1i}a_{i} 𝔼 ξ_{i} = (∑_{i=1}^{q} v_{1i}a_{i} 𝔼 ξ_{ij})_{j=1}^{q} = (∑_{i=1}^{q} v_{1i}A_{ji})_{j=1}^{q} = Av_{1} = λ_{1}v_{1},

(2.18)

the covariance matrix of ξ_{∗} is given by

B̂ := 𝔼 ξ_{∗}ξ′_{∗} − 𝔼 ξ_{∗} 𝔼 ξ′_{∗} = B − λ_{1}^{2}v_{1}v′_{1}.

Not surprisingly, this quantity will appear below. Indeed, since P_{λ}v_{1} = 0 when λ ≠ λ_{1}, we may replace B by B̂ in (2.15)–(2.17); this might be conceptually better, but we prefer B for computations.

We state our results on convergence of stochastic processes using the usual function space D of right-continuous functions with left-hand limits, always equipped with the Skorohod J_{1}-topology. Our processes will, however, be defined on several different intervals, so we will use several versions of D; we will also consider vector-valued processes. In general, for a finite-dimensional vector space E and any (open, closed or half-open) interval J ⊆ [−∞, ∞], we let D(J) = D(J, E) be the space of all right-continuous functions J → E with left-hand limits. We say that f_{n} → f in D if there exists strictly increasing continuous maps λ_{n} of J onto itself such that λ_{n} → ι (the identity map) and f_{n} ∘ λ_{n} → f uniformly on compact subsets of J. When f is continuous, this is equivalent to f_{n} → f uniformly on compact subsets of J. It is well-known that this topology is Polish, i.e. defined by some complete metric. (The case J = [0, 1] is discussed in detail in [15], and the case J = [0, ∞) in [37], [24]. See also [27].) Note that both the space and the topology are changed if we add or remove an endpoint of J. If Z_{n} →^{d} Z in D(J) for some processes Z_{n} and Z defined on J, and Z is a.s. continuous, then the restrictions to any subinterval J′ ⊂ J converge in D(J′).

3. Results

The basis of all our results for the branching process and generalized Pólya urns, is the following functional limit theorem for the branching process. (For previously known results, including parts of this theorem, see Section 6.) Since different normalizations (and different time scalings) are required for different components of X (t), the result is stated in terms of various projections of X (t); this is equivalent to stating results for scalar products η · X (t), where the normalization depends on η, as is done by several other authors. The proof is given in Section 9.

Theorem 3.1. Assume (A1)–(A6). Then, as t → ∞, e^{−λ_{1}t}X (t) →^{a.s.} W v_{1}, and, with joint convergence in distribution of all processes,

(i) e^{−λ_{1}(t+x)/2}P_{I}X (t + x) →^{d} W^{1/2}U_{I}(x) in D(−∞, ∞); equivalently,

e^{−λ_{1}t/2}P_{I}X (t + x) →^{d} W^{1/2}e^{λ_{1}x/2}U_{I}(x) in D(−∞, ∞);

(ii) for every λ ∈ Λ_{II} and k = 0, …, d_{λ},

t^{−(k+1/2)}e^{−λxt}Q_{λ,k}P_{λ}X (xt) →^{d} W^{1/2}U_{λ,k}(x) in D[0, ∞);

(iii) for every λ ∈ Λ_{III} and k = 0, …, d_{λ}, t^{−k}e^{−λt}Q_{λ,k}P_{λ}X (t) →^{a.s.} W_{λ,k}

and thus

(xt)^{−k}e^{−λxt}Q_{λ,k}P_{λ}X (xt) →^{d} W_{λ,k} in D(0, ∞).

Here, U_{I} and U_{λ,k} are continuous Gaussian vector-valued stochastic processes, with U_{I}(x) defined for −∞ < x < ∞ and U_{λ,k}(x) defined for 0 ≤ x < ∞, while W_{λ,k} are vector-valued random variables, also regarded as constant stochastic processes. Moreover, W is a non-negative random variable, related to W_{λ_{1}} := W_{λ_{1},0} by W = u_{1} · W_{λ_{1}} and W_{λ_{1}} = W v_{1}.

The process U_{I} is real, while the processes U_{λ,k} and variables W_{λ,k} are real for real λ but complex otherwise, with U_{λ̄,k} = \overline{U_{λ,k}} and W_{λ̄,k} = \overline{W_{λ,k}}. Furthermore, a.s., U_{I}(x) ∈ E_{I} := ⨁_{λ∈Λ_{I}} E_{λ}, U_{λ,k}(x) ∈ E_{λ,k}, and W_{λ,k} ∈ E_{λ,k}.

The process U_{I}, the families {U_{λ,k}}_{0≤k≤d_{λ}} for different λ ∈ Λ_{II} with Im λ ≥ 0, and the family {W_{λ,k}}_{λ∈Λ_{III},k≤d_{λ}} ∪ {W } are independent of each other.

The processes U_{I} and U_{λ,k} are characterized by being (jointly) Gaussian with mean 0 and covariances, for 0 ≤ x ≤ y,

𝔼(U_{I}(x)U_{I}(y)′) = Σ_{I}e^{(y−x)(A′−λ_{1}/2)},  (3.1)

and

𝔼(U_{λ,k}(x)U_{μ,l}(y)′) = c(k, l, x, y)Q_{λ,k}N_{λ}^{k}P_{λ}Σ_{II}P′_{μ}(N′_{μ})^{l}Q′_{μ,l}  (3.2)

(which vanishes unless μ = λ̄), where

\begin{aligned}c(k, l, x, y) & := \frac{1}{k! l!} ∫_{0}^{x} s^{k}(y − x + s)^{l} ds \\ & = \frac{1}{k! l!} ∑_{j=0}^{l} \binom{l}{j}(y − x)^{j} \frac{x^{k+l+1−j}}{k + l + 1 − j} . \end{aligned}  (3.3)

The results above holds also if we condition X and W, W_{λ,k} on W > 0, or (which is a.s. the same) on essential non-extinction.

Remark 3.2. When x = y, (3.3) simplifies to

c(k, l, x, x) = \frac{x^{k+l+1}}{k! l! (k + l + 1)}.

Remark 3.3. Taking k = d_{λ} in (ii) or (iii), we have Q_{λ,k} = I and we thus find the limit of P_{λ}X under appropriate normalization. The point of the variable k is that, when d_{λ} > 0 (i.e., when the nilpotent part N_{λ} does not vanish), some linear combinations of the components of P_{λ}X have smaller asymptotic variance than others, and thus they disappear in the normalization required for P_{λ}X. More precisely, if η ∈ E′_{λ} is such that (N′_{λ})^{k+1}η = 0, then η can be regarded as a linear functional on E_{λ,k}, and if, for example, λ ∈ Λ_{II}, (ii) shows that t^{−k−1/2}e^{−λ_{1}t/2}η · X (xt) converges to a Gaussian process.

Remark 3.4. If λ ∈ Λ_{II} and d_{λ} = 0, (ii) simplifies to t^{−1/2}e^{−λxt}P_{λ}X (xt) →^{d} W^{1/2}U_{λ}(x), with, from (3.2) and (3.3), for 0 ≤ x ≤ y,

𝔼(U_{λ}(x)U_{λ}^{∗}(y)) = xP_{λ}Σ_{II}P_{λ}^{∗} = xP_{λ}BP_{λ}^{∗}.

In this case, U_{λ} is a process with independent increments.

Remark 3.5. By (3.1), U_{I} is a stationary Gaussian process. It can be regarded as a multi-dimensional Ornstein–Uhlenbeck process.

Remark 3.6. If λ ∈ Λ_{II} with Im λ ≠ 0, it follows from (3.2) and (2.16) that 𝔼(U_{λ,k}(x)U_{λ,k}(x)′) = 0. Hence U_{λ,k}(x) is a vector-valued symmetric com-

plex Gaussian random variable, i.e. ωU_{λ,k}(x) =^{d} U_{λ,k}(x) for every complex number ω with |ω| = 1, see [28, Proposition 1.31].

Consequently, U_{λ,k} is either real (when λ is real) or symmetric complex.

Remark 3.7. We have no general description of the distributions of W_{λ,k} for λ ∈ Λ_{III}, and there seems to be no reason to expect any. They are (typically, at least) not normal, and not independent of each other. Moreover, their distributions (typically) depend on the initial state X (0), unlike U_{I} and U_{λ,k}.

Taking x = 0 in (i) and x = 1 in (ii) and (iii), we obtain as a corollary some standard results, cf. [9].

Corollary 3.8. Assume (A1)–(A6). Then, as t → ∞, with joint convergence,

(i) e^{−λ_{1}t/2}P_{I}X (t) →^{d} W^{1/2}U_{I};

(ii) for every λ ∈ Λ_{II} and k = 0, …, d_{λ},

t^{−(k+1/2)}e^{−λt}Q_{λ,k}P_{λ}X (t) →^{d} W^{1/2}U_{λ,k};

(iii) for every λ ∈ Λ_{III} and k = 0, …, d_{λ}, t^{−k}e^{−λt}Q_{λ,k}P_{λ}X (t) →^{d} W_{λ,k}

(and →^{a.s.}).

Here, U_{I}, U_{λ,k} and W_{λ,k} are vector-valued random variables with U_{I} and U_{λ,k} jointly Gaussian. The vector U_{I} is real, while U_{λ,k} and W_{λ,k} are real for real λ but complex otherwise, with U_{λ̄,k} = \overline{U_{λ,k}} and W_{λ̄,k} = \overline{W_{λ,k}}. Furthermore, a.s., U_{I} ∈ E_{I} := ⨁_{λ∈Λ_{I}} E_{λ}, U_{λ,k} ∈ E_{λ,k}, and W_{λ,k} ∈ E_{λ,k}.

The random vector U_{I}, the families {U_{λ,k}}_{0≤k≤d_{λ}} for different λ ∈ Λ_{II} with Im λ ≥ 0, and the family {W_{λ,k}}_{λ∈Λ_{III},k≤d_{λ}} ∪ {W } are independent of each other.

The Gaussian random vectors U_{I} and U_{λ,k} are characterized by mean 0 and covariances

𝔼(U_{I}U_{I}) = Σ_{I},
𝔼(U_{λ,k}U′_{μ,l}) = \frac{1}{(k + l + 1) k! l!}Q_{λ,k}N_{λ}^{k}P_{λ}Σ_{II}P′_{μ}(N′_{μ})^{l}Q′_{μ,l}

(which vanishes unless μ = λ̄). □

Limits at stopping times. The main interest in the functional limit theorem above is that it enables us to study X (t) at random times t. Our main interest is to let t be the n:th splitting time τ_{n}; as mentioned in the introduction, this yields results for the urn process X_{n}. As a preparation for this, we will first study another important example where we stop the process when we reach a given total number of balls, or a given number of balls of a given type. Somewhat more generally, let b ∈ ℝ^{q} be a fixed vector and define for z ≥ 0

τ_{b}(z) := min{t ≥ 0 : b · X (t) ≥ z} (with min ∅ = +∞).

We assume b · v_{1} > 0, which means that typically b · X (t) → ∞, as is shown by the following lemma. A more precise result is given in Lemma 10.1. The proof of this lemma and the following results are given in Section 10.

Lemma 3.9. Assume b · v_{1} > 0. Conditioned on essential non-extinction, we have a.s. b · X (t) → ∞ as t → ∞ and thus 0 ≤ τ_{b}(z) < ∞ for all z ≥ 0. Moreover, τ_{b}(z) → ∞ as z → ∞.

Theorem 3.10. Assume (A1)–(A6) and let b ∈ ℝ^{q} with b · v_{1} > 0. Conditioned on essential non-extinction we have as z → ∞,

z^{−1}X (τ_{b}(z)) →^{a.s.} (b · v_{1})^{−1}v_{1},  (3.4)

and, with joint convergence in distribution of all random vectors,

(i) z^{−1/2}P_{I}X (τ_{b}(z)) →^{d} (b · v_{1})^{−1/2}V_{I};

(ii) for every λ ∈ Λ_{II} and k = 0, …, d_{λ},

(z ln^{2k+1} z)^{−1/2}Q_{λ,k}P_{λ}X (τ_{b}(z)) →^{d} (b · v_{1})^{−1/2}V_{λ,k};

(iii) for every λ ∈ Λ_{III} and k = 0, …, d_{λ},

(ln z)^{−k}z^{−λ/λ_{1}}Q_{λ,k}P_{λ}X (τ_{b}(z)) →^{a.s.} (b · v_{1})^{−λ/λ_{1}} W̌_{λ,k};

(iv)

z^{−1/2} (P_{λ_{1}}X (τ_{b}(z)) − \frac{z}{b · v_{1}} v_{1} + ∑_{λ≠λ_{1}} \frac{b · P_{λ}X (τ_{b}(z))}{b · v_{1}}v_{1}) →^{a.s.} 0;

(v)

z^{−1/2} (X (τ_{b}(z)) − \frac{z}{b · v_{1}} v_{1} − (I − \frac{v_{1}b′}{b · v_{1}} ) ∑_{λ≠λ_{1}} P_{λ}X (τ_{b}(z))) →^{a.s.} 0.

Here, V_{I}, V_{λ,k} and W̌_{λ,k} are vector-valued random variables with V_{I} and V_{λ,k} jointly Gaussian. The vector V_{I} is real, while the vectors V_{λ,k} and W̌_{λ,k} are real for real λ but complex otherwise, with V_{λ̄,k} = \overline{V_{λ,k}} and W̌_{λ̄,k} = \overline{W̌_{λ,k}}. Furthermore, a.s., V_{I} ∈ E_{I} := ⨁_{λ∈Λ_{I}} E_{λ}, V_{λ,k} ∈ E_{λ,k}, and W̌_{λ,k} ∈ E_{λ,k}.

The random vector V_{I}, the families {V_{λ,k}}_{0≤k≤d_{λ}} for different λ ∈ Λ_{II} with Im λ ≥ 0, and the family {W_{λ,k}}_{λ∈Λ_{III},k≤d_{λ}} are independent of each other.

The Gaussian vectors V_{I} and V_{λ,k} are characterized by mean 0 and covariances

𝔼(V_{I}V′_{I}) = Σ_{I},  (3.5)
𝔼(V_{λ,k}V′_{μ,l}) = \frac{λ_{1}^{−k−l−1}}{(k + l + 1) k! l!}Q_{λ,k}N_{λ}^{k}P_{λ}Σ_{II}P′_{μ}(N′_{μ})^{l}Q′_{μ,l}  (3.6)

(which vanishes unless μ = λ̄).

The same results hold for X_{τ_{b}(z)} in the urn process.

We specialize to some important cases.

Corollary 3.11. Assume (A1)–(A6) and let b ∈ ℝ^{q} with b · v_{1} > 0. Suppose further Re λ_{2} < \frac{1}{2}λ_{1}. Conditioned on essential non-extinction we have as z → ∞,

z^{−1/2}(X (τ_{b}(z)) − \frac{z}{b · v_{1}} v_{1}) →^{d} N(0, Σ_{b}),

where the covariance matrix Σ_{b} is given by

\begin{aligned}Σ_{b} & = (b · v_{1})^{−1}(I − \frac{v_{1}b′}{b · v_{1}} )Σ_{I}(I − \frac{bv′_{1}}{b · v_{1}} ) \\ & = (b · v_{1})^{−1} ∫_{0}^{∞}(I − \frac{v_{1}b′}{b · v_{1}} )e^{sA}Be^{sA′}(I − \frac{bv′_{1}}{b · v_{1}} )e^{−λ_{1}s}ds. \end{aligned}  (3.7)

The same result holds for X_{τ_{b}(z)} in the urn process.

Corollary 3.12. Assume (A1)–(A6) and let b ∈ ℝ^{q} with b · v_{1} > 0. Suppose further Re λ_{2} = \frac{1}{2}λ_{1}, and let d := max\{d_{λ} : Re λ = \frac{1}{2}λ_{1}\}. Conditioned on essential non-extinction we have as z → ∞,

(z ln^{2d+1} z)^{−1/2}(X (τ_{b}(z)) − \frac{z}{b · v_{1}} v_{1}) →^{d} N(0, Σ_{b}),

where the covariance matrix Σ_{b} is given by

Σ_{b} = \frac{λ_{1}^{−2d−1}}{b · v_{1}} (I − \frac{v_{1}b′}{b · v_{1}} )Σ_{II,d}(I − \frac{bv′_{1}}{b · v_{1}} ).  (3.8)

The same result holds for X_{τ_{b}(z)} in the urn process.

Corollary 3.13. Assume (A1)–(A6) and let b ∈ ℝ^{q} with b · v_{1} > 0. Suppose further Re λ_{2} > \frac{1}{2}λ_{1}, and let d := max\{d_{λ} : Re λ = Re λ_{2}\}. Conditioned on essential non-extinction, the family of random variables

Y_{b}(z) := (ln z)^{−d}z^{− Re λ_{2}/λ_{1}}(X (τ_{b}(z)) − \frac{z}{b · v_{1}} v_{1}),

for z ≥ 2, say, is tight. More precisely, there exist complex random vectors W_{b,λ}, λ ∈ Λ′_{III} := {λ : Re λ = Re λ_{2}, Im λ ≥ 0 and d_{λ} = d}, such that, as z → ∞,

Y_{b}(z) − Re ∑_{λ∈Λ′_{III}} e^{−i(Im λ/λ_{1}) ln z}W_{b,λ} →^{a.s.} 0.  (3.9)

In particular, if λ_{2} is real and Λ′_{III} = {λ_{2}}, then as z → ∞, Y_{b}(z)−W_{b} →^{a.s.} 0 for some random vector W_{b}, and thus

Y_{b}(z) →^{d} W_{b}.  (3.10)

The same results hold for X_{τ_{b}(z)} in the urn process.

Remark 3.14. Unless all ξ_{i} a.s. belong to some fixed subspace of ℝ^{q}, the matrix B has full rank q, and thus the covariance matrix in (3.7) has rank q − 1; hence the limit distribution is concentrated on the hyperplane {v : b · v = 0} but not on any smaller subspace.

In contrast, the limits in Corollaries 3.12 and 3.13 are typically concentrated on a subspace of ℝ^{q} of low dimension (commonly 1 or 2).

Remark 3.15. Suppose that Re λ_{2} > \frac{1}{2}λ_{1} with Im λ_{2} > 0, and that, for simplicity, Λ′_{III} = {λ_{2}}. In that case, the sum in (3.9) contains only one term, and we see that Y_{b}(z) converges in distribution when z → ∞ along a subsequence where the fractional part of (Im λ_{2}/2πλ_{1}) ln z converges. The subsequence limits are of the form Re(e^{iα}W_{b}) for a complex random vector

W_{b} and real α. Hence, a limit distribution as in (3.10) exists only if ωW_{b} =^{d} W_{b} whenever |ω| = 1, i.e. if W_{b} is symmetric complex. This seems highly unlikely, and we conjecture that it never happens, so we do not expect that Y_{b} converges in distribution in this case (or, more generally, when Λ′_{III} contains a non-real λ). Unfortunately, we have not been able to show this conjecture in general, but in some particular cases [25], [16], non-convergence of Y_{b}(z) has been shown by computations of the first moments, showing that the moments oscillate.

Similarly, we conjecture that the (subsequence) limits in Corollary 3.13 never are normal; again we cannot prove this, but it can be verified in some particular cases by computations of moments [14], [49], [25], [16].

Limits for urns. In some urn processes we add a fixed number of balls each time (with varying types), say m; assuming as we may that we start with a non-random number of balls l, the total number of balls at time n then is deterministic mn + l and τ_{n} equals τ_{b}(mn + l) with b = (1, 1, …, 1)′ (or τ_{b}(n + l/m) with b = m^{−1}(1, 1, …, 1)′). In this case, Theorem 3.10 and its corollaries thus yield results for X (τ_{n}) and thus for X_{n}.

In general, we can change the setup by adding a dummy type q + 1 such that a_{q+1} = 0 (the dummy balls never split) and a dummy ball is added whenever a ball splits (see Section 10 for details). The dummy balls then count the number of splits, and Theorem 3.10 yields a description of the asymptotics of X_{n}. In order to give explicit expressions for the asymptotic variances, we define, for s ≥ 0,

ϕ(s, A) := ∑_{n=1}^{∞} \frac{s^{n}}{n!} A^{n−1} = ∫_{0}^{s} e^{tA}dt,  (3.11)
ψ(s, A) := e^{sA} − λ_{1}v_{1}a′ϕ(s, A).  (3.12)

The first result below is in [9, Section V.9.3] (under slightly different hypotheses), and is included for completeness.

Theorem 3.16. Assume (A1)–(A6). Conditioned on essential non-extinc-

tion, n^{−1}X_{n} →^{a.s.} λ_{1}v_{1} as n → ∞.

Theorem 3.17. Assume (A1)–(A6). Suppose further Re λ_{2} < \frac{1}{2}λ_{1}. Conditioned on essential non-extinction we have as n → ∞,

n^{−1/2}(X_{n} − nλ_{1}v_{1}) →^{d} N(0, Σ),

where the covariance matrix Σ is given by

Σ = ∫_{0}^{∞} ψ(s, A)Bψ(s, A)′e^{−λ_{1}s}λ_{1} ds − λ_{1}^{2}v_{1}v′_{1}.  (3.13)

Theorem 3.18. Assume (A1)–(A6). Suppose further Re λ_{2} = \frac{1}{2}λ_{1}, and let d := max\{d_{λ} : Re λ = \frac{1}{2}λ_{1}\}. Conditioned on essential non-extinction we have as z → ∞,

(n ln^{2d+1} n)^{−1/2}(X_{n} − nλ_{1}v_{1}) →^{d} N(0, Σ),

where the covariance matrix Σ is given by

Σ = λ_{1}^{−2d}(I − T)Σ_{II,d}(I − T′),  (3.14)

with T := ∑_{λ∈Λ_{II}} λ^{−1}λ_{1}v_{1}a′P_{λ}. If a ∈ Im(A′) and a = A′â, T can be replaced by T_{1} := λ_{1}v_{1}â′.

Theorem 3.19. Assume (A1)–(A6). Suppose further Re λ_{2} > \frac{1}{2}λ_{1}, and let d := max{d_{λ} : Re λ = Re λ_{2}}. Conditioned on essential non-extinction, the family of random variables

Y (n) := (n ln^{2d} n)^{−1/2}(X_{n} − nλ_{1}v_{1}), n ≥ 2,

is tight, and we have the same type of asymptotic behaviour as described in Corollary 3.13 and Remark 3.15. In particular, if λ_{2} is real and Re λ_{3} < λ_{2}, or more generally when Λ′_{III} = \{\frac{1}{2}λ_{1}\}, then Y (n) →^{d} W_{b} for some random vector W_{b}, but we conjecture that this fails otherwise.

Remark 3.20. In Theorem 3.18, the limit is typically concentrated on a subspace L of low dimension (as in Remark 3.14). For some vectors η ∈ ℝ^{q} with η ⊥ L (and thus ηΣη′ = 0) it is possible to obtain a non-degenerate limit in distribution of η · X_{n} with a different normalization from Theorem 3.10(i) or (ii); we leave the details to the reader. The same applies to Theorem 3.19 (and subsequence limits). Cf. [8], [9].

The drawn balls. Let N_{ni} be the number of drawn balls of type i in the first n draws, and let N_{n} := (N_{n1}, …, N_{nq}). (In the branching process, we would similarly study N_{i}(t), the number of deaths of balls of type i up to time t. We leave this case to the reader.)

If each ξ_{i} is deterministic and, for notational simplicity, each a_{i} = 1, then X_{n} = X_{0} + AN_{n}, so if further A is invertible, asymptotics for N_{n} = A^{−1}(X_{n} − X_{0}) follow from Theorems 3.16–3.19. In general, we can argue with dummy balls again, this time using q dummy types q + 1, …, 2q and adding a dummy ball of type q + i each time a ball of type i is drawn. This leads to the following theorem. (Explicit expressions when some a_{i} ≠ 1 can be found by this method too, but are left to the reader.)

Theorem 3.21. Assume (A1)–(A6). Conditioned on essential non-extinction we have, as n → ∞, n^{−1}N_{ni} →^{a.s.} ν_{i} := a_{i}v_{1i} and, furthermore,

(i) if Re λ_{2} < λ_{1}/2, then

n^{−1/2}(X_{n} − nλ_{1}v_{1}, N_{n} − nν) →^{d} (V, V̂),

where (V, V̂) is vector-valued Gaussian random variable with mean 0; if for simplicity each a_{i} = 1, and D_{v} is the diagonal matrix with entries D_{ii} = v_{1i},

𝔼 V̂V̂′ = ∫_{0}^{∞}(I − v_{1}a′)(ϕ(s, A)Bϕ(s, A)′  (3.15)
+ e^{sA}D_{v} + D_{v}e^{sA′} − D_{v})(I − av′_{1})e^{−λ_{1}s}λ_{1} ds.
𝔼 V V̂′ = ∫_{0}^{∞} ψ(s, A)(Bϕ(s, A)′ + AD_{v})(I − av′_{1})e^{−λ_{1}s}λ_{1} ds. (3.16)

(ii) if Re λ_{2} = λ_{1}/2 and d := max\{d_{λ} : Re λ = \frac{1}{2}λ_{1}\}, then

(n ln^{2d+1} n)^{−1/2}(X_{n} − nλ_{1}v_{1}, N_{n} − nν) →^{d} (V, V̂),

where (V, V̂) is a vector-valued Gaussian random variable with mean 0. If, for simplicity, each a_{i} = 1, then V = AV̂ and

𝔼(V̂V̂′) = λ_{1}^{−2d}T̂Σ_{II,d}T̂′,

with T̂ := (I − v_{1}a′) ∑_{λ∈Λ_{II}} λ^{−1}P_{λ}.

In the case of random ξ_{i}, the dummy ball method also shows asymptotic normality when λ_{2} ≤ λ_{1}/2 of the number of draws of a ball of type i leading to a specific set of balls being added. We leave the details to the reader.

Functional limit theorems for stopped processes and urns.

Theorem 3.22. Assume (A1)–(A6) and let b ∈ ℝ^{q} with b · v_{1} > 0. Conditioned on essential non-extinction we have as z → ∞, with joint convergence,

(i) z^{−1/2}P_{I}X (τ_{b}(xz)) →^{d} (b · v_{1})^{−1/2}V_{I}(x) in D[0, ∞);

(ii) for every λ ∈ Λ_{II} and k = 0, …, d_{λ}, in D[0, ∞),

(ln z)^{−k−1/2}(z^{x})^{−1/2−i Im λ/λ_{1}}Q_{λ,k}P_{λ}X (τ_{b}(z^{x})) →^{d} (b·v_{1})^{−1/2}λ_{1}^{−(k+1/2)}U_{λ,k}(x).

Here V_{I} and U_{λ,k} are continuous Gaussian vector-valued stochastic processes, defined on [0, ∞). U_{λ,k}(x) is as in Theorem 3.1, V_{I}(0) = 0, 𝔼 V_{I}(x) = 0 and

𝔼 V_{I}(x)V_{I}(y)′ = xΣ_{I}(\frac{y}{x})^{λ_{1}^{−1}A′},  0 < x ≤ y.  (3.17)

The process V_{I} and the families {U_{λ,k}}_{0≤k≤d_{λ}} for different λ ∈ Λ_{II} with Im λ ≥ 0 are independent.

The same results hold for X_{τ_{b}(xz)} and X_{τ_{b}(z^{x})} in the urn process.

Corollary 3.23. Assume (A1)–(A6) and let b ∈ ℝ^{q} with b · v_{1} > 0. Let z → ∞ and condition on essential non-extinction.

(i) If Re λ_{2} < \frac{1}{2}λ_{1}, then, in D[0, ∞),

z^{−1/2}(X (τ_{b}(xz)) − \frac{xz}{b · v_{1}} v_{1}) →^{d} V_{b}(x) := (b · v_{1})^{−1/2}(I − \frac{v_{1}b′}{b · v_{1}} )V_{I}(x).

(ii) If λ_{2} = \frac{1}{2}λ_{1} and Λ_{II} = \{λ_{2}\}, then, with d := d_{λ_{2}}, in D[0, ∞),

(ln z)^{−d−1/2}z^{−x/2}(X (τ_{b}(z^{x})) − \frac{z^{x}}{b · v_{1}} v_{1}) →^{d}

V_{b}(x) := λ_{1}^{−d−1/2}(b · v_{1})^{−1/2}(I − \frac{v_{1}b′}{b · v_{1}} )U_{λ_{2},d}(x).

The limit processes are Gaussian, vanish at 0, and have means 0 and covariances given by (3.17) for (i) and, using (3.3), for (ii)

𝔼 U_{λ_{2},d}(x)U_{λ_{2},d}(y)′ = c(d, d, x, y)(2d + 1)d!^{2}Σ_{II,d},  0 ≤ x ≤ y. (3.18)

The same results holds for the urn process.

With dummy balls as above, Corollary 3.23 leads to the corresponding result for urns.

Theorem 3.24. Assume (A1)–(A6). Let n → ∞ and condition on essential non-extinction.

(i) If Re λ_{2} < \frac{1}{2}λ_{1}, then, in D[0, ∞),

n^{−1/2}(X_{⌊xn⌋} − xnλ_{1}v_{1}) →^{d} V (x),

where V (x) is a continuous Gaussian vector-valued process with V (0) = 0, mean 𝔼 V (x) = 0 and, for 0 < x ≤ y,

𝔼 V (x)V′(y) = ∫_{−λ_{1}^{−1} ln x}^{∞} ψ(s + λ_{1}^{−1} ln x, A)Bψ(s + λ_{1}^{−1} ln y, A)′e^{−λ_{1}s}λ_{1} ds

− xλ_{1}^{2}v_{1}v′_{1}. (3.19)

(ii) If λ_{2} = \frac{1}{2}λ_{1} and Λ_{II} = \{λ_{2}\}, then, with d := d_{λ_{2}}, in D[0, ∞),

(ln n)^{−d−1/2}n^{−x/2}(X_{⌊n^{x}⌋} − n^{x}λ_{1}v_{1}) →^{d} V (x),

where V (x) is a continuous Gaussian vector-valued process with V (0) = 0, mean 𝔼 V (x) = 0 and, for 0 < x ≤ y,

𝔼 V (x)V (y)′ = c̃(d, x, y)λ_{1}^{−2d}(I − 2v_{1}a′P_{λ_{2}})Σ_{II,d}(I − 2P′_{λ_{2}}av′_{1}), (3.20)

where

c̃(d, x, y) := (2d + 1)d!^{2}c(d, d, x, y) = ∑_{j=0}^{d} \frac{2d + 1}{2d + 1 − j} \binom{d}{j}(y − x)^{j}x^{2d+1−j}.

Remark 3.25. If Re λ_{2} = \frac{1}{2}λ_{1} and Im λ_{2} ≠ 0 we obtain a more complicated behaviour. If, for simplicity, λ_{2} and \overline{λ_{2}} are the only eigenvalues with real part \frac{1}{2}λ_{1}, it follows from Theorem 3.22 that (X_{⌊n^{x}⌋} − n^{x}λ_{1}v_{1}) will, asymptotically, oscillate deterministically as a sine function with frequency (Im λ_{2}/2πλ_{1}) ln n, but with amplitude and phase drifting stochastically at a slower rate.

4. Some remarks, extensions and open problems

Remark 4.1. We assume throughout that the set of types is finite. However, in Examples 7.6 and 7.7 we consider two applications where the natural urn models have an infinite number of types (in this case ℕ); luckily it is possible in those applications to consider only a finite number of types at a time.

Another example where the results extend to an infinite space of types (in that case a compact group, for example 𝕋) is described in Example 7.10.

These examples suggest the possibility of (and desire for) an extension of the results in this paper to infinite sets of types (with suitable assumptions). Our matrix A would then be replaced by an operator acting in a suitable space, such as ℓ^{1}(ℕ) or L^{2}(𝕋). It is far from clear how such an extension should be formulated, and we have not pursued this.

Remark 4.2. In several applications, the assumptions (1.1), (1.2) are too restrictive; we want to allow the possibility of removing other balls than the drawn one. Several authors, following Bagchi and Pal [10], have studied the so-called tenable urn models where ξ_{ii} may be an arbitrary negative integer −d_{i}, but it is assumed that d_{i}|ξ_{ji} for all j and d_{i}|X_{0i}; hence X_{ni} is always a multiple of d_{i} and we can never be required to remove balls that do not exist. (We still assume ξ_{ji} ≥ 0 when j ≠ i. We let d_{i} = 1 if ξ_{ii} ≥ 0.)

Note that the corresponding continuous time process X (t) is well-defined, but (if some d_{i} ≥ 2) it is not a branching process because the balls do not evolve independently.

Tenable urns can, nevertheless, easily be reduced to the Athreya–Karlin setting, and thus studied by the results above, by replacing the balls with “superballs”; each superball of type i being equivalent to d_{i} ordinary balls and having activity d_{i}a_{i}. If D denotes the diagonal matrix with D_{ii} = d_{i}, this means that we consider the process X̃_{n} := D^{−1}X_{n}, which is of the type treated above with ξ̃_{i} = D^{−1}ξ_{i} and Ã = D^{−1}AD. Note that Ã and A have the same eigenvalues. If (A2)–(A6) hold, we can thus apply Theorem 3.1 and its consequences in Section 3 to the superball process. Returning to the original process (by multiplying with D), it can easily be verified that all results in Section 3 (including the variance formulas) hold for tenable urns too. See Example 7.5 for an example.

As mentioned in the introduction, the results can be extended even further. Let ξ_{ij} be arbitrary integers (or even real numbers), but assume that they are such that the urn process never can require the removal of balls that do not exist. (An example, from [39], is given in Example 7.8.) The corresponding continuous time process then is well-defined; it is a generalized branching process where the death of one ball may force the removal of others. In this case, we cannot use the Perron–Frobenius theory, so we add the assumption that A has a real eigenvalue λ_{1} > 0 with Re λ < λ_{1} for every other eigenvalue λ, and that there exist corresponding left and right eigenvectors u_{1} and v_{1} such that v_{1i} > 0 for every i and u_{1i} > 0 if i is dominating while u_{1i} = 0 otherwise; we also assume (A2)–(A6). Finally, we assume that Lemma 9.7(iii) holds (we have not been able to prove this in the present generality), for example because ℙ(W = 0) = 0. It may then be verified that the proofs in Sections 9 and 10 hold without modification; hence all results in Section 3 hold for such urn models and generalized branching processes too. (We conjecture that the results hold also if it is possible that the process stops by requiring some X_{ni} to become negative, provided that we condition on this not happening, but we have not pursued this.)

Remark 4.3. Assume for simplicity that essential extinction is impossible, so W > 0. By Theorem 3.1 or Corollary 3.8, the different sets of projections of X (t) in parts (i), (ii) and (iii) of these results, divided by the normalizing factor W^{1/2}, form three asymptotically independent families. This may seem surprising at first sight, but is explained by the three families being essentially determined by what happens in the end, the middle, and the beginning of the process, respectively. More formally, let ε(t) → 0 and ω(t) → ∞ as t → ∞. By Theorem 3.1 (and its proof), dividing everything by W^{1/2} and ignoring terms that are asymptotically negligible, the projections in (i) depend only on the random splits after time t−ω(t), the projections in (ii) depend only on the random splits in [ε(t)t, (1−ε(t))t], and the projections in (iii) depend only on the random splits before ω(t).

For the urn process in Theorems 3.17–3.19 we find, by the exponential relation between n and t, that the variables asymptotically depend only on the draws after εn, the draws in [n^{ε}, n^{1−ε}] and the draws before ω, respectively. The same holds if we consider the different cases together as in Remark 3.20.

Remark 4.4. It would be interesting to extend the results to cases where (A4) or (A6) does not hold. A typical case is when A is triangular: if, for example, q = 2 and A = ( \begin{matrix}α & 0 \\ β & δ \end{matrix} ) with β > 0, the conditions hold if α > δ but not otherwise. It might be possible to handle this case by combining the methods here with the ones in [34], where a detailed study is made in the related case of multitype Galton–Watson processes in discrete time. Some new phenomena will arise, however, see [34], [20] and [35].

Remark 4.5. Mahmoud [40] has initiated the study of urn models where several, say 2, balls are drawn at the same time, and balls are added depending on the drawn combination of types. It may be possible to study such models too by the methods of this paper, first considering the corresponding continuous time model, but we have not pursued this. (This case is substantially more complicated than the standard case treated here; for example, the continuous time model will explode in finite time.)

Remark 4.6. Asmussen [3] has proved laws of iterated logarithm for X (t) in the irreducible case. (The results in the cases Re λ_{2} < \frac{1}{2} Re λ_{1} and Re λ_{2} = \frac{1}{2} Re λ_{1} are different.) By the Athreya–Karlin embedding, this yields laws of iterated logarithm for the urn process X_{n} complementing the results above; we leave the details to the reader. Such laws for the urn process have been proved in a special case by Bai, Hu and Zhang [12].

Remark 4.7. Our methods give no information on the rate of convergence. Using other methods, Hwang [23] has found the rate of convergence to the limiting normal distribution for a specific variable in Example 7.9 below: the rate is n^{−γ} where γ = min(\frac{1}{2}, 3(\frac{1}{2} − Re λ_{2}/λ_{1})). It is tempting to guess that this might hold rather generally.

5. Variance calculations

In several of the theorems in Section 3, the variances and covariances of the limits are given as integrals of matrix functions. In any specific application, these integrals can be evaluated by first transforming the matrix A to Jordan normal form. (A computer algebra package is helpful, and can do the integration directly if q is not too large.) We will here give some simplifications in important special cases. We concentrate on the main cases; the reader may add further similar results. See the examples in Section 7 and for various applications of these results. See also [29] for more complicated applications.

First we consider the case when the replacement vectors ξ_{i} are deterministic. Let D be the diagonal q × q matrix with entries

D_{ii} := \begin{cases}v_{1i}/a_{i}, & a_{i} ≠ 0, \\ 0, & a_{i} = 0. \end{cases}  (5.1)

Lemma 5.1. (i) If each ξ_{i} is deterministic, then B = ADA′.

(ii) If furthermore Re λ_{2} < \frac{1}{2}λ_{1}, then the covariance matrix Σ in Theorem 3.17 is given by

Σ = ∫_{0}^{∞}(A − λ_{1}v_{1}a′)e^{sA}De^{sA′}(A − λ_{1}v_{1}a′)′e^{−λ_{1}s}λ_{1} ds.  (5.2)

Equivalently, if g_{i}(s) := (A − λ_{1}v_{1}a′)e^{sA}δ_{i}, where (δ_{i})_{j} = δ_{ij}, then

Σ = ∫_{0}^{∞} ∑_{a_{i}≠0} a_{i}^{−1}v_{1i}g_{i}(s)g_{i}(s)′e^{−λ_{1}s}λ_{1} ds.  (5.3)

If further u_{1} = a, then g_{i}(s) = Ae^{sA}(I − v_{1}u′_{1})δ_{i} = Ae^{sA}(δ_{i} − u_{1i}v_{1}).

Proof. We now have, by (2.13) and (2.1), when a_{i} ≠ 0,

(B_{i})_{jk} = 𝔼 ξ_{ij} 𝔼 ξ_{ik} = a_{i}^{−2}A_{ji}A_{ki},

and then by (2.14)

B_{jk} = ∑_{i=1}^{q} a_{i}v_{1i}(B_{i})_{jk} = ∑_{a_{i}≠0} a_{i}^{−1}v_{1i}A_{ji}A_{ki} = (ADA′)_{jk},

which proves (i).

For (ii), we have by (3.11) and (3.12) ϕ(s, A)A = e^{sA} − 1 and

ψ(s, A)A = e^{sA}A − λ_{1}v_{1}a′(e^{sA} − 1) = (A − λ_{1}v_{1}a′)e^{sA} + λ_{1}v_{1}a′.

Hence (i) now yields

ψ(s, A)Bψ(s, A)′ = (A − λ_{1}v_{1}a′)e^{sA}De^{sA′}(A − λ_{1}v_{1}a′)′

+λ_{1}v_{1}a′De^{sA′}(A − λ_{1}v_{1}a′)′ + λ_{1}(A − λ_{1}v_{1}a′)e^{sA}Dav′_{1}

+λ_{1}^{2}v_{1}a′Dav′_{1}.

It follows from (5.1) that (Da)_{i} = v_{1i} when a_{i} ≠ 0, and thus a′Da = a′v_{1} = 1 and ADa = Av_{1} = λ_{1}v_{1}. It follows that (A − λ_{1}v_{1}a′)e^{sA}(Da − v_{1}) = 0, and thus (A − λ_{1}v_{1}a′)e^{sA}Da = (A − λ_{1}v_{1}a′)e^{sA}v_{1} = 0. Hence, the second and third terms in the sum vanish. The fourth equals λ_{1}^{2}v_{1}v′_{1}, and thus (5.2) follows from (3.13). Since D = ∑_{i} D_{ii}δ_{i}δ′_{i}, (5.3) follows from (5.2).

In the special case u_{1} = a, we have (A − λ_{1}v_{1}a′)e^{sA} = Ae^{sA}(I − v_{1}u′_{1}), and the alternative formulas for g_{i} follow. □

Remark 5.2. In general, the argument above shows that if B_{i}^{∘} is the covariance matrix 𝔼(ξ_{i} − 𝔼 ξ_{i})(ξ_{i} − 𝔼 ξ_{i})′ and B^{∘} := ∑_{i=1}^{q} a_{i}v_{1i}B_{i}^{∘}, then

\begin{aligned}Σ & = ∫_{0}^{∞} ψ(s, A)B^{∘}ψ(s, A)′e^{−λ_{1}s}λ_{1} ds \\ & + ∫_{0}^{∞}(A − λ_{1}v_{1}a′)e^{sA}De^{sA′}(A − λ_{1}v_{1}a′)′e^{−λ_{1}s}λ_{1} ds, \end{aligned}

which separates the contributions to the asymptotic variance coming from the randomness in the ξ_{i} and the randomness in the draws.

Another simplifying case is when A is diagonalizable. In that case, there exist dual bases {u_{i}}_{i=1}^{q} and {v_{i}}_{i=1}^{q} of left and right eigenvectors of A, i.e. vectors such that u_{i}A = λ_{i}u_{i}, Av_{i} = λ_{i}v_{i} and u_{i} · v_{j} = δ_{ij} (where the λ_{i}, i = 1, …, q do not have to be distinct; we assume that the bases are ordered such that λ_{1} ≥ Re λ_{2} ≥ … as elsewhere).

Lemma 5.3. (i) If A is diagonalizable and {u_{i}}_{i=1}^{q} and {v_{i}}_{i=1}^{q} are dual bases of eigenvectors, then, with L_{I} := {i : λ_{i} ∈ Λ_{I}} and L_{II} := {i : λ_{i} ∈ Λ_{II}},

Σ_{I} = ∑_{j,k∈L_{I}} \frac{u′_{j}Bu_{k}}{λ_{1} − λ_{j} − λ_{k}} v_{j}v′_{k}  and  Σ_{II} = ∑_{j∈L_{II}}(u′_{j}B\overline{u_{j}})v_{j}v_{j}^{∗}.

(ii) If further each ξ_{i} is deterministic, then

Σ_{I} = ∑_{j,k∈L_{I}} \frac{λ_{j}λ_{k}u′_{j}Du_{k}}{λ_{1} − λ_{j} − λ_{k}} v_{j}v′_{k}  and  Σ_{II} = ∑_{j∈L_{II}} |λ_{j}|^{2}(u′_{j}D\overline{u_{j}})v_{j}v_{j}^{∗}.

(iii) If further Re λ_{2} < \frac{1}{2}λ_{1}, then the covariance matrix Σ in Theorem 3.17 is given by, with w_{j} := λ_{j}v_{j} − λ_{1}(a · v_{j})v_{1},

Σ = ∑_{j,k=2}^{q} \frac{λ_{1} u′_{j}Du_{k}}{λ_{1} − λ_{j} − λ_{k}} w_{j}w′_{k}.

(iv) If the assumptions in (i) and (ii) hold and Re λ_{2} = \frac{1}{2}λ_{1}, then the covariance matrix Σ in Theorem 3.18 is given by, with w_{j} as in (iii),

Σ = ∑_{j∈Λ_{II}}( u′_{j}Du_{j})w_{j}w_{j}^{∗}.

Proof. We have P_{I} = ∑_{j∈L_{I}} v_{j}u′_{j} and P_{I}e^{sA} = ∑_{j∈L_{I}} e^{sλ_{j}}v_{j}u′_{j}. Hence the first equality in (i) follows from the definition (2.15). The second follows similarly from (2.16).

For (ii) we use (i) and B = ADA′ from Lemma 5.1, recalling that u′_{j}A = λ_{j}u′_{j}.

In case (iii), when L_{I} = {2, …, q}, the result follows from (5.2) and

(A − λ_{1}v_{1}a′)e^{sA} = ∑_{j≥1}(A − λ_{1}v_{1}a′)e^{sA}v_{j}u′_{j} = ∑_{j≥2} w_{j}e^{sλ_{j}}u′_{j}.

Finally, since A is diagonalizable, d = 0 and (iv) follows from (3.14) and (ii). □

A third common simplifying case is when a fixed number of balls is added each time, i.e. each ∑_{j} ξ_{ij} = m is deterministic and independent of i, and further each a_{i} = 1. We state the result somewhat more generally.

Lemma 5.4. Suppose that a · 𝔼 ξ_{i} = m for some m > 0 and every i. (For example, this holds if exactly m balls are added each time, and each a_{i} = 1.) Then λ_{1} = m and u_{1} = a. The covariance matrix (3.13) in Theorem 3.17 equals mΣ_{I}. The covariance matrix (3.14) in Theorem 3.18 equals λ_{1}^{−2d}Σ_{II,d}.

Proof. We have, by (2.1) and assumption,

(a′A)_{k} = ∑_{j=1}^{q} a_{j}A_{jk} = ∑_{j=1}^{q} a_{j}a_{k} 𝔼 ξ_{kj} = a_{k}a · 𝔼 ξ_{k} = ma_{k}

so a′A = ma′ and a is a non-negative left eigenvector of A. This implies that m = λ_{1} and that a is a multiple of u_{1}; by our normalizations (2.2) and (2.3), a = u_{1}.

It follows that a′P_{λ_{1}} = a′ and a′P_{λ} = a′P_{λ_{1}}P_{λ} = 0 when λ ≠ λ_{1}; thus T = 0 in Theorem 3.18, which shows the claim about (3.14).

For (3.13), we may use algebraic manipulations as in the lemmas above, but it seems easier to proceed as follows. Conditioning on X_{0}, we may assume that X_{0} is fixed. Then a · X_{n} = mn + l, where l := a · X_{0}. Thus we can obtain X_{n} by stopping at τ_{b}(n + l/m), where b = m^{−1}a. We thus obtain the conclusion of Theorem 3.17 directly from Corollary 3.11, with Σ given by (3.7). Furthermore, since b′P_{λ} = m^{−1}a′P_{λ} = 0 for every λ ≠ λ_{1}, we have b′P_{I} = 0 and P′_{I}b = 0. Moreover, b · v_{1} = m^{−1}a · v_{1} = m^{−1}. Since Σ_{I} = P_{I}Σ_{I}P′_{I} by (2.15), b′Σ_{I} = 0 and Σ_{I}b = 0, and the middle expression in (3.7) equals mΣ_{I}. □

If further m = 1, we have in addition the following.

Lemma 5.5. If each |ξ_{i}| = 1, i.e. exactly one ball is added each time, and each a_{i} = 1, then B = D.

Proof. Since only one ξ_{ij} is non-zero at a time, B_{i} in (2.13) is diagonal, with (B_{i})_{jj} = ℙ(ξ_{ij} = 1) = 𝔼 ξ_{ij}. Thus B is diagonal with, see (2.18),

B_{jj} = ∑_{i=1}^{q} v_{1i}a_{i} 𝔼 ξ_{ij} = (λ_{1}v_{1})_{j} = v_{1j}.  □

We also give a calculation of the variances and covariances in Theorem 3.21.

Lemma 5.6. Suppose that a_{i} = 1 and 𝔼 ∑_{j} ξ_{ij} = m for every i, and that A is diagonalizable with dual bases of eigenvectors {u_{i}}_{i=1}^{q} and {v_{i}}_{i=1}^{q}. Then λ_{1} = m and, if Re λ_{2} < λ_{1}/2, the limit in Theorem 3.21 satisfies

𝔼 V̂V̂′ = ∑_{j,k=2}^{q}((\frac{1}{λ_{1} − λ_{j}} + \frac{1}{λ_{1} − λ_{k}} )\frac{u′_{j}Bu_{k}}{λ_{1} − λ_{j} − λ_{k}}

+ \frac{λ_{1}^{2} − λ_{j}λ_{k}}{(λ_{1} − λ_{j})(λ_{1} − λ_{k})}u′_{j}Du_{k})v_{j}v′_{k} (5.4)

𝔼 V V̂′ = ∑_{j,k=2}^{q}(\frac{λ_{1}}{(λ_{1} − λ_{j})(λ_{1} − λ_{j} − λ_{k})}u′_{j}Bu_{k} + \frac{λ_{1}λ_{j}}{λ_{1} − λ_{j}} u′_{j}Du_{k})v_{j}v′_{k}

+ ∑_{k=2}^{q} \frac{u′_{1}Bu_{k}}{λ_{1} − λ_{k}} v_{1}v′_{k} (5.5)

Proof. By Lemma 5.4, λ_{1} = m and u_{1} = a. Thus, I − v_{1}a′ = I − v_{1}u′_{1} = ∑_{j=2}^{q} v_{j}u′_{j}, and it follows from (3.15) that, with ϕ(s, λ) = ∫_{0}^{s} e^{tλ}dt,

𝔼 V̂V̂′ = ∑_{j,k=2}^{q} v_{j} ∫_{0}^{∞} u′_{j}(ϕ(s, λ_{j})Bϕ(s, λ_{k})

+ e^{sλ_{j}}D + De^{sλ_{k}} − D)u_{k}e^{−λ_{1}s}λ_{1} ds v′_{k}

By the definition of ϕ(s, λ) and changes of order of integration,

∫_{0}^{∞} ϕ(s, λ_{j})ϕ(s, λ_{k})e^{−λ_{1}s}λ_{1} ds = ∫∫∫_{0<t,u<s} e^{λ_{j}t+λ_{k}u−λ_{1}s}λ_{1} ds du dt

= ∫∫_{0<t<u} e^{λ_{j}t+(λ_{k}−λ_{1})u} du dt + ∫∫_{0<u<t} e^{λ_{k}u+(λ_{j}−λ_{1})t} dt du

= \frac{1}{λ_{1} − λ_{k}} \frac{1}{λ_{1} − λ_{k} − λ_{j}} + \frac{1}{λ_{1} − λ_{j}} \frac{1}{λ_{1} − λ_{j} − λ_{k}} .

The integrals of the other terms are easily computed and a summation yields (5.4).

Similarly, since (3.12) and a = u_{1} imply u′_{1}ψ(s, A) = u′_{1} and u′_{j}ψ(s, A) = e^{sλ_{j}}u′_{j} for j ≠ 1, (3.16) yields

𝔼 V V̂′ = ∑_{j,k=2}^{q} v_{j} ∫_{0}^{∞} u′_{j}e^{sλ_{j}}(Bϕ(s, λ_{k}) + λ_{j}D)u_{k}e^{−λ_{1}s}λ_{1} ds v′_{k}

+ ∑_{k=2}^{q} v_{1} ∫_{0}^{∞} u′_{1}(Bϕ(s, λ_{k}) + λ_{1}D)u_{k}e^{−λ_{1}s}λ_{1} ds v′_{k}

Again, the integrals are easily computed, and (5.5) follows; note that u′_{1}Du_{k} = v′_{1}u_{k} = 0 for k ≠ 1. □

Remark 5.7. These results can be extended to the covariances in the process limits. For example, under the assumptions of Lemma 5.3(iii), an argument as in the proofs of Lemmas 5.1 and 5.3 shows that in (3.19)

𝔼 V (x)V (y)′ = ∑_{j,k≥2} x^{1−λ_{k}/λ_{1}}y^{λ_{k}/λ_{1}} \frac{λ_{1} u′_{j}Du_{k}}{λ_{1} − λ_{j} − λ_{k}} w_{j}w′_{k}, 0 ≤ x ≤ y.

Similarly, under the assumptions of Lemma 5.4 and with Re λ_{2} < \frac{1}{2}λ_{1}, we have in (3.19), by Corollary 3.23 with b := m^{−1}a and (3.17),

𝔼 V (x)V (y)′ = mxΣ_{I}(y/x)^{m^{−1}A′},  0 < x ≤ y.  (5.6)

6. Relations to previously known results

Large parts of Theorem 3.1 are known since Athreya’s thesis, at least in the irreducible case: the a.s. convergence is in [4] and [9, Theorem V.7.2], and the limits in Corollary 3.8 are proved in [5], [6], see also [9, §V.8]. Our results give more explicit formulas for the asymptotic variances and covariances, and the extension to stochastic processes. The independence between the limit processes in (i) and (ii) seems new too.

Also Theorems 3.16–3.19 for urn models are basically due to Athreya and Karlin [8], see also [9, §V.9], but it is not evident how to obtain explicit formulas for asymptotic variances from their paper.

One of the purposes of this paper is to draw attention to the embedding method in [8], adding some details and making the results simpler to apply. In our opinion, this method has been neglected for too long. Several authors have, however, derived similar results for Pólya urns, in more or less general situations, either by calculating moments by recursion formulas or by martingale methods. (Of course, the embedding method uses martingales too, for the branching process. It is thus not really a question of using martingales or not; the main difference is rather whether to use discrete time martingales directly or to first randomize the splitting times by the continuous time branching process and then use continuous time martingales.) It seems that the (discrete time) martingale methods works fine when the number of added balls is fixed, but the extra randomization in the embedding method makes it much easier to handle the general case.

Some important papers with general limit theorems for urn processes that contain special cases of our results are the following. (Our description is brief; for exact conditions and results, see the cited papers. See also further references in these papers.)

Bagchi and Pal [10] gave, using the method of moments, limit theorems when q = 2 and the number of added balls is fixed, see Example 7.3 below. Gouet [20] gave functional limit theorems in the same case, using a martingale central limit theorem.

Smythe [51] used martingale methods to establish asymptotic normality of X_{n} and joint asymptotic normality of X_{n} and N_{n} (thus inspiring our Theorem 3.21) when λ_{2} < \frac{1}{2}λ_{1}, allowing removals (in the tenable case) and assuming that the expected number of added balls does not depend on the type of the drawn ball, i.e. that the column sums of A are the same, and thus equal to λ_{1}. (He also assumed some technical simplifications, for example that A is diagonalizable). However, no general variance formulas were given, although some examples are given and it is stated that asymptotic variances and covariances in principle can be computed more generally. (See also the special case in [44] where asymptotic variances are given.) It is noteworthy that only the expected number of added balls is assumed constant in [51]; most papers applying martingale methods use the stronger assumption that the actual number of added balls is constant.

Bai and Hu [11] used similar martingale methods to establish asymptotic normality of X_{n}. They assumed that the number of added balls is constant, usually 1, and that Re λ_{2} ≤ \frac{1}{2}λ_{1} = \frac{1}{2}. A new feature of [11] is that they allow time-dependent transition probabilities (converging to a limit), a case not studied in this paper.

Bai, Hu and Zhang [12] studied the case of two types in more detail, using martingale methods and the Skorohod embedding theorem. Again they allow some time-dependency, and assume that the number of added balls is 1 (with some randomness allowed in the time dependent case), and that Re λ_{2} ≤ \frac{1}{2}λ_{1} = \frac{1}{2}. They results include functional limit theorems, laws of iterated logarithm, and estimates of rates of convergence.

Remark 6.1. A method that has, as far as we know, not yet been used to study this type of urn models is to use a general limit theorem for Markov processes such as [31, Theorem 19.28]. This seems to have the potential of giving comparatively simple proofs of several results in this paper, and could probably be used to attack some of the extensions mentioned in Section 4.

7. Examples and applications to random trees

We give several examples of urn models that illustrate the results above. We concentrate on already studied models and show how several previosly known results follow from our theorems by routine calculations; we encourage the reader to compare the methods. We also give some new results.

Unless otherwise stated, all activities a_{i} = 1, the urn is irreducible, (A1)– (A6) hold, and (essential) extinction is impossible; this can in each case easily be verified using Lemma 2.1. We sometimes omit minor details, such as specifying X_{0}; similarly, when convenient we shift the indices and start the process with X_{1}.

Example 7.1. First a trivial example. If ξ_{1}, …, ξ_{q} are random with the same distribution, the drawn types do not matter and X_{n} − X_{0} is a sum of n i.i.d. random vectors. Thus X_{n} is asymptotically normal by the central limit theorem. In this case, A has rank 1 so λ_{2} = 0 with multiplicity q − 1, and there are no further eigenvalues. It can be verified that Theorem 3.17 indeed yields the normal limit given by the central limit theorem. Similarly, Theorem 3.24 yields the same result as Donsker’s theorem.

If b = δ_{1} = (1, 0, …, 0)′, then X_{τ_{b}(z)} is the vector obtained by summing i.i.d. copies of ξ_{1} until the first component is at least z. Corollary 3.11 yields asymptotic normality, as shown in [22], [21].

Example 7.2 (Friedman’s urn). A classic example is Friedman’s urn [19],

[18] (studied already by Bernstein [13], [14] and Savkevitch [49]), where q = 2, ξ_{1} = (α, β)′ and ξ_{2} = (β, α)′ for some integers α and β. We assume α + β > 0. If β = 0, we have the original Pólya urn [17], [47] which is reducible; (A4), (A5), (A6) fail so our results do not apply, and it is well known that X_{n}/n converges to a Beta distribution instead of a constant [47], [30].

If β > 0, our assumptions hold. We have λ_{1} = α + β and λ_{2} = α − β, so λ_{2} < \frac{1}{2}λ_{1} is equivalent to α < 3β; in this case Theorem 3.17 yields by Lemma 5.3(iii), with v_{1} = \frac{1}{2}(1, 1)′, v_{2} = \frac{1}{2}(1, −1)′, u_{1} = (1, 1)′, u_{2} = (1, −1)′,

n^{−1/2} (X_{n} − n\frac{α + β}{2} (\begin{matrix}1 \\ 1 \end{matrix})) →^{d} N (0, \frac{(α + β)(α − β)^{2}}{4(3β − α)} (\begin{matrix}1 & −1 \\ −1 & 1 \end{matrix}))

as shown by Bernstein [13], [14] and Freedman [18].

Similarly, if α = 3β > 0, Theorem 3.18 and Lemma 5.3(iv) yield [13, 18]

(n ln n)^{−1/2} (X_{n} − n\frac{α + β}{2} (\begin{matrix}1 \\ 1 \end{matrix})) →^{d} N (0, \frac{(α − β)^{2}}{4} (\begin{matrix}1 & −1 \\ −1 & 1 \end{matrix})) .

Example 7.3 (general 2-type urn). More generally, consider the case q = 2 with nonrandom ξ_{1} = (α, β)′ and ξ_{2} = (γ, δ)′. We assume that β, γ > 0 so that the urn is irreducible, and that λ_{1} > 0. If α, δ ≥ −1, (A1)–(A6) hold and extinction is impossible; by Remark 4.2 we can also allow other negative values for α and δ under suitable conditions, for example in the tenable case.

The case α + β = γ + δ has been studied by several authors, including Bernstein [13], [14] and (also for the tenable case) Bagchi and Pal [10], who proved asymptotic normality when λ_{2} ≤ λ_{1}/2. Gouet [20] gave functional limit theorems.

We extend their results as follows. We write, for notational convenience, κ := (α + δ)/2, ε := (α − δ)/2, ρ := \sqrt{ε^{2} + βγ} > 0. Thus

A = (\begin{matrix}α & γ \\ β & δ \end{matrix}) = (\begin{matrix}κ + ε & γ \\ β & κ − ε \end{matrix}) .

Simple calculations yield the eigenvalues λ_{±} = κ ± ρ (the indices +, − are more convenient than 1, 2 in this example) and the dual bases of eigenvectors (v_{+}, v_{−}) and (u_{+}, u_{−}) with

v_{±} = c_{±}^{−1} (\begin{matrix}ε ± ρ \\ β \end{matrix}) ,  u_{±} = \frac{c_{±}}{(ε ± ρ)^{2} + βγ} (\begin{matrix}ε ± ρ \\ γ \end{matrix}) ,

where c_{±} are normalization constants. We take c_{+} = β + ρ + ε so that a · v_{+} = 1, and choose c_{−} = 1.

If λ_{−} < \frac{1}{2}λ_{+}, or, equivalently, κ < 3ρ, Theorem 3.17 yields n^{−1/2}(X_{n} −

n(κ + ρ)v_{+}) →^{d} N(0, Σ), where Lemma 5.3(iii) yields

Σ = \frac{λ_{+}}{λ_{+} − 2λ_{−}} (u′_{−}Du_{−})ww′,

with w := λ_{−}v_{−} −λ_{+}(a·v_{−})v_{+}. By elementary calculations, using (ρ−ε)^{2} + βγ = 2ρ(ρ − ε) and (ρ − ε)(β + ρ + ε) = β(γ + ρ − ε), we find

u′_{−}Du_{−} = \frac{γ}{4ρ^{2}(ρ − ε)}  and  w = \frac{2ρβ}{β + ρ + ε} (\begin{matrix}γ − α \\ δ − β \end{matrix})

and

Σ = \frac{(κ + ρ)βγ}{(3ρ − κ)(β + ρ + ε)(γ + ρ − ε)} (\begin{matrix}γ − α \\ δ − β \end{matrix}) (γ − α, δ − β) . (7.1)

(It is easy to see that X_{n1} and X_{n2} are linearly dependent, which explains why Σ has rank 1. It is thus sufficient to study only X_{n1}.)

In the special case α + β = γ + δ = m studied in [10] this simplifies: now λ_{+} = κ + ρ = m, ρ = (β + γ)/2 and 3ρ − κ = m + 2ρ − 2κ = m + 2(γ − α), so

Σ = \frac{mβγ}{(m + 2(γ − α))(β + γ)^{2}} (γ − α)^{2} (\begin{matrix}1 & −1 \\ −1 & 1 \end{matrix})

in accordance with [10].

When κ = 3ρ, λ_{−} = \frac{1}{2}λ_{+} and Theorem 3.18 yields (n ln n)^{−1/2}(X_{n} −

n(κ+ρ)v_{+}) →^{d} N(0, Σ), where by Lemma 5.3(iv) and the calculations above

Σ = \frac{βγ}{(β + ρ + ε)(γ + ρ − ε)} (\begin{matrix}γ − α \\ δ − β \end{matrix}) (γ − α, δ − β) .  (7.2)

In the special case α + β = γ + δ (when α = β + 2γ, δ = 2β + γ) this simplifies to βγ ( \begin{matrix}1 & −1 \\ −1 & 1 \end{matrix} ) as given by [10].

Theorem 3.24 yields functional convergence in D[0, ∞) to a Gaussian process when κ ≤ 3ρ. (In the special case α + β = γ + δ, this was proved by Gouet [20].) If κ < 3ρ, then

n^{−1/2}(X_{⌊xn⌋} − xn(κ + ρ)v_{+}) →^{d} V (x),

where 𝔼 V (x) = 0 and, by Remark 5.7, with Σ given in (7.1),

𝔼 V (x)V′(y) = x^{1−λ_{−}/λ_{+}}y^{λ_{−}/λ_{+}}Σ, 0 ≤ x ≤ y.

If κ = 3ρ, then, instead,

(ln n)^{−1/2}n^{−x/2}(X_{⌊n^{x}⌋} − n^{x}(κ + ρ)v_{+}) →^{d} V (x),

where 𝔼 V (x) = 0 and, with Σ given in (7.2),

𝔼 V (x)V′(y) = xΣ, 0 ≤ x ≤ y.

Example 7.4 (randomized play-the-winner). The calculations in Example 7.3 can be extended to random ξ_{1} and ξ_{2}. We consider for simplicity only the case when exactly one ball is added each time. This gives the randomized play-the-winner rule for clinical trials introduced by Wei and Durham [53], see also [54]: If the drawn ball has type i (i = 1 or 2), we add a ball with the same type with probability p_{i} and a ball with the opposite type with probability q_{i} = 1 − p_{i}. Here p_{1} and p_{2} are given numbers with 0 ≤ p_{i} < 1. We have

A = (\begin{matrix}p_{1} & q_{2} \\ q_{1} & p_{2} \end{matrix}) .

It is easily seen that, see Lemma 5.4, λ_{1} = 1, λ_{2} = p_{1} + p_{2} − 1 = p_{1} − q_{2},

u_{1} = (\begin{matrix}1 \\ 1 \end{matrix}) , v_{1} = \frac{1}{q_{1} + q_{2}} (\begin{matrix}q_{2} \\ q_{1} \end{matrix}) , u_{2} = (\begin{matrix}q_{1} \\ −q_{2} \end{matrix}) , v_{2} = \frac{1}{q_{1} + q_{2}} (\begin{matrix}1 \\ −1 \end{matrix}) .

Lemma 5.5 yields

B = D = \frac{1}{q_{1} + q_{2}} (\begin{matrix}q_{2} & 0 \\ 0 & q_{1} \end{matrix}) .

When λ_{2} < 1/2, Theorem 3.17 applies with, by Lemmas 5.4 and 5.3(i),

Σ = Σ_{I} = \frac{u′_{2}Bu_{2}}{1 − 2λ_{2}} v_{2}v′_{2} = \frac{q_{1}q_{2}}{(1 − 2λ_{2})(q_{1} + q_{2})^{2}} (\begin{matrix}1 & −1 \\ −1 & 1 \end{matrix}) .

(Note that q_{1} + q_{2} = 1 − λ_{2}.) Moreover, by Theorem 3.21(i), we have joint asymptotic normality of the numbers of balls of different types after n draws and the numbers of drawn balls of different types. Considering, as we may, only type 1, we have

n^{−1/2}(X_{n1} − nq_{2}/(q_{1} + q_{2}), N_{n1} − nq_{2}/(q_{1} + q_{2})) →^{d} N(0, Σ_{1}),

where Lemma 5.6 easily gives

Σ_{1} = \frac{q_{1}q_{2}}{(1 − 2λ_{2})(q_{1} + q_{2})^{2}} (\begin{matrix}1 & 1 + 2λ_{2} \\ 1 + 2λ_{2} & 3 + 2λ_{2} \end{matrix}) .

Similarly, when λ_{2} = 1/2, Theorem 3.18 applies with d = 0 and, by Lemmas 5.4 and 5.3(i),

Σ = Σ_{II,d} = Σ_{II} = \frac{q_{1}q_{2}}{(q_{1} + q_{2})^{2}} (\begin{matrix}1 & −1 \\ −1 & 1 \end{matrix}) = 4q_{1}q_{2} (\begin{matrix}1 & −1 \\ −1 & 1 \end{matrix}) .

Furthermore, Theorem 3.21(ii) applies, and V̂ = A^{−1}V = 2V.

For earlier proofs of these results and some extensions, see [52], [51], [11], [12].

Example 7.5 (random 2-3 trees). Bagchi and Pal [10] applied their general result to random 2-3 trees. In such trees, all internal nodes have 2 or 3 children, and all external nodes (leaves) are at the same distance from the root. Keys are associated either with the leaves or with the internal nodes [10], [2]. We define, following [10], the type of an internal node to be W if it has 2 children and B if it has 3; an external node has the same type as its parent. When the tree is grown randomly, a new external node is inserted adjacent to a randomly chosen old one; this either transforms a W -type node at the lowest level to B, or splits a B-type node at the lowest level into two nodes of type W (possibly inducing further splits higher up). Bagchi and Pal [10] study the types of the external nodes as an urn process, with ξ_{1} = (−2, 3)′ and ξ_{2} = (4, −3)′. This does not satisfy (1.2) so, as noted in [10], the Athreya–Karlin embedding is not immediately available. However, this example is of the “tenable” type described in Remark 4.2, where we can use the superball trick. In this case, the superball method is very natural; it means that we consider the lowest level internal nodes instead of the external nodes. This “internal” urn model for 2-3 trees was one of the examples considered by Aldous, Flannery and Palacios [2], who noted that the Athreya–Karlin embedding works.

For the internal version, we have the urn model with ξ_{1} = (−1, 1)′ and ξ_{2} = (2, −1)′ and the activities a_{1} = 2, a_{2} = 3; hence A = ( \begin{matrix}−2 & 6 \\ 2 & −3 \end{matrix} ). Straightforward calculations yield λ_{1} = 1, λ_{2} = −6 and, by Theorem 3.17 and Lemma 5.3(iii), for example n^{−1/2}(W_{n} − \frac{2}{7}n) →^{d} N(0, \frac{108}{637}), if W_{n} is the number of W -type internal nodes in the lowest level when we have n nodes. For w_{n}, the number of W -type external nodes, we have w_{n} = 2W_{n} and thus

n^{−1/2}(w_{n} − \frac{4}{7}n) →^{d} N(0, \frac{432}{637}),  (7.3)

as shown by other methods in [10].

Alternatively, we can obtain this directly by using Theorem 3.17 on the external version in [10], although (A1) is not satisfied; as remarked in Remark 4.2, this is allowed for tenable urns (and for some other urns too). In this direct approach, A = ( \begin{matrix}−2 & 4 \\ 3 & −3 \end{matrix} ), and (7.3) follows by simple computations using Lemma 5.3(iii), or directly by (7.1). Note that, as always with the superball trick, A differs for the two versions, but the eigenvalues are the same, see Remark 4.2.

Example 7.6 (random recursive trees). Mahmoud and Smythe [43] used a generalized Pólya urn to study random recursive trees and obtained the asymptotic normal distribution of the nodes of outdegrees 0, 1 and 2. They indicated that the results in principle extend to higher degrees; we can now do this.

The distribution of outdegrees is the same as the distribution of types in a generalized Pólya urn with infinitely many types {0, 1, 2, … } and the rule that if a ball with type i is drawn, it is removed and replaced by a ball of type i + 1 and a ball of type 0, see [43]. In our notation, ξ_{ij} = −δ_{ij} + δ_{0j} + δ_{i+1,j}. Our theorems assume that the number of types is finite, but luckily we can in this application truncate and lump all high degrees together. Thus, let M ≥ 1 be an integer and use the types {0, 1, …, M} only (thus q = M + 1), where now type M represents all outdegrees ≥ M. The replacement vectors ξ_{i} are as in the infinite model when i < M, while now ξ_{M} = (1, 0, …, 0)′.

Exactly one ball is added each time, so λ_{1} = 1 and u_{1} = a = (1, 1, …, 1)′ by Lemma 5.4. It is easily verified that v_{1} = (1/2, 1/4, …, 2^{−M}, 2^{−M})′, i.e. v_{1i} = 2^{−i−1} for 0 ≤ i < M and v_{1M} = 2^{−M}. In particular, Theorem 3.16 shows that X_{ni}/n →^{a.s.} 2^{−i−1} for every i ≥ 0 (by taking M > i); the weaker statement X_{ni}/n →^{p} 2^{−i−1} was shown by Meir and Moon [45].

It can be shown, see [29], that A has besides λ_{1} = 1 only the eigenvalue −1 (with multiplicity q − 1 = M). (Moreover, d_{−1} = M − 1, so A is not diagonalizable when M ≥ 2.) Since thus λ_{2} = −1, Theorem 3.17 applies for every M, and the vector (n^{−1/2}(X_{ni}−2^{−i−1}n))_{i=0}^{M−1} converges in distribution to a Gaussian vector. Since M is arbitrary, this is the same as convergence of the infinite vector (n^{−1/2}(X_{ni} − 2^{−i−1}n))_{i=0}^{∞} in ℝ^{∞}, see [15, p. 19]. In other

words, n^{−1/2}(X_{ni} − 2^{−i−1}n) →^{d} V_{i}, jointly for all i ≥ 0, as n → ∞, where the V_{i} are jointly Gaussian variables with means 𝔼 V_{i} = 0. The (co)variances Σ_{jk} := Cov(V_{j}, V_{k}) are calculated in [29] using Lemma 5.1.

Similarly, Theorem 3.24 yields a functional limit theorem: n^{−1/2}(X_{⌊xn⌋,i}−

2^{−i−1}xn) →^{d} V_{i}(x) in D[0, ∞), where the V_{i}(x) are continuous Gaussian processes with 𝔼 V_{i}(x) = 0. Again, see [29] for covariances.

Example 7.7 (random plane recursive trees). Mahmoud, Smythe and Szymański [44] studied random plane recursive trees and obtained (among other results) the asymptotic normal distribution of the number of nodes of outdegrees 0, 1 and 2. The outdegrees can be modelled using a generalized Pólya urn with infinitely many colours as in Example 7.6; the ξ_{ij} are the same, but now the activity a_{i} = i + 1. In this case it is advantageous to use the reverse of the superball trick: we replace each ball of type i by i + 1 balls of the same type. (The new balls can be interpreted as external vertices as in [44].) This yields a new generalized Pólya urn with infinitely many types, all activities 1, and the transitions given by

ξ_{ij} = −(j + 1)δ_{ij} + δ_{0j} + (j + 1)δ_{i+1,j} = −(i + 1)δ_{ij} + δ_{0j} + (i + 2)δ_{i+1,j}.

Again we truncate and use the M +1 types 0, …, M only, with ξ_{Mj} changed to δ_{0j} + δ_{Mj}. (Note that such truncation does not work in the original urn model representing internal nodes.)

In this case the eigenvalues are λ_{1} = 2 and −1, −2, …, −M, see [29]. Theorem 3.17 applies for every M, and extends the joint asymptotic normality found by Mahmoud, Smythe and Szymański [44] to all degrees. Theorem 3.24 yields a functional limit theorem. The (co)variances are computed, using Lemmas 5.4 and 5.3, in [29].

Example 7.8 (rotations in a binary tree). Mahmoud [39] modelled rotations in the construction of a fringe-balanced binary tree by a generalized Pólya urn with three types, with X_{0} = (2, 0, 0)′ and ξ_{1} = (−2, 1, 2)′, ξ_{2} = ξ_{3} = (4, −1, −2)′. The number R_{n} of rotations in the n first insertions in the binary tree then equals the number of times a ball of type 3 is drawn, i.e. R_{n} = N_{n3}.

Note that ξ_{23}, ξ_{32} < 0, so (1.1) is violated.

Nevertheless, Mahmoud [39] observed that the proof of asymptotic normality in Smythe [51] holds for this urn too, and after finding exact formulas for the mean and variance he obtained

n^{−1/2}(R_{n} − \frac{2}{7}n) →^{d} N(0, \frac{66}{637}).  (7.4)

Although (1.1), and thus (A1), does not hold, we can derive this result from our theorems in several ways; since these methods may be useful in other applications where some ξ_{ij} are negative, we sketch three different approaches, leaving simple calculations to the reader. (It is instructive to compare the different calculations leading to the same result.)

First, note that in this urn, X_{n1} is even and X_{n3} = 2X_{n2}, which guarantees that we are never required to remove balls that do not exist. Moreover, λ_{1} = 1 with eigenvectors u_{1} = (1, 1, 1)′ and v_{1} = \frac{1}{7}(4, 1, 2)′. Thus, as asserted in Remark 4.2, our theorems hold for this urn too, and (7.4) follows by Theorem 3.21 and Lemma 5.6.

Secondly, since ξ_{2} = ξ_{3}, we may combine types 2 and 3 and consider the urn with two types and ξ_{1} = (−2, 3)′, ξ_{2} = (4, −3)′; to obtain R_{n} we add a dummy ball with probability 2/3 each time we draw a ball of type 2. This is a tenable urn (the same as in Example 7.5), and the result can be obtained by applying Corollary 3.11 and Lemma 5.3(i) to this urn (with dummy balls), stopping when X_{n1} + X_{n2} = n + 2.

Thirdly, consider the same 2-type urn again, but instead of adding dummy balls at random, observe that given the number N_{n2} of type 2 draws, R_{n} ∼ Bi(N_{n2}, 2/3). By Theorem 3.21, or by Example 7.5 and N_{n2} = \frac{1}{6}(3n − X_{n2}) = \frac{1}{6}(2n + w_{n} − w_{0}),

n^{−1/2}(N_{n2} − \frac{3}{7}n) →^{d} N(0, \frac{12}{637}).  (7.5)

In particular, N_{n2}/n →^{p} 3/7. Hence, the central limit theorem for the binomial distribution implies n^{−1/2}(R_{n} − \frac{2}{3}N_{n2}) →^{d} N(0, \frac{3}{7} · \frac{2}{3} · \frac{1}{3}); moreover, this holds jointly with (7.5), with independent limits. Thus

n^{−1/2}(R_{n}−\frac{2}{7}n) = n^{−1/2}(R_{n}−\frac{2}{3}N_{n2})+\frac{2}{3}n^{−1/2}(N_{n2}−\frac{3}{7}n) →^{d} N(0, \frac{2}{21}+(\frac{2}{3})^{2} \frac{12}{637}).

This also shows that of the variance in (7.4), only a fraction (\frac{2}{3})^{2} \frac{12}{66} = \frac{8}{99} comes from the random variation in the urn, i.e. from the shape of the tree.

In the second and third methods, we reduce to the tenable urn for external nodes in Example 7.5. We can replace this urn by the urn for internal nodes in Example 7.5, thus reducing the problem further to an urn that satisfies (A1). This is equivalent to reducing the original 3-type urn by an extension of the superball trick, where the superballs may combine balls of different types. In this case we have two types of superballs: the first represents two balls of type 1, while the second represents 3 balls, 1 of type 2 and 2 of type 3. (This yields two further ways of deriving (7.4).)

Example 7.9 (random m-ary search tree). In an m-ary search tree, where m ≥ 2 is a fixed integer, a node may contain up to m − 1 keys. The tree is constructed recursively, starting with an empty root node. Incoming keys are added to the root node until it is full; it then get m daughters, initially empty, and further keys are passed on to one of the daughters, where the procedure repeats. See e.g. [38, Chapter 3].

Let us say that a node containing i keys has type i. With random input, the number of nodes of different types is modelled by a generalized Pólya urn with m types 0, …, m − 1. A ball of type i with i < m − 2 has activity i + 1; if drawn, it is removed and replaced by a ball of type i + 1. A ball of type m − 2 has activity m − 1; if drawn, it is removed and replaced by a ball of type m − 1 and m balls of type 0. A ball of type m − 1 has activity 0. (Since balls of type m − 1 are dead, we can ignore them when studying the other types.)

Alternatively, we can study external vertices; each (internal) vertex of type i ≤ m − 2 has i + 1 external vertices, which we label with the same type i. The external vertices evolve as an urn with m − 1 types 0, …, m − 2, all activities 1, and the replacement rules ξ_{i} = −(i + 1)δ_{i} + (i + 2)δ_{i+1}, 0 ≤ i ≤ m − 3, and ξ_{m−2} = −(m − 1)δ_{m−2} + mδ_{0}. For example, for m = 4 the external version has the matrix

A = ( \begin{matrix}−1 & 0 & 4 \\ 2 & −2 & 0 \\ 0 & 3 & −3 \end{matrix} )

For the external version, it is easily seen that A has the characteristic polynomial ϕ_{m}(λ) := ∏_{i=1}^{m−1}(λ + i) − m!; the largest real root is λ_{1} = 1, see Lemma 5.4. (For the internal version, we have λϕ_{m}(λ) with an additional root λ = 0.)

A detailed study [42] shows that Re λ_{2} → 1 = λ_{1} as m → ∞, and that Re λ_{2} < \frac{1}{2} for m ≥ 26, but Re λ_{2} > \frac{1}{2} for m > 26. Hence, the numbers of nodes of different types have an asymptotic normal distribution when m ≤ 26, but, as rigorously shown by Chern and Hwang [16], not for larger m. This has earlier been shown by other methods [42], [36], [16], and by urns as here by [41].

Example 7.10 (a branched random walk). Let G be a finite group and ξ a random element of G, with some distribution μ. We define an urn process where the types are the elements of G and ξ_{gh} = δ_{gξ,h}. We thus draw a ball, replace it and add a new ball with a type shifted according to μ. This is a special type of a branching random walk on G.

Since exactly one ball is added each time, λ_{1} = 1. By the symmetry, u_{1} = (1, …, 1)′ and v_{1} = q^{−1}(1, …, 1)′. The matrix A operates by convolution on G: Av = v∗μ for v ∈ C^{q} = ℓ^{2}(G). Hence, if G is commutative, the characters of G are eigenvectors of A and the eigenvalues are the Fourier coefficients of μ: Aχ = μ̂(χ)χ, χ ∈ Ĝ. In particular, we see that n^{−1/2}(X_{n} − nv_{1}) converges to a Gaussian limit N(0, Σ) if and only if Re μ̂(χ) < 1/2 for all χ ≠ 1. Lemma 5.5 yields B = q^{−1}I (for any μ). Hence, Lemmas 5.4 and 5.3(i) together with the orthogonality of the characters yield

Σ = Σ_{I} = q^{−2} ∑_{χ≠1} \frac{1}{1 − 2 Re μ̂(χ)}χ\overline{χ}′.

For a concrete example, let G be the cyclic group ℤ_{q} and suppose that only nearest-neighbour shifts are allowed, i.e. μ is supported on ±1. (For example, we may always shift one step forward, or make a symmetric random choice each time.) Then Re μ̂(χ_{k}) = cos(2πk/q), k = 0, …, q − 1, and Re λ_{2} = cos(2π/q). Hence, if q ≤ 5, Re λ_{2} < 1/2 and X_{n} is asymptotically normal with variance of the order n; if q = 6, X_{n} is still asymptotically normal (but more degenerate) but the variance is of order n log n, and if q > 6, the variance is of larger order and X_{n} is not asymptotically normal.

For non-commutative G, we obtain similar results by considering the irreducible representations of G.

These results were proved in [25] by a different method (moment calculations). Moreover, [25] treats also infinite compact groups, obtaining the same results there. This suggests that the results in this paper may have generalizations to infinite sets of types, see Remark 4.1.

8. A lemma

We will later use a lemma on joint convergence in distribution. The lemma is a simple exercise in measure theory, but since we do not know any good reference, and the notation makes it look more complicated than it really is, we give a detailed statement and proof. We begin with a simpler version.

Lemma 8.1. Suppose that (η_{n}, ζ_{n}) are pairs of random variables with values in S_{1} × S_{2} for some separable metric spaces S_{1} and S_{2}, and that η and ζ are random variables with values in S_{1} and S_{2}, respectively. Suppose further:

(i) η_{n} →^{d} η as n → ∞;

(ii) for every measurable set A ⊆ S_{1} such that lim inf_{n→∞} ℙ(η_{n} ∈ A) > 0, it holds that, conditioned on η_{n} ∈ A, ζ_{n} →^{d} ζ.

Then we have joint convergence (η_{n}, ζ_{n}) →^{d} (η′, ζ′) as n → ∞, with η′ and ζ′ independent copies of η and ζ, respectively.

Proof. Suppose that A and B are measurable sets in S_{1} and S_{2} with ℙ(η ∈ ∂A) = ℙ(ζ ∈ ∂B) = 0. Then, by (i), ℙ(η_{n} ∈ A) → ℙ(η ∈ A), see e.g. [15, Theorem 2.1]. If ℙ(η ∈ A) > 0, we thus have by (ii) that, conditioned on η_{n} ∈ A, ζ_{n} →^{d} ζ, and thus ℙ(ζ_{n} ∈ B | η_{n} ∈ A) → ℙ(ζ ∈ B). Consequently, when ℙ(η ∈ A) > 0,

\begin{aligned}ℙ((η_{n}, ζ_{n}) ∈ A × B) & = ℙ(ζ_{n} ∈ B | η_{n} ∈ A) ℙ(η_{n} ∈ A) \\ & → ℙ(ζ ∈ B) ℙ(η ∈ A) = ℙ((η′, ζ′) ∈ A × B). \end{aligned} (8.1)

The same holds trivially if ℙ(η ∈ A) = 0, and (η_{n}, ζ_{n}) →^{d} (η′, ζ′) follows by [15, Theorem 3.1]. □

We will need the following extension. (The index n may be replaced by a continuous parameter t, since it suffices to consider sequences t_{n} → ∞.)

Lemma 8.2. Suppose that η_{n}, ζ_{n}, η, ζ, S_{1} and S_{2} are as in Lemma 8.1. Let further E_{1} ⊆ E_{2} ⊆ … be an increasing sequence of measurable subsets of S_{1}, let E = ⋃_{1}^{∞} E_{m}, and suppose that:

(i) η_{n} →^{d} η as n → ∞;

(ii) ℙ(η ∈ ∂E_{m}) = 0 for every m = 1, 2, …; (iii) ℙ(η_{n} ∈ E) → ℙ(η ∈ E) > 0 as n → ∞; (iv) for every m = 1, 2, … and every measurable set A ⊆ E_{m} such that

lim inf_{n→∞} ℙ(η_{n} ∈ A) > 0, it holds that, conditioned on η_{n} ∈ A,

ζ_{n} →^{d} ζ.

Then, conditioned on η_{n} ∈ E, we have ζ_{n} →^{d} ζ and joint convergence

L((η_{n}, ζ_{n}) | η_{n} ∈ E) →^{d} L((η′, ζ′) | η′ ∈ E),

with η′ and ζ′ independent copies of η and ζ, respectively.

Proof. Suppose again that A and B are measurable sets in S_{1} and S_{2} with ℙ(η ∈ ∂A) = ℙ(ζ ∈ ∂B) = 0. Let ε > 0 and choose m such that ℙ(η ∈ E \ E_{m}) < ε. By (i), (ii) and (iii) we have

\begin{aligned}ℙ(η_{n} ∈ E \ E_{m}) & = ℙ(η_{n} ∈ E) − ℙ(η_{n} ∈ E_{m}) \\ & → ℙ(η ∈ E) − ℙ(η ∈ E_{m}) = ℙ(η ∈ E \ E_{m}), \end{aligned}

and thus ℙ(η_{n} ∈ E \ E_{m}) < ε too for large n.

By (ii), ℙ(η ∈ ∂(A ∩ E_{m})) ≤ ℙ(η ∈ ∂A) + ℙ(η ∈ ∂E_{m}) = 0. Hence, (iv) implies as in (8.1)

ℙ((η_{n}, ζ_{n}) ∈ (A ∩ E_{m}) × B) → ℙ((η′, ζ′) ∈ (A ∩ E_{m}) × B)

and thus

|ℙ((η_{n}, ζ_{n}) ∈ (A ∩ E) × B) − ℙ((η′, ζ′) ∈ (A ∩ E) × B)|
= |ℙ((η_{n}, ζ_{n}) ∈ (A ∩ (E \ E_{m})) × B)
+ ℙ((η_{n}, ζ_{n}) ∈ (A ∩ E_{m}) × B) − ℙ((η′, ζ′) ∈ (A ∩ E_{m}) × B)
− ℙ((η′, ζ′) ∈ (A ∩ (E \ E_{m})) × B)|
≤ ℙ(η_{n} ∈ E \ E_{m})
+ |ℙ((η_{n}, ζ_{n}) ∈ (A ∩ E_{m}) × B) − ℙ((η′, ζ′) ∈ (A ∩ E_{m}) × B)|
+ ℙ(η ∈ E \ E_{m})
< 3ε,

provided n is large enough. Consequently, as n → ∞,

ℙ((η_{n}, ζ_{n}) ∈ (A ∩ E) × B) → ℙ((η′, ζ′) ∈ (A ∩ E) × B).

Dividing by ℙ(η_{n} ∈ E) → ℙ(η ∈ E) = ℙ(η′ ∈ E), we find

ℙ((η_{n}, ζ_{n}) ∈ A × B | η_{n} ∈ E) → ℙ((η′, ζ′) ∈ A × B | η′ ∈ E),

and the result follows by [15, Theorem 3.1]. □

Note that (iii) follows from (i) if ℙ(η ∈ ∂E) = 0. However, this stronger condition does not always hold in our applications.

9. Proof of Theorem 3.1

The proof of Theorem 3.1 is based on martingale theory, in particular a martingale convergence theorem by Jacod and Shiryaev [24]. The theorem uses the quadratic variation [X, X]_{t} of a martingale X defined on [0, ∞), and its bilinear extension [X, Y]_{t} to two martingales X and Y. For a general definition see e.g. [24] or [48]; for us it will suffice to know that, if X and Y are (real or complex) martingales of finite variation, then

[X, Y]_{t} = ∑_{0<s≤t} ΔX(s)ΔY (s),  (9.1)

where ΔX(s) := X(s) − X(s−) is the jump of X at s and, similarly, ΔY (s) := Y (s) − Y (s−). (A martingale X is said to be of finite variation if it is so pathwise, i.e. if t ↦ X(t) a.s. has bounded variation on each finite interval.) The sum in (9.1) is formally uncountable, but in reality countable since there is only a countable number of jumps; in the applications below, the sum will be finite. (There is some disagreement in the literature on the definition of [X, Y] in the case X(0)Y (0) ≠ 0; we have chosen the version with [X, Y]_{0} = 0.) For martingales of infinite variation (such as Brownian motion), (9.1) fails, but we have always the inequality

∑_{0<s≤t} |ΔX(s)|^{2} ≤ [X, X̄]_{t}.  (9.2)

For vector-valued martingales X = (X_{i})_{i=1}^{m} and Y = (Y_{j})_{j=1}^{n}, we define the square bracket [X, Y] to be the m × n matrix ([X_{i}, Y_{j}])_{i,j}.

For a real-valued martingale X, the quadratic variation [X, X]_{t} is a nonnegative and non-decreasing process. A real-valued martingale X(s) on [0, t] is an L^{2}-martingale if and only if 𝔼[X, X]_{t} < ∞ and 𝔼 |X(0)|^{2} < ∞, and then

𝔼 |X(t)|^{2} = 𝔼[X, X]_{t} + 𝔼 |X(0)|^{2}.  (9.3)

(For complex-valued martingales one has to consider [X, X̄]_{t}.) There is also a corresponding bilinear formula, which extends to (real or complex) vectorvalued L^{2}-martingales in the form

𝔼 X(t)Y′(t) = 𝔼[X, Y]_{t} + 𝔼 X(0)Y′(0).  (9.4)

We will use the following general result based on [24]; see [26] and [27] for similar versions. Again, n may be replaced by a continuous parameter t.

Proposition 9.1. (i) Assume that for each n, M_{n}(x) = (M_{ni}(x))_{i=1}^{q} is a real q-dimensional martingale on [0, ∞) with M_{n}(0) = 0, and that Σ(x), x ≥ 0, is a (non-random) continuous matrix-valued function such that for every fixed x ≥ 0,

[M_{n}, M_{n}]_{x} →^{p} Σ(x) as n → ∞,  (9.5)
sup_{n} 𝔼 |M_{n}(x)|^{2} < ∞.  (9.6)

Then M_{n} →^{d} M as n → ∞, in D[0, ∞), where M is a continuous qdimensional Gaussian process with 𝔼 M(x) = 0 and covariances

𝔼 M(x)M′(y) = Σ(x), 0 ≤ x ≤ y < ∞.  (9.7)

(ii) The same holds for complex M_{n} (with M complex), provided (9.5) is supplemented by

[M_{n}, \overline{M}_{n}]_{x} →^{p} Σ^{†}(x) as n → ∞,  (9.8)

for some continuous matrix-valued function Σ^{†}(x), and then further

𝔼 M(x)\overline{M}′(y) = Σ^{†}(x), 0 ≤ x ≤ y < ∞.

Proof. (i): Note first that (9.5) implies that if x ≤ y, then Σ(y) − Σ(x) is positive semidefinite, so there exists a q-dimensional Gaussian process M with independent increments such that 𝔼 M(x) = 0 and 𝔼 M(x)M′(y) = Σ(x), x ≤ y, see e.g. [24, Theorem II.5.2]. M is a martingale, and it is continuous because each component is a (deterministic) time change of a Brownian motion.

Since, by (9.2) and (9.3),

𝔼 sup_{y≤x} |ΔM_{ni}(y)|^{2} ≤ 𝔼 ∑_{y≤x} |ΔM_{ni}(y)|^{2} ≤ 𝔼[M_{ni}, M_{ni}]_{x} = 𝔼 |M_{ni}(x)|^{2},

it follows from (9.6) that, for each fixed x > 0, the sequence sup_{y≤x} |ΔM_{ni}(y)| is uniformly integrable for each i, and thus sup_{y≤x} |ΔM_{n}(s)| is uniformly integrable. The result now follows from [24, Theorem VIII.3.12, (ii)⇒(i)].

(ii): This follows from (i) by considering the real 2q-dimensional martingales (Re M_{n}, Im M_{n}). Note that (9.5) for these (with the appropriate right-hand sides) follows from (9.5) and (9.8). □

In order to apply this result to our process X (t), we have to first define a suitable martingale, and then estimate its quadratic variation. The martingale is a standard one in branching process theory, and the estimates will be derived by standard methods too, although the details will take some time. We proceed with a series of lemmas. (Some of these are known, but included here for completeness and because our conditions are slightly more general than the standard ones; see also Remark 4.2.) We make the definition

Y(t) := e^{−tA}X (t)  (9.9)

and begin with a fundamental well-known result, cf. [9, Theorem V.8.1].

Lemma 9.2. Y(t) is a martingale for t ≥ 0. In particular,

𝔼 X (t) = e^{tA} 𝔼 X (0)  (9.10)

and thus

𝔼 X (t) = O(e^{λ_{1}t}).  (9.11)

Proof. It follows from the definitions of X and A that

\frac{d}{dt} 𝔼 X_{j}(t) = ∑_{i=1}^{q} a_{i} 𝔼 X_{i}(t) 𝔼 ξ_{ij} = (A 𝔼 X (t))_{j}  (9.12)

and hence \frac{d}{dt} 𝔼 X (t) = A 𝔼 X (t). This yields (9.10) by integration, and the martingale property follows from (9.10) and the Markov property. Finally, (2.12) implies (9.11). □

Let 0 < τ_{i1} < τ_{i2} < … denote the times a ball of type i splits, and let N_{i}(t) := #{k : τ_{ik} ≤ t} be the number of such splits up to time t. Since the martingale Y(t) := e^{−tA}X (t) has finite variation, its quadratic variation is by (9.1) given by its jumps

\begin{aligned}[Y, Y]_{t} & = ∑_{i=1}^{q} ∑_{k:τ_{ik}≤t} ΔY(τ_{ik})ΔY′(τ_{ik}) \\ & = ∑_{i=1}^{q} ∑_{k:τ_{ik}≤t} e^{−τ_{ik}A}ΔX (τ_{ik})ΔX′(τ_{ik})e^{−τ_{ik}A′}. \end{aligned}  (9.13)

The main part of the proof consists of estimating this sum, and components of it. It will be convenient to state a general lemma for sums of this type, for simplicity considering a single i.

Lemma 9.3. Fix i ∈ {1, …, q}. Let f_{1} and f_{2} be continuous matrix-valued functions defined on [0, ∞) and let g be a matrix-valued function on ℝ^{q} such that 𝔼 |g(ξ_{i})| < ∞. Suppose further that the dimensions of f_{1}, g, and f_{2} are such that the product f_{1}gf_{2} is defined. Let

Z(t) := ∑_{k:τ_{ik}≤t} f_{1}(τ_{ik})g(ΔX (τ_{ik}))f_{2}(τ_{ik}) = ∫_{0}^{t} f_{1}(s)g(ΔX (s))f_{2}(s) dN_{i}(s)

(9.14)

and

Z̃(t) := Z(t) − ∫_{0}^{t} f_{1}(s)(𝔼 g(ξ_{i}))f_{2}(s)a_{i}X_{i}(s) ds.  (9.15)

Then Z̃(t), t ≥ 0, is a (matrix-valued) martingale; in particular

𝔼 Z(t) = ∫_{0}^{t} f_{1}(s) 𝔼 g(ξ_{i})f_{2}(s)a_{i} 𝔼 X_{i}(s) ds.  (9.16)

Proof. In the special case when f_{1} = f_{2} = g = 1, Z(t) = N_{i}(t) and Z̃(t) = Ñ_{i}(t) := N_{i}(t) − ∫_{0}^{t} a_{i}X_{i}(s) ds. The fact that Ñ_{i}(t) is a martingale is a wellknown simple consequence of the assumption that the balls have independent exponential lifetimes. The present extension can be proved in the same way, because ΔX (τ_{ik}) is independent of the previous history; for example, a straightforward calculation shows that Z̃(t ∧ τ_{n}) − Z̃(t ∧ τ_{n−1}) (with τ_{0} := 0) is a martingale for each n, and the result follows by summing over n. We omit the details. □

Note that (9.16) with |f_{1}|, |f_{2}| and |g|, implies, using (9.11), that 𝔼 |Z(t)| < ∞ and that Z̃(t) is a uniformly integrable martingale on any finite interval [0, T].

From (9.4), (9.13), (2.13) and Lemma 9.3 follows

\begin{aligned}𝔼 Y(t)Y′(t) & = 𝔼[Y, Y]_{t} + 𝔼 Y(0)Y′(0) \\ & = ∑_{i=1}^{q} ∫_{0}^{t} e^{−sA}B_{i}e^{−sA′}a_{i} 𝔼 X_{i}(s) ds + 𝔼 Y(0)Y′(0). \end{aligned} (9.17)

By (9.9), this yields the following formula from [9, §V.7.3].

\begin{aligned}𝔼 X (t)X′(t) & = ∑_{i=1}^{q} ∫_{0}^{t} e^{(t−s)A}B_{i}e^{(t−s)A′}a_{i} 𝔼 X_{i}(s) ds + e^{tA} 𝔼 X (0)X′(0)e^{tA′}. \\ & (9.18) \end{aligned}

We apply the generalized eigenspace decomposition given by {P_{λ}} to (9.17) and obtain the following estimate. (See also Lemma 10.6 below.)

Lemma 9.4. (i) If Re λ ≤ λ_{1}/2, then

𝔼 |P_{λ}Y(t)|^{2} = O((1 + t)^{2d_{λ}+1}e^{(λ_{1}−2 Re λ)t}).

(ii) If Re λ > λ_{1}/2, then

𝔼 |P_{λ}Y(t)|^{2} = O(1).

Proof. By (9.17), (9.11) and (2.10),

\begin{aligned}|𝔼(P_{λ}Y(t)Y′(t)P′_{λ})| & ≤ C_{1} ∫_{0}^{t}|P_{λ}e^{−sA}| e^{λ_{1}s} |e^{−sA′}P′_{λ}| ds + C_{2} \\ & ≤ C_{3} ∫_{0}^{t}(1 + s)^{2d_{λ}}e^{(λ_{1}−2 Re λ)s} ds + C_{2}. \end{aligned}

The result follows by integration. □

An immediate consequence of Lemmas 9.2 and 9.4(ii) and the martingale convergence theorem for L^{2}-bounded martingales is the following [9, Theorem V.8.2]. We let F_{t} denote the σ-field generated by X (s), 0 ≤ s ≤ t.

Lemma 9.5. If Re λ > λ_{1}/2, then there exists a random vector W̃_{λ} ∈ E_{λ} such that P_{λ}Y(t) → W̃_{λ} as t → ∞, a.s. and in L^{2}. Moreover, P_{λ}Y(t) = 𝔼(W̃_{λ} | F_{t}). □

We can now prove part (iii) of Theorem 3.1. By (9.9) and (2.8),

e^{−tλ}Q_{λ,k}P_{λ}X (t) = e^{−tλ}Q_{λ,k}P_{λ}e^{tA}Y(t) = ∑_{j=0}^{k} \frac{t^{j}}{j!}Q_{λ,k}N_{λ}^{j}P_{λ}Y(t).

The first part of (iii) of Theorem 3.1 now follows from Lemma 9.5, with W_{λ,k} := (k!)^{−1}Q_{λ,k}N_{λ}^{k}W̃_{λ}. The second part (convergence in distribution to a constant process) follows because the first part implies uniform convergence a.s. for a ≤ x ≤ b, on every compact interval [a, b] ⊂ (0, ∞).

Before proving parts (i) and (ii), we observe that the limit in (iii) is nontrivial.

Lemma 9.6. If λ ∈ Λ_{III} and 0 ≤ k ≤ d_{λ}, then W_{λ,k} is non-degenerate. More precisely, ℙ(W_{λ,k} = 0) < 1 and ℙ(W_{λ,k} = w) = 0 for every w ≠ 0.

Proof. If W_{λ,k} = 0 a.s., then N_{λ}^{k}W̃_{λ} ∈ N_{λ}^{k+1}E_{λ} and thus N_{λ}^{d_{λ}} W̃_{λ} = 0 a.s.; hence, by Lemma 9.5, for any t,

N_{λ}^{d_{λ}}P_{λ}Y(t) = 𝔼(N_{λ}^{d_{λ}} W̃_{λ} | F_{t}) = 0 a.s.

and thus

N_{λ}^{d_{λ}}P_{λ}X (t) = e^{tA}N_{λ}^{d_{λ}}P_{λ}Y(t) = 0 a.s.  (9.19)

Considering first rational t and using the right-continuity, we see that (9.19) a.s. holds for all t ≥ 0. It follows that a.s. N_{λ}^{d_{λ}}P_{λ}ΔX (t) = 0 for all t > 0, and thus N_{λ}^{d_{λ}}P_{λ}ξ_{i} = 0 a.s. for every i with a_{i} > 0, since every such transition occurs with positive probability. Taking the expectation we find, since a_{i} 𝔼 ξ_{i} equals the i:th column of A by definition, N_{λ}^{d_{λ}}P_{λ}A = 0. By (2.5) this yields

0 = N_{λ}^{d_{λ}}P_{λ}A = λN_{λ}^{d_{λ}}P_{λ}.

Since N_{λ}^{d_{λ}}P_{λ} ≠ 0 (it is P_{λ} when d_{λ} = 0 and N_{λ}^{d_{λ}} otherwise), this implies λ = 0, contradicting λ ∈ Λ_{III}. Consequently, ℙ(W_{λ,k} = 0) < 1.

Next, let τ be the time of the first death, and let X̂(t) := X (τ + t). Then X̂ is a branching process with the same transitions as X but a different (random) initial state X̂(0) = X (τ); moreover, X̂ is independent of τ. Letting Ŵ_{λ,k} denote the limit corresponding to W_{λ,k} for X̂, one easily finds W_{λ,k} = e^{−λτ} Ŵ_{λ,k} a.s. Since τ and Ŵ_{λ,k} are independent and τ has a continuous distribution, ℙ(W_{λ,k} = w) = 0 follows for w ≠ 0 by conditioning on Ŵ_{λ,k}. □

Because λ_{1} > 0 by our assumption (A3), Lemma 9.5 applies in particular to λ = λ_{1}. Note that W_{λ_{1}} = W_{λ_{1},0} = W̃_{λ_{1}}. We write W := u_{1}·W_{λ_{1}} = u′_{1}W_{λ_{1}} and have by (2.7)

W_{λ_{1}} = P_{λ_{1}}W_{λ_{1}} = v_{1}u′_{1}W_{λ_{1}} = v_{1}W = W v_{1},  (9.20)

as asserted in Theorem 3.1.

Some well-known properties of W are collected in the next lemma, see [9, Theorems V.6.2 and V.7.2].

Lemma 9.7. We have W ≥ 0 a.s., ℙ(W > 0) > 0 and ℙ(W = w) = 0 for every w ≠ 0. Moreover,

(i) As t → ∞, u_{1} · Y(t) = e^{−λ_{1}t}u_{1} · X (t) →^{a.s.} W.

(ii) If u_{1} · Y(t) = 0 for some t then a.s. W = 0.

(iii) Conversely, ℙ(u_{1} · Y(t) > 0 and W = 0) → 0 as t → ∞. Consequently, a.s. W = 0 if and only if X becomes essentially extinct.

Proof. Since P′_{λ_{1}}u_{1} = u_{1}, we have by Lemma 9.5 and (9.20), (2.3)

u_{1} · Y(t) = P′_{λ_{1}}u_{1} · Y(t) = u_{1} · P_{λ_{1}}Y(t) →^{a.s.} u_{1} · W_{λ_{1}} = u_{1} · v_{1}W = W.

Moreover, A′u_{1} = λ_{1}u_{1} and thus

u_{1} · Y(t) = u_{1} · e^{−tA}X (t) = e^{−tA′}u_{1} · X (t) = e^{−λ_{1}t}u_{1} · X (t),

proving (i). Since u_{1} · X (t) ≥ 0, W ≥ 0 follows. If W = 0 a.s., then W_{λ_{1}} = 0 a.s. by (9.20), which contradicts Lemma 9.6. Similarly, if w ≠ 0 and ℙ(W = w) > 0, then ℙ(W_{λ_{1}} = wv_{1}) > 0, which again contradicts Lemma 9.6.

For the second part, recall that u_{1i} > 0 for i ∈ C_{1} and u_{1i} = 0 for i ∉ C_{1}. Hence, u_{1} · Y(t) = e^{−λ_{1}t}u_{1} · X (t) = 0 if and only if X_{i}(t) = 0 for all i ∈ C_{1}. (In other words, X is essentially extinct by time t.) Letting X̃ denote the branching process obtained from X by ignoring all balls of types not in C_{1}, this means that X̃(t) = 0. As a consequence, for all x ≥ t, X̃(x) = 0 and thus u_{1} · Y(x) = 0; hence the limit W = 0.

For the converse, we again consider the process X̃(t). This is an irreducible continuous time branching process, and by [9, Theorem V.7.2], W = 0 a.s. implies extinction, X̃(t) = 0 for large t, and thus u_{1} · Y(t) = 0 for large t. The result follows. □

We may further improve this result to the first claim in Theorem 3.1.

Lemma 9.8. As t → ∞, e^{−λ_{1}t}X (t) →^{a.s.} W v_{1}.

Proof. This is Theorem V.7.2 in [9], see also [4], but since our setting is somewhat more general, we give a complete proof.

Fix an eigenvalue λ ≠ λ_{1} and let δ := λ_{1} − Re λ > 0. Let ε > 0 and let E_{n} be the event sup_{t∈[n−1,n]}|e^{−λ_{1}t}P_{λ}X (t)| > ε.

If E_{n} occurs, then, by (2.10), for some t ∈ [n − 1, n],

εe^{λ_{1}t} < |P_{λ}X (t)| = |P_{λ}e^{tA}P_{λ}Y(t)| ≤ Cn^{d_{λ}}e^{t(λ_{1}−δ)}|P_{λ}Y(t)|

and thus |P_{λ}Y(t)| ≥ cεn^{−d_{λ}}e^{tδ}. Consequently, using Doob’s inequality and Lemma 9.4,

ℙ(E_{n}) ≤ ℙ(sup_{t≤n} |P_{λ}Y(t)| ≥ cεn^{−d_{λ}}e^{(n−1)δ})
≤ Cε^{−2}n^{2d_{λ}}e^{−2nδ} 𝔼 |P_{λ}Y(n)|^{2}
≤ Cε^{−2}n^{4d_{λ}+1}e^{−2nδ}e^{max(0,2δ−λ_{1})n}.

Hence ∑_{n=1}^{∞} ℙ(E_{n}) < ∞, so by the Borel–Cantelli lemma a.s. lim sup_{t→∞} |e^{−λ_{1}t}P_{λ}X (t)| ≤ ε. In other words,

e^{−λ_{1}t}P_{λ}X (t) → 0 a.s.,  λ ≠ λ_{1}.  (9.21)

Moreover, P_{λ_{1}}e^{−tA} = e^{−λ_{1}t}P_{λ_{1}} and thus, by Lemma 9.5 and (9.20),

e^{−λ_{1}t}P_{λ_{1}}X (t) = P_{λ_{1}}Y(t) →^{a.s.} W̃_{λ_{1}} = W_{λ_{1}} = W v_{1}.  (9.22)

The result follows by (9.21), (9.22) and (2.4). □

Lemma 9.9. Let, for t, y ≥ 0, Z_{y}(t) be defined as in Lemma 9.3 using matrix functions f_{1y}, f_{2y} and g. Suppose further that, as t → ∞ and for any fixed T,

sup_{t>1} ∫_{0}^{t} |f_{1t}(s)| |f_{2t}(s)|e^{λ_{1}s} ds < ∞,  (9.23)
∫_{0}^{t} |f_{1t}(s)|^{2}|f_{2t}(s)|^{2}e^{λ_{1}s} ds → 0,  (9.24)
∫_{0}^{T} |f_{1t}(s)| |f_{2t}(s)| ds → 0.  (9.25)

Then, as t → ∞,

Z_{t}(t) − a_{i}v_{1i}W ∫_{0}^{t} f_{1t}(s) 𝔼 g(ξ_{i})f_{2t}(s)e^{λ_{1}s} ds →^{p} 0.

Proof. We begin by showing that, defining Z̃_{t} as in (9.15),

𝔼 |Z̃_{t}(t)| → 0.  (9.26)

By considering the components separately, and taking real and imaginary parts, we may assume that f_{1y}, f_{2y}, g and Z_{y} are real-valued; we write f_{y} = f_{1y}f_{2y}.

Assume first that 𝔼 |g(ξ_{i})|^{2} < ∞. Since the martingale Z̃_{y}(t), t ≥ 0, has finite variation, its quadratic variation is by (9.1) given by

[Z̃_{y}, Z̃_{y}]_{t} = ∑_{s≤t}(ΔZ_{y}(t))^{2} = ∑_{k:τ_{ik}≤t} f_{y}(τ_{ik})^{2}g(ΔX (τ_{ik}))^{2}.

This is a sum of the same type as Z, with f_{1}, g and f_{2} replaced by f_{y}^{2}, g^{2} and 1. Hence Lemma 9.3 yields, together with (9.3) and (9.11),

𝔼 Z̃_{y}(t)^{2} = 𝔼[Z̃_{y}, Z̃_{y}]_{t} = ∫_{0}^{t} f_{y}(s)^{2} 𝔼 g(ξ_{i})^{2}a_{i} 𝔼 X_{i}(s) ds

≤ C ∫_{0}^{t} f_{y}(s)^{2}e^{λ_{1}s} ds = C ∫_{0}^{t} f_{1y}(s)^{2}f_{2y}(s)^{2}e^{λ_{1}s} ds

and thus, by (9.24), 𝔼 |Z̃_{t}(t)|^{2} → 0, which proves (9.26).

In the case 𝔼 |g(ξ_{i})|^{2} = ∞ we truncate, defining g_{1}(x) = g(x)1[|g(x)| ≤ M] and g_{2}(x) = g(x)−g_{1}(x) for a constant M. Thus g = g_{1} +g_{2} and there is a corresponding decomposition Z̃_{y}(t) = Z̃_{1y}(t)+Z̃_{2y}(t), where 𝔼 |Z̃_{1t}(t)| → 0 by the case just proved. Applying Lemma 9.3 with |f_{1y}|, |f_{2y}| and |g_{2}|,

\begin{aligned}𝔼 |Z̃_{2y}(t)| & ≤ 𝔼 |Z_{2y}(t)| + ∫_{0}^{t} |f_{1y}(s)| 𝔼 |g_{2}(ξ_{i})||f_{2y}(s)|a_{i} 𝔼 X_{i}(s) ds \\ & ≤ 2 ∫_{0}^{t} |f_{1y}(s)| 𝔼 |g_{2}(ξ_{i})||f_{2y}(s)|a_{i} 𝔼 X_{i}(s) ds. \end{aligned}

Hence, by (9.11) and (9.23), for t ≥ 1,

𝔼 |Z̃_{2t}(t)| ≤ C_{1} 𝔼 |g_{2}(ξ_{i})| ∫_{0}^{t} |f_{1t}(s)| |f_{2t}(s)|e^{λ_{1}s} ds ≤ C_{2} 𝔼 |g_{2}(ξ_{i})|

and thus

lim sup_{t→∞} 𝔼 |Z̃_{t}(t)| = lim sup_{t→∞} 𝔼 |Z̃_{2t}(t)| ≤ C_{2} 𝔼 |g_{2}(ξ_{i})|.

Since 𝔼 |g_{2}(ξ_{i})| can be made arbitrarily small by choosing M large, (9.26) follows.

Next, for every ε > 0, there is by Lemma 9.8 a.s. a random T such that |e^{−λ_{1}t}X_{i}(t) − v_{1i}W | < ε for t ≥ T. Hence, for t ≥ T, for some random K independent of t and again using (9.23),

|∫_{0}^{t} f_{1t}(s) 𝔼 g(ξ_{i})f_{2t}(s)a_{i}(X_{i}(s) − v_{1i}W e^{λ_{1}s})ds|

≤ C_{3} ∫_{0}^{T} |f_{1t}(s)||f_{2t}(s)|(X_{i}(s) + v_{1i}W e^{λ_{1}s})ds + C_{4} ∫_{T}^{t} |f_{1t}(s)||f_{2t}(s)|εe^{λ_{1}s} ds

≤ K ∫_{0}^{T} |f_{1t}(s)||f_{2t}(s)| ds + C_{5}ε.

It follows by (9.25) and the arbitrariness of ε that

∫_{0}^{t} f_{1t}(s) 𝔼 g(ξ_{i})f_{2t}(s)a_{i}(X_{i}(s) − v_{1i}W e^{λ_{1}s})ds → 0 a.s.,  (9.27)

and thus in probability, as t → ∞. The result follows from (9.15), (9.26) and (9.27). □

Lemma 9.9 and (9.13) will give us the estimates of [Y, Y]_{t} that we need in order to apply Proposition 9.1. We begin with P_{I}Y, the part corresponding to eigenvalues λ with Re λ < λ_{1}/2. For convenience we set Ã:= P_{I}A−λ_{1}/2, and note that each eigenvalue μ of Ã satisfies Re μ < 0. Hence, for some δ > 0 and every t ≥ 0,

|e^{tÃ}| = O(e^{−δt}).  (9.28)

Define, recalling (9.9),

Y_{y}(t) := e^{yÃ}P_{I}Y(t) = e^{−λ_{1}y/2+(y−t)A}P_{I}X (t).  (9.29)

For fixed y, Y_{y}(t) is a martingale by Lemma 9.2, and, by (9.13),

[Y_{y}, Y_{y}]_{t} = e^{yÃ}P_{I}[Y, Y]_{t}P′_{I}e^{yÃ′} = ∑_{i=1}^{q} ∑_{k:τ_{ik}≤t} f_{1y}^{(i)}(τ_{ik})g(ΔX (τ_{ik}))f_{2y}^{(i)}(τ_{ik}),

(9.30)

where g(ξ) = ξξ′ and

f_{1y}^{(i)}(s) = f_{2y}^{(i)}′(s) = e^{yÃ}P_{I}e^{−sA} = e^{−λ_{1}s/2+(y−s)Ã}P_{I}.  (9.31)

The inner sum in (9.30) is of the type studied in Lemmas 9.3 and 9.9. We apply Lemma 9.9 for each i separately. By (9.31) and (9.28),

|f_{1t}^{(i)}(s)| = |f_{2t}^{(i)}(s)| = O(e^{−λ_{1}s/2−(t−s)δ}),  0 ≤ s ≤ t,  (9.32)

and the conditions (9.23)–(9.25) follow.

Moreover, 𝔼 g(ξ_{i}) = B_{i} by (2.13), and by (9.31)

∫_{0}^{t} f_{1t}^{(i)}(s) 𝔼 g(ξ_{i})f_{2t}^{(i)}(s)e^{λ_{1}s} ds = ∫_{0}^{t} e^{(t−s)Ã}P_{I}B_{i}P′_{I}e^{(t−s)Ã′} ds

= P_{I} ∫_{0}^{t} e^{rÃ}B_{i}e^{rÃ′} dr P′_{I} → P_{I} ∫_{0}^{∞} e^{rÃ}B_{i}e^{rÃ′} dr P′_{I},

where the integral converges by (9.28). Consequently, (9.30) and Lemma 9.9 yield, as t → ∞, see (2.15),

[Y_{t}, Y_{t}]_{t} →^{p} ∑_{i=1}^{q} a_{i}v_{1i}W P_{I} ∫_{0}^{∞} e^{sÃ}B_{i}e^{sÃ′} ds P′_{I} = W Σ_{I}.  (9.33)

For any fixed real x, Y_{t}(s) = e^{−xÃ}Y_{t+x}(s) and thus, as t → ∞,

[Y_{t}, Y_{t}]_{t+x} = e^{−xÃ}[Y_{t+x}, Y_{t+x}]_{t+x}e^{−xÃ′} →^{p} W e^{−xÃ}Σ_{I}e^{−xÃ′}.  (9.34)

Moreover, by (9.30) and Lemma 9.3

𝔼[Y_{t}, Y_{t}]_{t} = ∑_{i=1}^{q} ∫_{0}^{t} f_{1t}^{(i)}(s)B_{i}f_{2t}^{(i)}(s)a_{i} 𝔼 X_{i}(s) ds

which by (9.32) and (9.11) is bounded for t ≥ 0. Consequently, for t+x ≥ 0,

𝔼[Y_{t}, Y_{t}]_{t+x} = e^{−xÃ} 𝔼[Y_{t+x}, Y_{t+x}]_{t+x}e^{−xÃ′} = O(|e^{−xÃ}|^{2}).

We have further

|Y_{t}(0)| = |e^{tÃ}P_{I}Y(0)| = O(|e^{tÃ}|).

and thus, using (9.3) on each component of Y_{t},

𝔼 |Y_{t}(t + x)|^{2} = |Y_{t}(0)|^{2} + Tr 𝔼[Y_{t}, Y_{t}]_{t+x} = O(|e^{tÃ}|^{2} + |e^{−xÃ}|^{2}). (9.35)

In particular, for any fixed x, by (9.28),

sup_{t≥−x} 𝔼 |Y_{t}(t + x)|^{2} < ∞.  (9.36)

We are close to applying Proposition 9.1, but we still have the obstacle that the right-hand side of (9.34) is random. To simply divide Y_{t} by \sqrt{W} would destroy the martingale property, but we can do something similar: Let, for t ≥ 1, h = h(t) := t^{1/2} (any function of t that increases slowly to ∞ would do), and define (with the convention 0^{−1/2} = 0)

Ỹ_{t}(x) := \begin{cases}(u_{1} · Y(t − h))^{−1/2}(Y_{t}(t + ln x) − Y_{t}(t − h)), & x ≥ e^{−h}, \\ 0, & x < e^{−h}. \end{cases} (9.37)

Clearly, Ỹ_{t} is a martingale on [0, ∞), and this is still true conditioned on any event E_{t} ∈ F_{t−h}.

For fixed x > 0 we have from (9.1), for t so large that −h(t) < ln x,

[Ỹ_{t}, Ỹ_{t}]_{x} = (u_{1} · Y(t − h))^{−1}([Y_{t}, Y_{t}]_{t+ln x} − [Y_{t}, Y_{t}]_{t−h}).  (9.38)

Further, using Doob’s inequality, (9.35) and (9.28),

𝔼 sup_{x≤t−h} |Y_{t}(x)|^{2} ≤ 4 𝔼 |Y_{t}(t − h)|^{2} = O(e^{−2δt} + e^{−2δh}) = o(1), (9.39)

as t → ∞. It follows by (9.3) applied to the components Y_{ti} that, for i = 1, …, q, 𝔼[Y_{ti}, Y_{ti}]_{t−h} → 0. The same holds for the nondiagonal entries of 𝔼[Y_{t}, Y_{t}]_{t−h} by the Cauchy–Schwarz inequality and the Kunita–Watanabe inequality |[Y_{ti}, Y_{tj}]_{s}| ≤ [Y_{ti}, Y_{ti}]_{s}^{1/2}[Y_{tj}, Y_{tj}]_{s}^{1/2}, which in this setting is the Cauchy–Schwarz inequality applied to (9.1). Hence,

𝔼[Y_{t}, Y_{t}]_{t−h} → 0.  (9.40)

Since t − h → ∞, Lemma 9.7 yields u_{1} · Y(t − h) →^{a.s.} W. Combined with (9.38), (9.34) and (9.40), this implies that, on the event {W > 0}, for every fixed x > 0,

[Ỹ_{t}, Ỹ_{t}]_{x} →^{p} e^{−(ln x)Ã}Σ_{I}e^{−(ln x)Ã′}.  (9.41)

Now assume that ε > 0 and that, for large t, E_{t} is any event in F_{t−h} with ℙ(E_{t}) ≥ ε and u_{1} · Y(t − h) ≥ ε on E_{t}. Since ℙ(E_{t} ∩ {W = 0}) → 0 by Lemma 9.7, it follows from (9.41), (9.36) and (9.39) that Ỹ_{t} conditioned on E_{t} satisfies the conditions of Proposition 9.1 with Σ(x) := e^{−(ln x)Ã}Σ_{I}e^{−(ln x)Ã′}. (Σ(0) = 0; note that Σ is continuous at 0 by (9.28).) Consequently, conditioned on E_{t}, Ỹ_{t}(x) →^{d} Ũ_{I}(x) in D[0, ∞) as t → ∞, where Ũ_{I} is the continuous vector-valued Gaussian stochastic process with mean 0 and covariances

𝔼 Ũ_{I}(x)Ũ_{I}(y)′ = 𝔼 Ũ_{I}(x)Ũ_{I}(x)′ = Σ(x) = e^{−(ln x)Ã}Σ_{I}e^{−(ln x)Ã′}, x ≤ y.

(9.42)

Next, from the definition of Ỹ_{t} and (9.39), on E_{t},

sup_{x≥0}|Ỹ_{t}(x)−(u_{1}·Y(t−h))^{−1/2}Y_{t}(t+ln x)| = (u_{1}·Y(t−h))^{−1/2} sup_{y≤t−h} |Y_{t}(y)| →^{p} 0

so we may replace Ỹ_{t}(x) by (u_{1} · Y(t − h))^{−1/2}Y_{t}(t + ln x) in this limit result. Changing variables x → e^{x} we find, in D[−∞, ∞) and thus in D(−∞, ∞), conditioned on E_{t}, as t → ∞,

(u_{1} · Y(t − h))^{−1/2}Y_{t}(t + x) →^{d} Ũ_{I}(e^{x}).

Since, by (9.29), Y_{t}(t + x) = e^{−λ_{1}t/2−xA}P_{I}X (t + x), we obtain by multiplication by e^{xÃ}, as t → ∞ and conditioned on E_{t}, in D(−∞, ∞),

(u_{1} · Y(t − h))^{−1/2}e^{−λ_{1}(t+x)/2}P_{I}X (t + x) →^{d} U_{I}(x) := e^{xÃ}Ũ_{I}(e^{x}).

U_{I} is a Gaussian process, and the covariances (3.1) follow from (9.42). We now apply Lemma 8.2 taking η_{t} = u_{1} · Y(t − h), ζ_{t} = (u_{1} · Y(t −

h))^{−1/2}e^{−λ_{1}(t+x)/2}P_{I}X (t + x), η = W, ζ = U_{I}, S_{1} = ℝ, S_{2} = D(−∞, ∞), E = {r ∈ ℝ : r > 0}, and E_{m} = {r ∈ ℝ : r > 1/m}.

We have just proved assumption (iv) in Lemma 8.2, and assumptions (i), (ii) and (iii) follow by Lemma 9.7. Hence Lemma 8.2 shows that, conditioned on u_{1} · Y(t − h) > 0, (η_{t}, ζ_{t}) →^{d} (W̃, U_{I}), where W̃ has the distribution of W conditioned on being positive and with U_{I} independent of W̃. Consequently, conditioned on u_{1} · Y(t − h) > 0,

e^{−λ_{1}(t+x)/2}P_{I}X (t + x) = η_{t}^{1/2}ζ_{t} →^{d} W̃^{1/2}U_{I}(x) in D(−∞, ∞). (9.43)

If ℙ(W = 0) = 0, then W̃ = W and u_{1} · Y(t − h) > 0 a.s. by Lemma 9.7, and part (i) of Theorem 3.1 is proved.

If ℙ(W = 0) > 0, we define the martingale

Ŷ_{t}(x) := \begin{cases}Y_{t}(t + ln x) − Y_{t}(t − h), & x ≥ e^{−h}, \\ 0, & x < e^{−h}. \end{cases}

Since u_{1} · Y(t − h) = 0 implies W = 0 a.s. by Lemma 9.7, (9.34) implies that, conditioned on u_{1} · Y(t − h) = 0,

[Ŷ_{t}, Ŷ_{t}]_{x} →^{p} 0

for every fixed x. Together with (9.36) and (9.39), this shows that we can apply Proposition 9.1 to Ŷ_{t}(x) conditioned on u_{1} · Y(t − h) = 0, now with Σ(x) = 0 and thus M = 0. Hence, conditioned on u_{1} · Y(t − h) = 0, we find first Ŷ_{t}(x) →^{d} 0 in D[0, ∞), then by (9.39) and a change of variable Y_{t}(t + x) →^{d} 0 in D(−∞, ∞), and then, after multiplication by e^{xÃ′},

e^{−λ_{1}(t+x)/2}P_{I}X (t + x) →^{d} 0 = W^{1/2}U_{I}(x) in D(−∞, ∞).

This complements (9.43) and together they imply part (i) of Theorem 3.1, with U_{I} independent of W.

A simple modification of the argument shows that (i) and (iii) hold jointly. (When we talk about (iii) holding jointly with other assertions, here and later in the proof, we only mean the part on convergence in distribution.) Indeed, let Y_{III,t}(xt) denote the collection of the left-hand sides of (iii) in Theorem 3.1 (for λ ∈ Λ_{III} and 0 ≤ k ≤ d_{λ}), regarded as a large vector-valued process, and let W_{III} denote the corresponding collection of the right-hand sides. The a.s. convergence in (iii) implies that Y_{III,t}(xt) − Y_{III,t}(h) →^{p} 0 in D(0, ∞) as functions of x. Hence the convergence in distribution in (iii) is equivalent to Y_{III,t}(h) →^{d} W_{III}.

By enlarging η_{t} (and S_{1}) in the application of Lemma 8.2 above, changing η_{t} into (u_{1} · Y(t − h), Y_{III,t}(h)), it now follows that (i) and (iii) hold jointly, with U_{I} independent of {W_{λ,k}}_{λ∈Λ_{III},k≤d_{λ}}.

For λ ∈ Λ_{II} we argue similarly, now putting

Y_{y,λ,k}(t) = y^{−k−1/2}Q_{λ,k}P_{λ}Y(t).  (9.44)

Again, this is a martingale (for fixed y, λ, k). If also μ ∈ Λ_{II} and 0 ≤ l ≤ d_{μ}, then [Y_{y,λ,k}, Y_{y,μ,l}]_{t} is by (9.13) given by a sum as in (9.30), where now by (2.8)

f_{1y}^{(i)}(s) = y^{−k−1/2}Q_{λ,k}P_{λ}e^{−sA} = y^{−k−1/2}e^{−sλ} ∑_{j=0}^{k} \frac{(−s)^{j}}{j!} Q_{λ,k}N_{λ}^{j}P_{λ}

and similarly for f_{2y}^{(i)}(s). Instead of (9.32) we have by (2.11), for t ≥ 1,

|f_{1t}^{(i)}(s)| = O(t^{−k−1/2}e^{−λ_{1}s/2}(1 + s)^{k}) = O(t^{−1/2}e^{−λ_{1}s/2}), 0 ≤ s ≤ t,

and similarly for f_{2y}^{(i)}(s). In this case too, (9.23)–(9.25) follow, so Lemma 9.9 applies for each i. We now have

\begin{matrix}∫_{0}^{t} f_{1t}^{(i)}(s) 𝔼 g(ξ_{i})f_{2t}^{(i)}(s)e^{λ_{1}s} ds \\ = ∫_{0}^{t} t^{−k−l−1}e^{s(λ_{1}−λ−μ)} ∑_{j=0}^{k} ∑_{m=0}^{l} \frac{(−s)^{j+m}}{j! m!} Q_{λ,k}N_{λ}^{j}P_{λ}B_{i}P′_{μ}(N′_{μ})^{m}Q′_{μ,l} ds. \end{matrix} ; (9.45)

If μ = λ̄, then λ_{1} − λ − μ = λ_{1} − 2 Re λ = 0, and the integral (9.45) equals

∑_{j=0}^{k} ∑_{m=0}^{l}(−1)^{j+m} \frac{t^{j+m+1−k−l−1}}{(j + m + 1) j! m!}Q_{λ,k}N_{λ}^{j}P_{λ}B_{i}P′_{μ}(N′_{μ})^{m}Q′_{μ,l}

→ \frac{(−1)^{k+l}}{(k + l + 1) k! l!}Q_{λ,k}N_{λ}^{k}P_{λ}B_{i}P′_{μ}(N′_{μ})^{l}Q′_{μ,l},

as t → ∞. If μ ≠ λ̄, then λ_{1} − λ − μ is imaginary and nonzero. It is easily seen, by integration by parts, that then ∫_{0}^{t} e^{s(λ_{1}−λ−μ)}s^{r} ds = O(t^{r}) for each r ≥ 0, which implies that the integral (9.45) tends to 0.

Hence, in both cases, Lemma 9.9 and (2.16) imply that, as t → ∞,

[Y_{t,λ,k}, Y_{t,μ,l}]_{t} →^{p} W · \frac{(−1)^{k+l}}{(k + l + 1) k! l!}Q_{λ,k}N_{λ}^{k}P_{λ}Σ_{II}P′_{μ}(N′_{μ})^{l}Q′_{μ,l}

(with the right-hand side 0 by (2.16) unless μ = λ̄). Consequently, for every fixed x > 0,

\begin{matrix}[Y_{t,λ,k}, Y_{t,μ,l}]_{xt} = x^{k+l+1}[Y_{xt,λ,k}, Y_{xt,μ,l}]_{xt} \\ →^{p} W · x^{k+l+1} \frac{(−1)^{k+l}}{(k + l + 1) k! l!}Q_{λ,k}N_{λ}^{k}P_{λ}Σ_{II}P′_{μ}(N′_{μ})^{l}Q′_{μ,l}. \end{matrix}  (9.46)

Clearly, this holds for x = 0 too. Moreover, it follows easily from (9.44), (9.17) and (2.11) that sup_{t≥1} 𝔼 |Y_{t,λ,k}(xt)|^{2} < ∞ for every x ≥ 0.

Let Y_{II,t}(s) be the vector obtained by combining all vectors Y_{t,λ,k}(s) for λ and k as in part (ii) of Theorem 3.1; we write this as Y_{II,t} = (Y_{t,λ,k})_{λ,k}. We thus see from (9.46) that [Y_{II,t}, Y_{II,t}]_{xt} converges. Moreover, since \overline{Y_{t,λ,k}} = Y_{t,λ̄,k}, \overline{Y_{II,t}} equals Y_{II,t} with a certain permutation of the components. Hence [Y_{II,t}, \overline{Y_{II,t}}]_{xt} too converges.

We this time define, with h := t^{1/2} as above,

Ỹ_{t}(x) := \begin{cases}(u_{1} · Y(h))^{−1/2}(Y_{II,t}(xt) − Y_{II,t}(h)), & x ≥ h/t, \\ 0, & x < h/t. \end{cases}

Repeating the argument after (9.37) above (replacing t − h by t/h = h and using Proposition 9.1(ii) if some λ ∈ Λ_{II} is non-real), we obtain, in D[0, ∞),

Y_{II,t}(xt) →^{d} W^{1/2}Ũ_{II}(x)  (9.47)

or, jointly for all λ ∈ Λ_{II} and 0 ≤ k ≤ d_{λ},

Y_{t,λ,k}(xt) →^{d} W^{1/2}Ũ_{λ,k}(x),  (9.48)

where Ũ_{II} = (Ũ_{λ,k})_{λ,k} is Gaussian and, for 0 ≤ x ≤ y,

𝔼(Ũ_{λ,k}(x)Ũ_{μ,l}(y)′) = c̃(k, l, x)Q_{λ,k}N_{λ}^{k}P_{λ}Σ_{II}P′_{μ}(N′_{μ})^{l}Q′_{μ,l}.  (9.49)

with

c̃(k, l, x) := \frac{(−1)^{k+l}}{(k + l + 1)k! l!}x^{k+l+1} = \frac{(−1)^{k+l}}{k! l!} ∫_{0}^{x} s^{k+l} ds.  (9.50)

In the case d_{λ} = 0 for λ ∈ Λ_{II}, this completes the proof of (ii). In general, we have by (9.9), (2.8), (2.6) and (9.44)

t^{−(k+1/2)}e^{−λxt}Q_{λ,k}P_{λ}X (xt) = t^{−k−1/2}e^{−λxt}Q_{λ,k}P_{λ}e^{xtA}Y(xt)

= t^{−k−1/2}Q_{λ,k}e^{xtN_{λ}}P_{λ}Y(xt) = t^{−k−1/2} ∑_{j=0}^{k} \frac{(xt)^{j}}{j!} Q_{λ,k}N_{λ}^{j}P_{λ}Y(xt)

= ∑_{j=0}^{k} \frac{x^{j}}{j!} N_{λ}^{j}Y_{t,λ,k−j}(xt),

which by (9.48) yields part (ii) of the theorem with

U_{λ,k}(x) := ∑_{j=0}^{k} \frac{x^{j}}{j!} N_{λ}^{j}Ũ_{λ,k−j}(x).

By (9.49) we have (3.2), with, using (9.50) and the binomial theorem,

c(k,l, x, y) = ∑_{j=0}^{k} ∑_{m=0}^{l} \frac{x^{j}}{j!} \frac{y^{m}}{m!} c̃(k − j, l − m, x)

= ∑_{j=0}^{k} ∑_{m=0}^{l} \frac{x^{j}}{j!} \frac{y^{m}}{m!} \frac{1}{(k − j)! (l − m)!} ∫_{0}^{x}(−s)^{k−j+l−m} ds

= \frac{1}{k! l!} ∫_{0}^{x}(x − s)^{k}(y − s)^{l} ds,

which yields (3.3) by a change of variable and the binomial theorem again.

Again, by Lemma 8.2 with an enlarged η_{t} as above, (ii) and (iii) hold jointly, with {U_{λ,k}}_{λ∈Λ_{II}, k≤d_{λ}} independent of {W_{λ,k}}_{λ∈Λ_{III},k≤d_{λ}}.

We have proved (i), (ii) and (iii) separately, but only partly proved the joint convergence and the asserted independence.

First, note that the families {U_{λ,k}}_{0≤k≤d_{λ}} for different λ ∈ Λ_{II} with Im λ ≥ 0 are independent because they are jointly Gaussian and all covariances of real or imaginary parts vanish by (3.2) and (2.16).

Next, consider the stopped processes Y_{II,t}^{T}(s) := Y_{II,t}(s∧T) and Ũ_{II}^{T}(x) := Ũ_{II}(x ∧ T). Stopping at x = 1 we obtain from (9.47)

Y_{II,t}^{t}(xt) →^{d} W^{1/2}Ũ_{II}^{1}(x).  (9.51)

Moreover, we have shown that this holds jointly with (iii) and with Ũ_{II}^{1} independent of {W_{λ,k}}_{λ∈Λ_{III},k≤d_{λ}}. It is another consequence of (9.47) that

sup_{x}|Y_{II,t}^{t}(xt) − Y_{II,t}^{t−h}(xt)| = sup_{1−h/t≤x≤1} |Y_{II,t}(xt) − Y_{II,t}(t − h)| →^{p} 0. (9.52)

Hence (9.51) is equivalent to

Y_{II,t}^{t−h}(xt) →^{d} W^{1/2}Ũ_{II}^{1}(x).  (9.53)

Next, redo the application above of Lemma 8.2 to (i) once again, now further enlarging η_{t} (and S_{1}) to contain also Y_{II,t}^{t−h}(xt) (which is F_{t−h}- measurable). This shows that (9.53), (iii) and (i) hold jointly, with U_{I}, Ũ_{II}^{1} and {W_{λ,k}}_{λ∈Λ_{III},k≤d_{λ}} independent.

Finally, stop the left-hand side of (i) at x = h; this does not affect the limit because h → ∞. Consider again the argument for (9.47), but applied to Y_{II,t} − Y_{II,t}^{t+h}. Applying Lemma 8.2 as above but with η_{t} containing Y_{II,t}^{t−h}, Y_{III,t}(h) and the left-hand side of (i) stopped at x = h, we find

Y_{II,t}(xt) − Y_{II,t}^{t+h}(xt) →^{d} W^{1/2}(Ũ_{II}(x) − Ũ_{II}^{1}(x)) in D[0, ∞), (9.54)

jointly with (i), (iii) and (9.53), with Ũ_{II} − Ũ_{II}^{1} independent of U_{I}, Ũ_{II}^{1} and {W_{λ,k}}_{λ∈Λ_{III},k≤d_{λ}}.

By (9.52) and the analogous sup_{x} |Y_{II,t}^{t+h}(xt) − Y_{II,t}^{t}(xt)| →^{p} 0, we can replace Y_{II,t}^{t±h} in (9.53) and (9.54) by Y_{II,t}^{t}. Together, these yield (9.47), and hence (ii), now with joint convergence with (i) and (iii) and the asserted independence. □

10. Remaining proofs

Proof of Lemma 3.9. By Lemma 9.8,

e^{−λ_{1}t}b · X (t) →^{a.s.} (b · v_{1})W  as t → ∞,  (10.1)

which shows that b · X (t) → ∞ a.s. when W > 0. The remaining properties are immediate, observing that sup_{t≤T} b · X (t) < ∞ for every finite T. □

The next lemma extends a result by Athreya and Karlin [7], [9, Theorem V.7.3].

Lemma 10.1. Assume b · v_{1} > 0. As z → ∞, (with ln 0 = −∞)

τ_{b}(z) − \frac{1}{λ_{1}} ln z →^{a.s.} −\frac{1}{λ_{1}} (ln W + ln(b · v_{1})).

Proof. Let E be the event b · X (t) → ∞. On E, τ_{b}(z) < ∞ for all z and τ_{b}(z) → ∞ as z → ∞. Hence (10.1) implies that, a.s. on E,

e^{−λ_{1}τ_{b}(z)}b · X (τ_{b}(z)) → (b · v_{1})W  as z → ∞.  (10.2)

By the right-continuity of X (t), b · X (τ_{b}(z)) ≥ z. On the other hand, if ε(z) = 1/z, say, and z is so large that τ_{b}(z) > ε(z), then b·X (τ_{b}(z)−ε(z)) < z, and, again by (10.1),

e^{−λ_{1}τ_{b}(z)}b · X (τ_{b}(z) − ε(z)) → (b · v_{1})W  as z → ∞.  (10.3)

Combining (10.2) and (10.3) we find that, a.s. on E,

e^{−λ_{1}τ_{b}(z)}z → (b · v_{1})W  as z → ∞.  (10.4)

On the complement of E, W = 0 a.s. by Lemma 3.9 and τ_{b}(z) = ∞ for large z; hence (10.4) holds trivially. In other words, (10.4) holds a.s. Taking logarithms, we obtain the lemma. □

Lemma 10.2. As t → ∞, e^{−λ_{1}t/2} sup_{s≤t} |ΔX (s)| →^{a.s.} 0.

Proof. Let Δ^{∗}X (t) := sup_{s≤t} |ΔX (s)| and let M > 0. Define

Z(t) := ∑_{i=1}^{q} ∑_{k:τ_{ik}≤t} 1[|ΔX (τ_{ik})| > M],  (10.5)

the number of jumps larger than M until time t. Clearly,

ℙ(Δ^{∗}X (t) > M) = ℙ(Z(t) ≥ 1) ≤ 𝔼 Z(t).

Moreover, applying Lemma 9.3 with f_{1}(s) = f_{2}(s) = 1 and g(x) = 1[|x| > M] to the inner sum in (10.5) and using (9.11),

𝔼 Z(t) = ∑_{i=1}^{q} ∫_{0}^{t} ℙ(|ξ_{i}| > M)a_{i} 𝔼 X_{i}(s) ds = ∑_{i=1}^{q} ℙ(|ξ_{i}| > M) · O(e^{λ_{1}t}).

Now let ε > 0 and let E_{n} be the event sup_{t∈[n,n+1]} e^{−λ_{1}t/2}Δ^{∗}X (t) > ε. Then, taking M = εe^{λ_{1}n/2},

\begin{aligned}ℙ(E_{n}) & ≤ ℙ(Δ^{∗}X (n + 1) > M) ≤ 𝔼 Z(n + 1) ≤ C ∑_{i=1}^{q} ℙ(|ξ_{i}| > M)e^{λ_{1}(n+1)} \\ & = C ∑_{i=1}^{q} ℙ(ε^{−2}|ξ_{i}|^{2} > e^{λ_{1}n})e^{λ_{1}(n+1)}. \end{aligned}

Summing over n we find

∑_{n} ℙ(E_{n}) ≤ C_{1} ∑_{i=1}^{q} 𝔼 ∑_{n} e^{λ_{1}(n+1)}1[e^{λ_{1}n} < ε^{−2}|ξ_{i}|^{2}] ≤ C_{2} ∑_{i=1}^{q} 𝔼 ε^{−2}|ξ_{i}|^{2} < ∞,

and the Borelli–Cantelli lemma completes the proof. □

Lemma 10.3. Conditioned on W > 0, z^{−1/2} sup_{s≤τ_{b}(z)} |ΔX (s)| →^{a.s.} 0 as z → ∞.

Proof. Conditioned on W > 0, we have τ_{b}(z) →^{a.s.} ∞ by Lemma 3.9 and thus e^{−λ_{1}τ_{b}(z)/2} sup_{s≤τ_{b}(z)} |ΔX (s)| →^{a.s.} 0 by Lemma 10.2. Moreover, by Lemma 10.1 (or (10.4)), z^{−1/2}e^{λ_{1}τ_{b}(z)/2} →^{a.s.} (b · v_{1})^{−1/2}W^{−1/2}, and the result follows by multiplication. □

Proof of Theorem 3.10. We condition on W > 0, which by Lemma 9.7 is the same as essential non-extinction, and let z → ∞.

First, τ_{b}(z) → ∞ and thus by Lemma 9.8 e^{−λ_{1}τ_{b}(z)}X (τ_{b}(z)) →^{a.s.} W v_{1}. Dividing by (10.4), we find (3.4).

Secondly, let t = λ_{1}^{−1} ln z and x = τ_{b}(z) − t. Then t → ∞ and, by Lemma 10.1, x → x_{0} := −λ_{1}^{−1}(ln W + ln(b · v_{1})) a.s. Thus Theorem 3.1(i) yields

e^{−λ_{1}t/2}P_{I}X (t + x) →^{d} e^{λ_{1}x_{0}/2}W^{1/2}U_{I}(x_{0}) = (b · v_{1})^{−1/2}U_{I}(x_{0}).

Here x_{0} is random, but U_{I} is a stationary process independent of W and thus of x_{0}, and thus U_{I}(x_{0}) =^{d} U_{I}(0). Hence (i) follows with V_{I} = U_{I}(0).

Next, let t = λ_{1}^{−1} ln z and x = τ_{b}(z)/t. Thus Lemma 10.1 yields xt − t → −λ_{1}^{−1}(ln W + ln(b · v_{1})) and x → 1 a.s. Hence Theorem 3.1(ii) implies for λ ∈ Λ_{II}, and thus Re λ = λ_{1}/2,

t^{−(k+1/2)}e^{−λ_{1}t/2−i Im λxt}Q_{λ,k}P_{λ}X (τ_{b}(z)) →^{d} (b · v_{1})^{−1/2}U_{λ,k}(1). (10.6)

If λ = λ_{1}/2, (10.6) yields the limit in (ii) with V_{λ,k} = λ_{1}^{−k−1/2}U_{λ,k}(1). If Im λ ≠ 0, we have a further factor e^{−i Im λxt} of modulus 1 on the left-hand side of (10.6). In this case, however, the distribution of U_{λ,k}(1) is symmetric complex Gaussian (Remark 3.6), and it follows that (ii) follows in this case too. (Consider subsequences where e^{−i Im λt} converges.)

Similarly, (iii) follows from Theorem 3.1(iii), with W̌_{λ,k} := λ_{1}^{−k}W^{−λ/λ_{1}}W_{λ,k}.

Next, |b·X (τ_{b}(z))−z| ≤ |b||ΔX (τ_{b}(z))| and thus z^{−1/2}(b·X (τ_{b}(z))−z) →^{a.s.}

0 by Lemma 10.3. We multiply this by (b · v_{1})^{−1}v_{1}, which yields (iv) if we write X = ∑_{λ} P_{λ}X and observe that for every v ∈ Im P_{λ_{1}} = E_{λ_{1}}, (b · v/b · v_{1})v_{1} = v. (v) is a consequence of (iv) and X = ∑_{λ} P_{λ}X.

The result for urn processes follows by the embedding argument by Athreya and Karlin [8] discussed in the introduction. □

Proof of Corollary 3.11. By assumption, Λ = Λ_{I} ∪ {λ_{1}}. Hence Theorem 3.10(v) yields, with some R(z) →^{a.s.} 0,

z^{−1/2}(X (τ_{b}(z)) − \frac{z}{b · v_{1}} v_{1}) = z^{−1/2}(I − \frac{v_{1}b′}{b · v_{1}} )P_{I}X (τ_{b}(z)) + R(z).

Theorem 3.10(i) shows that this converges to the Gaussian limit

(b · v_{1})^{−1/2}(I − \frac{v_{1}b′}{b · v_{1}} )V_{I}

with the covariance matrix (3.7), by (3.5) and (2.15); note that (I − (b · v_{1})^{−1}v_{1}b′)P_{I} = I − (b · v_{1})^{−1}v_{1}b′ because P_{I} = I − P_{λ_{1}} = I − v_{1}u′_{1} and (I − (b · v_{1})^{−1}v_{1}b′)v_{1} = 0. □

Proof of Corollary 3.12. This is similar to the proof of the preceding corollary; we now use Theorem 3.10(v), (i) and (ii), and find the Gaussian limit

(b · v_{1})^{−1/2}(I − \frac{v_{1}b′}{b · v_{1}} ) ∑_{λ∈Λ_{II}} V_{λ,d}

(where V_{λ,d} = 0 when d_{λ} < d), which yields (3.8) by (3.6) and (2.17). □

Proof of Corollary 3.13. By Theorem 3.10(v),

Y_{b}(z) − (I − \frac{v_{1}b′}{b · v_{1}} ) ∑_{λ≠λ_{1}}(ln z)^{−d}z^{− Re λ_{2}/λ_{1}}P_{λ}X (τ_{b}(z)) →^{a.s.} 0.

By Theorem 3.10(i), (ii) and (iii) (with k = d_{λ}), the terms in the sum with λ ∉ Λ̃′_{III} := {λ : Re λ = Re λ_{2}, d_{λ} = d} tend to 0 in probability; this can be improved to almost surely for λ ∈ Λ_{III} by Theorem 3.10(iii) and for λ ∈ Λ_{I} ∪ Λ_{II} by the argument in the proof of Lemma 9.8; we omit the details. Moreover, for λ ∈ Λ̃′_{III}, Theorem 3.10(iii) implies

(ln z)^{−d}z^{− Re λ_{2}/λ_{1}}P_{λ}X (τ_{b}(z)) − z^{i Im λ/λ_{1}}(b · v_{1})^{−λ/λ_{1}} W̌_{λ,k} →^{a.s.} 0.

Hence,

Y_{b}(z) − ∑_{λ∈Λ̃′_{III}} z^{i Im λ/λ_{1}} W̃_{b,λ} →^{a.s.} 0,

where W̃_{b,λ} := (I − \frac{v_{1}b′}{b·v_{1}} )(b·v_{1})^{−λ/λ_{1}} W̌_{λ,k}, and (3.9) follows with W_{b,λ} = W̃_{b,λ} when Im λ = 0 and W_{b,λ} = 2W̃_{b,λ} when Im λ > 0. Tightness follows, as does (3.10) in the special case when Λ_{III} = {λ_{2}} with λ_{2} real. □

Proof of Theorem 3.16. As said in Section 3, we change the rules and add an additional dummy ball of a new type q +1 whenever a ball splits; dummy balls have a_{q+1} = 0 and ξ_{q+1} = 0 and thus never split. This does not affect the process of the balls of types 1, …, q, but adds a count of the number of splits. Note that (A1)–(A6) hold for the modified process too. (However, it is not irreducible even if the original process is; this is our main reason for not assuming irreducibility in this paper.)

We write over the symbols to denote the modified process and various quantities defined for it. Writing vectors and matrices in block form, corresponding to a split ℝ^{q} × ℝ, we see that, for i = 1, …, q,

ξ̃_{i} = (\begin{matrix}ξ_{i} \\ 1 \end{matrix}) ,  ã = (\begin{matrix}a \\ 0 \end{matrix}) ,  Ã = (\begin{matrix}A & 0 \\ a′ & 0 \end{matrix}) .

Consequently, the eigenvalues Λ̃ = Λ ∪ {0}, in particular, λ̃_{1} = λ_{1}. It is easily verified that Ã has the corresponding eigenvectors

ũ_{1} = (\begin{matrix}u_{1} \\ 0 \end{matrix}) ,  ṽ_{1} = (\begin{matrix}v_{1} \\ λ_{1}^{−1} \end{matrix}) .  (10.7)

Finally, we choose b = ( \begin{matrix}0 \\ 1 \end{matrix} ), which means that τ_{b}(n) is the first time we have n dummy balls, i.e. the n:th split time τ_{n}, and thus X̃_{n} = X̃(τ_{b}(n)). By (10.7), we have b · ṽ_{1} = λ_{1}^{−1}.

We may now apply Theorem 3.10 to X̃_{n} = X̃(τ_{b}(n)), and the result follows from (3.4). □

Lemma 10.4. (i) We have ψ(s, A)v_{1} = v_{1}.

(ii) For s ≥ 1 and λ ∈ Λ,

ψ(s, A)P_{λ} = \begin{cases}P_{λ_{1}}, & λ = λ_{1}, \\ \frac{s^{d_{λ}}}{d_{λ}!} e^{λs}(I − λ^{−1}λ_{1}v_{1}a′)N_{λ}^{d_{λ}}P_{λ} + O(s^{d_{λ}−1}e^{λs}), & 0 < Re λ < λ_{1}, \\ O(s^{d_{λ}+1}), & Re λ ≤ 0. \end{cases}

Proof. By (3.11) and (3.12),

ϕ(s, A)v_{1} = ∫_{0}^{s} e^{tA}v_{1} dt = ∫_{0}^{s} e^{tλ_{1}}v_{1} dt = λ_{1}^{−1}(e^{sλ_{1}} − 1)v_{1},

and thus

ψ(s, A)v_{1} = e^{sλ_{1}}v_{1} − λ_{1}v_{1}a′λ_{1}^{−1}(e^{sλ_{1}} − 1)v_{1} = e^{sλ_{1}}v_{1} − v_{1}(e^{sλ_{1}} − 1) = v_{1}.

Hence, (2.7) implies ψ(s, A)P_{λ_{1}} = P_{λ_{1}}.

Similarly, by (2.8), when Re λ > 0 and s ≥ 1,

\begin{aligned}ϕ(s, A)P_{λ} & = ∫_{0}^{s} e^{tA}P_{λ} dt = ∑_{j=0}^{d_{λ}} \frac{N_{λ}^{j}}{j!} P_{λ} ∫_{0}^{s} t^{j}e^{λt} dt \\ & = (d_{λ}!)^{−1}λ^{−1}s^{d_{λ}}e^{λs}N_{λ}^{d_{λ}}P_{λ} + O(s^{d_{λ}−1}e^{λs}). \end{aligned}  (10.8)

For Re λ ≤ 0, the same argument yields ϕ(s, A)P_{λ} = O(s^{d_{λ}+1}). The result follows from this, (3.12) and (2.8). □

Proof of Theorem 3.17. We continue the argument from the proof of Theorem 3.16. We have, see (2.13) and (2.14), for i = 1, …, q,

B̃_{i} = 𝔼 (\begin{matrix}ξ_{i} \\ 1 \end{matrix}) (ξ′_{i} 1) = 𝔼 (\begin{matrix}ξ_{i}ξ′_{i} & ξ_{i} \\ ξ′_{i} & 1 \end{matrix}) = (\begin{matrix}B_{i} & 𝔼 ξ_{i} \\ 𝔼 ξ′_{i} & 1 \end{matrix}) ,

and thus, using (2.14), (2.18) and a · v = 1,

B̃ = ∑_{i=1}^{q+1} ṽ_{1i}ã_{i}B̃_{i} = ∑_{i=1}^{q} v_{1i}a_{i}B̃_{i} = (\begin{matrix}B & λ_{1}v_{1} \\ λ_{1}v′_{1} & 1 \end{matrix}) .  (10.9)

Further, Re λ̃_{2} = max(Re λ_{2}, 0) < \frac{1}{2} λ̃_{1}, so we may apply Corollary 3.11 to X̃_{n} = X̃(τ_{b}(n)), obtaining a Gaussian limit Ṽ_{b}. Ignoring the dummy balls,

we obtain n^{−1/2}(X_{n} − λ_{1}nv_{1}) →^{d} V, with V = (I, 0)Ṽ_{b}.

The covariance matrix of Ṽ_{b} is given by (3.7), with added everywhere. We have, by induction,

Ã^{n} = (\begin{matrix}A^{n} & 0 \\ a′A^{n−1} & 0 \end{matrix}) , n ≥ 1,

and thus

e^{sÃ} = ∑_{n=0}^{∞} \frac{s^{n}}{n!} Ã^{n} = (\begin{matrix}e^{sA} & 0 \\ a′ϕ(s, A) & 1 \end{matrix}) .  (10.10)

Hence

\begin{aligned}(I − \frac{ṽ_{1}b′}{b · ṽ_{1}} )e^{sÃ} & = ((\begin{matrix}I & 0 \\ 0 & 1 \end{matrix}) − λ_{1} (\begin{matrix}v_{1} \\ λ_{1}^{−1} \end{matrix}) (0 1)) e^{sÃ} \\ & = (\begin{matrix}I & −λ_{1}v_{1} \\ 0 & 0 \end{matrix}) (\begin{matrix}e^{sA} & 0 \\ a′ϕ(s, A) & 1 \end{matrix}) \\ & = (\begin{matrix}ψ(s, A) & −λ_{1}v_{1} \\ 0 & 0 \end{matrix}) . \end{aligned}  (10.11)

and thus, from (3.7), (10.9) and Lemma 10.4(i),

𝔼 V V′ = λ_{1} ∫_{0}^{∞} (ψ(s, A), −λ_{1}v_{1}) B̃ (ψ(s, A), −λ_{1}v_{1})′ e^{−λ_{1}s}ds
= λ_{1} ∫_{0}^{∞}(ψ(s, A)Bψ(s, A)′ − λ_{1}^{2}v_{1}v′_{1}ψ(s, A)′ − λ_{1}^{2}ψ(s, A)v_{1}v′_{1}
+ λ_{1}^{2}v_{1}v′_{1})e^{−λ_{1}s}ds
= λ_{1} ∫_{0}^{∞}(ψ(s, A)Bψ(s, A)′ − λ_{1}^{2}v_{1}v′_{1})e^{−λ_{1}s}ds,

which yields (3.13). □

We need another lemma, which expresses the covariance matrix in (3.8) as a limit of the integral in (3.7) over finite intervals, suitably renormalized.

Lemma 10.5. If Re λ_{2} = λ_{1}/2 and d := max{d_{λ} : Re λ = λ_{1}/2}, then, as t → ∞,

t^{−2d−1} ∫_{0}^{t}(I − P_{λ_{1}})e^{sA}Be^{sA′}(I − P′_{λ_{1}})e^{−λ_{1}s} ds → Σ_{II,d},  (10.12)

and thus, for any b with b · v_{1} ≠ 0,

t^{−2d−1} ∫_{0}^{t}(I − \frac{v_{1}b′}{b · v_{1}} )e^{sA}Be^{sA′}(I − \frac{bv′_{1}}{b · v_{1}} )e^{−λ_{1}s} ds

→ (I − \frac{v_{1}b′}{b · v_{1}} )Σ_{II,d}(I − \frac{bv′_{1}}{b · v_{1}} ). (10.13)

Proof. The left-hand side of (10.12) equals

\begin{aligned}∑_{λ,μ≠λ_{1}} & t^{−2d−1} ∫_{0}^{t} P_{λ}e^{sλ+sN_{λ}}Be^{sμ+sN′_{μ}}P′_{μ}e^{−λ_{1}s}ds \\ & = ∑_{λ,μ≠λ_{1}} t^{−2d−1} ∑_{k=0}^{d} ∑_{l=0}^{d} ∫_{0}^{t} \frac{s^{k+l}}{k! l!} e^{(λ+μ−λ_{1})s}ds P_{λ}N_{λ}^{k}B(N′_{μ})^{l}P′_{μ}. \end{aligned}

Here Re λ, Re μ ≤ λ_{1}/2. For such λ and μ, the integral is O(t^{2d}) unless Re λ = λ_{1}/2, μ = λ̄ and k = l = d, cf. the argument after (9.45). In the remaining case, the integral is (2d + 1)^{−1}d!^{−2}t^{2d+1}, and (10.12) follows.

Multiplying (10.12) to the left by (I − (b · v_{1})^{−1}v_{1}b′) and to the right by the transpose, we obtain (10.13) because (I − (b · v_{1})^{−1}v_{1}b′)P_{λ_{1}} = 0. □

Proof of Theorem 3.18. We argue as in the proof of Theorem 3.17, now applying Corollary 3.12. We have Re λ̃_{2} = Re λ_{2} = \frac{1}{2} λ̃_{1}. Moreover, it follows from (2.9) that |e^{sA}(I − P_{λ_{1}})| grows as s^{d}e^{Re λ_{2}s} as s → ∞; it follows easily from (10.10) and (3.11) that |e^{sÃ}(I − P̃_{λ_{1}})| grows at the same rate, and thus d̃= d.

Hence, it follows from Corollary 3.12 that (n ln^{2d+1} n)^{−1/2}(X_{n}−nλ_{1}v_{1}) →^{d} V := (I, 0)Ṽ_{b}, where Ṽ_{b} is Gaussian and, by (3.8) and Lemma 10.5,

\begin{aligned}𝔼 Ṽ_{b}Ṽ′_{b} & = λ_{1}^{−2d}(I − \frac{ṽ_{1}b′}{b · ṽ_{1}} )Σ̃_{II,d}(I − \frac{bv′_{1}}{b · ṽ_{1}} ) \\ & = λ_{1}^{−2d} lim_{t→∞} t^{−2d−1} ∫_{0}^{t}(I − \frac{ṽ_{1}b′}{b · ṽ_{1}} )e^{sÃ}B̃e^{sÃ′}(I − \frac{bv′_{1}}{b · ṽ_{1}} )e^{−λ_{1}s} ds. \end{aligned} ; (10.14)

This implies, using (10.11), (10.9) and Lemma 10.4(i),

\begin{aligned}𝔼 V V′ & = λ_{1}^{−2d} lim_{t→∞} t^{−2d−1} ∫_{0}^{t} (ψ(s, A), −λ_{1}v_{1}) B̃ (ψ(s, A), −λ_{1}v_{1})′ e^{−λ_{1}s}ds \\ & = λ_{1}^{−2d} lim_{t→∞} t^{−2d−1} ∫_{0}^{t}(ψ(s, A)Bψ(s, A)′ − λ_{1}^{2}v_{1}v′_{1})e^{−λ_{1}s}ds. \end{aligned}

We write ψ(s, A) = ∑_{λ} ψ(s, A)P_{λ} and use Lemma 10.4(ii); (3.14) follows. For the final claim, note that

T P_{λ}N_{λ}^{d} = λ^{−1}λ_{1}v_{1}a′P_{λ}N_{λ}^{d} = λ^{−1}λ_{1}v_{1}â′AP_{λ}N_{λ}^{d} = T_{1}P_{λ}N_{λ}^{d}  □

Proof of Theorem 3.19. Follows directly by applying Corollary 3.13 to X̃_{n}. □

Proof of Theorem 3.21. We add dummy balls with a_{i} = 0, i = q + 1, …, 2q, as described in Section 3, and argue as in the proofs of Theorems 3.16–3.18. For the modified process we now have, with vectors and matrices in block form corresponding to a split ℝ^{q} × ℝ^{q},

X̃_{n} = (\begin{matrix}X_{n} \\ N_{n} \end{matrix}) ,  ξ̃_{i} = (\begin{matrix}ξ_{i} \\ δ_{i} \end{matrix}) ,  ã = (\begin{matrix}a \\ 0 \end{matrix}) ,  Ã = (\begin{matrix}A & 0 \\ D_{a} & 0 \end{matrix}) ,

where (δ_{i})_{j} = δ_{ij} and (D_{a})_{ij} = a_{i}δ_{ij}. Again, the non-zero eigenvalues are the same: Λ̃ = Λ ∪ {0}. Ã has the eigenvectors

ũ_{1} = (\begin{matrix}u_{1} \\ 0 \end{matrix}) ,  ṽ_{1} = (\begin{matrix}v_{1} \\ λ_{1}^{−1}D_{a}v_{1} \end{matrix}) .

We take b = ( \begin{matrix}0 \\ 1 \end{matrix} ), where 1 denotes the vector (1, …, 1)′, and note that b · ṽ_{1} = λ_{1}^{−1}a · v_{1} = λ_{1}^{−1}. The a.s. convergence follows from (3.4) and the Gaussian limits in (i) and (ii) follow by Corollaries 3.11 and 3.12.

For the explicit forms of the covariance matrices in (i) we use (3.7) and compute (omitting the details)

e^{sÃ} = ∑_{n=0}^{∞} \frac{s^{n}}{n!} Ã^{n} = (\begin{matrix}e^{sA} & 0 \\ D_{a}ϕ(s, A) & I \end{matrix}) ,

(I − \frac{ṽ_{1}b′}{b · ṽ_{1}} )e^{sÃ} = (\begin{matrix}ψ(s, A) & −λ_{1}v_{1}1′ \\ D_{a}(I − v_{1}a′)ϕ(s, A) & I − D_{a}v_{1}1′ \end{matrix}) ,  (10.15)
B̃ = (\begin{matrix}B & AD_{v} \\ D_{v}A′ & D_{a}D_{v} \end{matrix}) .  (10.16)

Assuming a = 1, and thus D_{a} = I, the result now follows from (3.7) by simple calculations, using D_{v}1 = v_{1} and ϕ(s, A)A = e^{sA} − I.

For (ii), we use (3.8) and Lemma 10.5 and obtain (10.14) as in the proof of Theorem 3.18. We extract the leading terms in the integral in (10.14) using (10.15) and (10.16) together with Lemma 10.4 and (10.8), cf. the proof of Lemma 10.5. In the case a = 1, this yields, after some calculations, with T as in (3.14),

𝔼 (V ; V̂) (V′ V̂′) = λ_{1}^{−2d} (I − T ; T̂ ) Σ_{II,d} (I − T′ T̂′) .

This yields the (co)variances. Since AT̂P_{λ}N_{λ}^{d} = (I − T)P_{λ}N_{λ}^{d}, λ ∈ Λ_{II}, it also implies that 𝔼(V − AV̂)(V − AV̂)′ = 0, and thus V = AV̂ a.s. □

Lemma 10.6. For all t ≥ 0,

𝔼 sup_{s≤t} |P_{I}X (s)|^{2} ≤ Ce^{λ_{1}t}

and, if λ ∈ Λ_{II} and k ≥ 0,

𝔼 sup_{s≤t} |Q_{λ,k}P_{λ}X (s)|^{2} ≤ C(t + 1)^{2k+1}e^{λ_{1}t}.

Proof. It follows easily from (9.18), (9.11), (2.10) and (2.11) that, for λ ∈ Λ_{I},

𝔼 |P_{λ}X (t)|^{2} = O(e^{λ_{1}t})

and, for λ ∈ Λ_{II},

𝔼 |Q_{λ,k}P_{λ}X (t)|^{2} = O((t + 1)^{2k+1}e^{λ_{1}t}).

For each u, e^{A(u−t)}P_{λ}X (t) = e^{Au}P_{λ}Y(t) is a martingale, and thus by Doob’s inequality

𝔼 sup_{u−1≤s≤u} |P_{λ}X (s)|^{2} ≤ C_{1} 𝔼 sup_{u−1≤s≤u} |e^{A(u−s)}P_{λ}X (s)|^{2} ≤ C_{2} 𝔼 |P_{λ}X (u)|^{2},

and similarly for Q_{λ,k}P_{λ}X, and the results follows by summing over all integers u less than t + 1. □

Proof of Theorem 3.22. We condition on W > 0.

(i): Let t = λ_{1}^{−1} ln z. By Lemma 10.1,

τ_{b}(xz) − t → ϕ(x) := λ_{1}^{−1}(ln x − ln W − ln(b · v_{1})) a.s.,

uniformly for x ∈ K for any compact interval K ⊂ (0, ∞). This together with Theorem 3.1(i) yields, in D(0, ∞),

e^{−λ_{1}t/2}P_{I}X (τ_{b}(xz)) →^{d} W^{1/2}e^{λ_{1}ϕ(x)/2}U_{I}(ϕ(x)) = x^{1/2}(b · v_{1})^{−1/2}U_{I}(ϕ(x)).

Since U_{I} is translation invariant and independent of W, the processes U_{I}(ϕ(x)) and U_{I}(ϕ(x) − ϕ(1)) = U_{I}(λ_{1}^{−1} ln x) have the same distribution, and thus

z^{−1/2}P_{I}X (τ_{b}(xz)) = e^{−λ_{1}t/2}P_{I}X (τ_{b}(xz)) →^{d} (b · v_{1})^{−1/2}V_{I}(x) in D(0, ∞),

with V_{I}(x) := x^{1/2}U_{I}(λ_{1}^{−1} ln x). V_{I} is a continuous Gaussian process because U_{I} is, and (3.1) implies (3.17).

It remains to improve the result to convergence in D[0, ∞), with V_{I}(0) = 0. For this, it suffices to show that for every ε > 0,

lim sup_{z→∞} ℙ( sup_{0≤x≤h} z^{−1/2}|P_{I}X (τ_{b}(xz))| > ε) → 0

as h → 0, see e.g. [27, Proposition 2.4], and this is an easy consequence of Lemmas 10.6 and 10.1. Continuity of V_{I} at 0 follows, or is proved directly by standard methods [31, Theorem 3.23].

(ii): Convergence in D(0, ∞) follows as in (i), now using Theorem 3.1(ii) and Lemma 10.1, observing as in the proof of Theorem 3.10(ii) that the constant factor (b·v_{1}W)^{i Im λ/λ_{1}} of modulus one can be ignored. The convergence extends to D[0, ∞) as in part (i), using Lemma 10.6(ii). □

Proof of Corollary 3.23. A simple consequence of Theorems 3.22 and 3.10(v); (3.18) follows from (3.2) and (2.16)–(2.17). □

Proof of Theorem 3.24. This follows from Corollary 3.23 by the arguments in the proofs of Theorems 3.17 and 3.18. In particular, the (co)variances in (3.19) follow easily using (10.11) and a change of variables; for (3.20) we observe that Corollary 3.23 shows that the covariance matrix depends on x and y only through the numerical factor c(d, d, x, y), and the result follows by comparison with (3.14) (the case x = y = 1), where now T = 2v_{1}a′P_{λ_{2}}, together with (3.3). □

References

[1] D. Aldous, Asymptotic fringe distributions for general families of random trees. Ann. Appl. Probab. 1 (1991), no. 2, 228–266.

[2] D. Aldous, B. Flannery & J.L. Palacios, Two applications of urn processes. Probab. Engineering Inform. Sci. 2 (1988), 293–307.

[3] S. Asmussen, Almost sure behavior of linear functionals of supercritical branching processes. Trans. Amer. Math. Soc. 231 (1977), no. 1, 233–248.

[4] K.B. Athreya, Some results on multitype continuous time Markov branching processes. Ann. Math. Statist. 39 (1968), 347–357.

[5] K.B. Athreya, Limit theorems for multitype continuous time Markov branching processes. I. The case of an eigenvector linear functional. Z. Wahrsch. Verw. Gebiete 12 (1969), 320–332.

[6] K.B. Athreya, Limit theorems for multitype continuous time Markov branching processes. II. The case of an arbitrary linear functional. Z. Wahrsch. Verw. Gebiete 13 (1969), 204–214.

[7] K.B. Athreya & S. Karlin, Limit theorems for the split times of branching processes. J. Math. Mech. 17 (1967), 257–277.

[8] K.B. Athreya & S. Karlin, Embedding of urn schemes into continuous time Markov branching processes and related limit theorems. Ann. Math. Statist. 39 (1968), 1801– 1817.

[9] K.B. Athreya & P.E. Ney, Branching Processes. Springer-Verlag, Berlin, 1972.

[10] A. Bagchi & A.K. Pal, Asymptotic normality in the generalized Pólya–Eggenberger urn model, with an application to computer data structures. SIAM J. Algebraic Discrete Methods 6 (1985), no. 3, 394–405.

[11] Z.D. Bai & F. Hu, Asymptotic theorems for urn models with nonhomogeneous generating matrices. Stochastic Process. Appl. 80 (1999), no. 1, 87–101.

[12] Z.D. Bai, F. Hu & L.-X. Zhang, Gaussian approximation theorems for urn models and their applications. Ann. Appl. Probab. 12 (2002), no. 4, 1149–1173.

[13] S. Bernstein, Nouvelles applications des grandeurs aléatoires presqu’indépendantes. (Russian) Izv. Akad. Nauk SSSR Ser. Mat. 4 (1940), 137–150.

[14] S. Bernstein, Sur un problème du schéma des urnes à composition variable. C. R. (Doklady) Acad. Sci. URSS (N.S.) 28 (1940), 5–7.

[15] P. Billingsley, Convergence of Probability Measures. Wiley, New York, 1968.

[16] H.-H. Chern & H.-K. Hwang, Phase changes in random m-ary search trees and generalized quicksort. Random Structures Algorithms 19 (2001), no. 3-4, 316–358.

[17] F. Eggenberger & G. Pólya, Über die Statistik verketteter Vorgänge. Zeitschrift Angew. Math. Mech. 3 (1923), 279–289.

[18] D.A. Freedman, Bernard Friedman’s urn. Ann. Math. Statist. 36 (1965), 956–970.

[19] B. Friedman, A simple urn model. Comm. Pure Appl. Math. 2 (1949), 59–70.

[20] R. Gouet, Martingale functional central limit theorems for a generalized Pólya urn. Ann. Probab. 21 (1993), no. 3, 1624–1639.

[21] A. Gut, Stopped Random Walks. Springer-Verlag, New York, 1988.

[22] A. Gut & S. Janson, The limiting behaviour of certain stopped sums and some applications. Scand. J. Statistics 10 (1983), 281–292.

[23] H.-K. Hwang, Second phase changes in random m-ary search trees and generalized quicksort: convergence rates. Preprint, 2001.

[24] J. Jacod & A.N. Shiryaev, Limit Theorems for Stochastic Processes. Springer-Verlag, Berlin, 1987.

[25] S. Janson, Limit theorems for certain branching random walks on compact groups and homogeneous spaces. Ann. Probab. 11 (1983), 909-930.

[26] S. Janson, A functional limit theorem for random graphs with applications to subgraph count statistics, Random Structures Algorithms 1 (1990), 15–37.

[27] S. Janson, Orthogonal Decompositions and Functional Limit Theorems for Random Graph Statistics. Mem. Amer. Math. Soc., vol. 111, no. 534, American Mathematical Society, Providence, R.I., 1994.

[28] S. Janson, Gaussian Hilbert Spaces. Cambridge Univ. Press, Cambridge, 1997.

[29] S. Janson, Asymptotic degree distribution in random recursive trees. To appear.

   Available from http://www.math.uu.se/~svante/

[30] N.L. Johnson & S. Kotz, Urn models and their application. Wiley, New York, 1977.

[31] O. Kallenberg, Foundations of modern probability. 2nd ed., Springer-Verlag, New York, 2002.

[32] S. Karlin, A First Course in Stochastic Processes. Academic Press, New York, 1969.

[33] H. Kesten & B.P. Stigum, Additional limit theorems for indecomposable multidimensional Galton–Watson processes. Ann. Math. Statist. 37 (1966), 1463–1481.

[34] H. Kesten & B.P. Stigum, Limit theorems for decomposable multi-dimensional Galton–Watson processes. J. Math. Anal. Appl. 17 (1967), 309–338.

[35] S. Kotz, H. Mahmoud & P. Robert, On generalized Pólya urn models. Statist. Probab. Lett. 49 (2000), no. 2, 163–173.

[36] W. Lew & H.M. Mahmoud, The joint distribution of elastic buckets in multiway search trees. SIAM J. Comput. 23 (1994), no. 5, 1050–1074.

[37] T. Lindvall, Weak convergence of probability measures and random functions in the function space D(0, ∞). J. Appl. Probab. 10 (1973), 109–121.

[38] H.M. Mahmoud, Evolution of Random Search Trees. Wiley, New York, 1992.

[39] H.M. Mahmoud, On rotations in fringe-balanced binary trees. Inform. Process. Lett. 65 (1998), no. 1, 41–46.

[40] H.M. Mahmoud, Urn models evolving by drawing multisets of balls. Preprint, 2000.

[41] H.M. Mahmoud, The size of random bucket trees via urn models. Acta Inform. 38 (2002), no. 11-12, 813–838.

[42] H.M. Mahmoud & B. Pittel, Analysis of the space of search trees under the random insertion algorithm. J. Algorithms 10 (1989), no. 1, 52–75.

[43] H.M. Mahmoud & R.T. Smythe, Asymptotic joint normality of outdegrees of nodes in random recursive trees. Random Structures Algorithms 3 (1992), no. 3, 255–266.

[44] H.M. Mahmoud, R.T. Smythe & J. Szymański, On the structure of random planeoriented recursive trees and their branches. Random Structures Algorithms 4 (1993), no. 2, 151–176.

[45] A. Meir & J.W. Moon, Recursive trees with no nodes of out-degree one. Congr. Numer. 66 (1988), 49–62.

[46] K. Nomizu, Fundamentals of linear algebra. 2nd ed. Chelsea Publishing Co., New York, 1979.

[47] G. Pólya, Sur quelques points de la théorie des probabilités. Ann. Inst. Poincaré 1 (1931), 117–161.

[48] P. Protter, Stochastic integration and differential equations. A new approach. Applications of Mathematics, 21, Springer-Verlag, Berlin, 1990.

[49] V. Savkevitch, Sur le schéma des urnes à composition variable. C. R. (Doklady) Acad. Sci. URSS (N.S.) 28 (1940), 8–12.

[50] E. Seneta, Nonnegative matrices and Markov chains. 2nd ed. Springer-Verlag, New York, 1981.

[51] R.T. Smythe, Central limit theorems for urn models. Stochastic Process. Appl. 65 (1996), no. 1, 115–137.

[52] R.T. Smythe & W.F. Rosenberger, Play-the-winner designs, generalized Pólya urns, and Markov branching processes. In Adaptive designs, eds. N Flournoy & W.F. Rosenberger, IMS Lecture Notes Monogr. Ser. 25, Inst. Math. Statist., Hayward, CA, 1995, 13–22.

[53] L.J. Wei & S. Durham, The randomized play-the-winner rule in medical trials. J. Amer. Statist. Assoc. 73 (1978), issue 364, 840–843.

[54] L.J. Wei, R.T. Smythe, D.Y. Lin & T.S. Park, Statistical inference with datadependent treatment allocation rules. J. Amer. Statist. Assoc. 85 (1990), issue 409, 156–162.

Department of Mathematics, Uppsala University, PO Box 480, S-751 06 Uppsala, Sweden

E-mail address: svante.janson@math.uu.se
URL: http://www.math.uu.se/~svante/