Cyril Banderier

Markus Kuba

Michael Wallner

Phase transitions of composition schemes: Mittag-Leffler and mixed Poisson distributions

Keywords: Primary 60C05; secondary 60E07, 60G50, 05A15 Mixed Poisson distributions, Mittag-Leffler distributions, stable laws, Boltzmann distributions, analytic combinatorics, generating functions, singularity analysis, critical composition schemes

1. Introduction. Many combinatorial structures are an assemblage of more basic building blocks, and this situation is ubiquitous in many different fields, such as combinatorics, probability theory, and statistical mechanics. It appears for example in permutations, random walks, random mappings, random forests, parking functions, Pólya-Eggenberger urn models, (Bienaymé)-Galton-Watson processes, destruction procedures in simply generated trees, inversions in labelled tree families, generalized plane-oriented recursive trees (scale-free trees), set or integer partitions with some constraints, sequences of words, tilings, different families of graphs or maps, etc.; see, e.g., [5-7, 11-13, 28, 29, 37, 40-42, 49, 83, 99, 100, 102, 112, 118]. In the language of generating functions, one then has a functional composition scheme such as

F (z) = G H(z) .
Let us illustrate this composition scheme with some examples (each of them being in fact the starting point of many theorems in the literature): A random forest is a set of random trees; a permutation is a set of cycles; a bridge (a random walk on Z) is a sequence of arches; functional mappings are cycles of Cayley trees; supertrees are trees in which each leaf is replaced by another family of trees; an integer partition is a sequence of parts; the factorization of a polynomial in a finite field is a multiset of irreducible factors; following the work of Tutte, several important families of planar maps can be seen as a "simple core" in which each node is replaced by some "simple map", etc. The reader can find many other examples illustrating the universality of the scheme F (z) = G H(z) in the wonderful book by Flajolet and Sedgewick on analytic combinatorics [START_REF] Flajolet | Analytic Combinatorics[END_REF]. Structurally, this composition scheme is at the heart of many fascinating phase transition phenomena (analytically corresponding, e.g., to coalescing saddle points or to confluence of singularities). More precisely, let G(z) = n≥0 g n z n and H(z) = n≥0 h n z n be analytic functions at the origin with nonnegative coefficients and H(0) = 0. Let ρ G and ρ H be the radii of convergence of G(z) and H(z), respectively. Then, following [START_REF] Banderier | Random maps, coalescing saddles, singularity analysis, and Airy phenomena[END_REF][START_REF] Flajolet | Analytic Combinatorics[END_REF], we focus on critical composition schemes. DEFINITION 1.1 (Critical composition scheme). The composition scheme

F (z) = G H(z) is critical if it satisfies H(ρ H ) = ρ G .
In other words, G(z) and H(z) are concomitantly singular. We will assume throughout this work that we are always in the critical case (the asymptotic analysis is straightforward otherwise). Note that this terminology is a generalization of the notion of critical/supercritical/subcritical Galton-Watson processes, initially popularized by Harris for neutron branching processes [START_REF] Edward | The Theory of Branching Processes[END_REF].

Often, G(z) and H(z) are the counting series of certain combinatorial families G and H such that F = G(H). We refer to the first part of [START_REF] Flajolet | Analytic Combinatorics[END_REF] for a more detailed presentation of this combinatorial approach, starting from the so-called atoms, and then assembling them into more elaborate structured blocks via combinatorial constructors. Some important subclasses of such structures were also subject of more probabilistic approaches; see, e.g., [START_REF] Arratia | Logarithmic Combinatorial Structures: a Probabilistic Approach[END_REF][START_REF] Goldschmidt | Stable graphs: distributions and line-breaking construction[END_REF][START_REF] Valentin | Random Graphs[END_REF]. Now, our goal is to analyse probabilistic properties of critical compositions like

F (z, u) = G uH(z) ,
where u marks each occurrence of objects of H. From a combinatorial perspective, F (z, u) enumerates F -structures of size n made of k "building blocks from H" (also simply called H-components), i.e., G-structures made of k blocks, where each block is then replaced by an H-block (which is itself a structured set of atoms). For any combinatorial structure in F , its corresponding G-structure is sometimes called its core1 (or "skeleton", or "backbone").

A first natural question is what is the typical size of this core, i.e., what is the typical number of H-components? Such insight helps, for example, to make many algorithms on combinatorial structures more efficient, as Knuth shows in [START_REF] Knuth | Sorting and searching[END_REF] (see also [START_REF] Banderier | Random maps, coalescing saddles, singularity analysis, and Airy phenomena[END_REF], where this insight is used to design faster random generation algorithms). To answer this question, one considers the discrete random variable X n associated with this core size in a uniformly chosen object of size n. Its probability mass function is obtained by extraction of coefficients:

P{X n = k} = [z n u k ]F (z, u) [z n ]F (z, 1) = g k [z n ]H(z) k f n ,
where [z n ] denotes the extraction of coefficient operator:

[z n ] n f n z n = f n .
As H(z) has typically a singular expansion of the type

H(z) = τ H + c H 1 - z ρ H λH + . . . ,
this implies that the asymptotic behaviour of P{X n = k} depends on the exponent λ H (which is called the singular exponent of H). Actually, this exponent even plays a key rôle, as it entails four types of asymptotic behaviour for X n :

• For λ H > 2 the limit law is related to a Gaussian law.

• For 1 < λ H < 2 the limit law is related to a stable law of parameter λ H (this distribution is supported on R and possesses moments up to order λ H ; e.g., for λ H = 3/2 this gives the map-Airy distribution). • For 0 < λ H < 1 we will show that the limit law is related to a stable law of parameter λ H , or more precisely to a generalized Mittag-Leffler distribution (this distribution is supported on R + and has moments of any order; e.g., for λ = 1/2, this gives the Rayleigh distribution). • For λ H < 0 the scheme is not critical because the function H(z) diverges at z = ρ H , and thus leads to a singularity of G(H(z)) at some z < ρ H . Such a scheme is called supercritical and typically leads to a Gaussian limit law.

The case λ H < 0 is analysed by Flajolet and Sedgewick [START_REF] Flajolet | Analytic Combinatorics[END_REF], building on the seminal work of Bender [START_REF] Bender | Central and local limit theorems applied to asymptotic enumeration[END_REF]. The three cases 0 < λ H < 1, 1 < λ H < 2, and λ H > 2 were partially analysed by Flajolet et al. [START_REF] Banderier | Random maps, coalescing saddles, singularity analysis, and Airy phenomena[END_REF], but without a precise statement for the limit laws, the right renormalizations, etc. It is partly due to the fact that the initial motivation of the authors of [START_REF] Banderier | Random maps, coalescing saddles, singularity analysis, and Airy phenomena[END_REF] was to analyse the core of planar maps, so they focused on the subcase λ H = 3 2 , which corresponds to the map-Airy distribution. Thus, a more complete analysis of the composition scheme in these three regions of λ H remained to be done.

For sure, we expected that the different possible analytic behaviours of G introduce further subcases, but we were surprised that the detailed analyses were much more challenging than expected: As we shall see, they require several new ingredients. Our identification of the limit laws involves moment-tilted distributions, product distributions, and Boltzmann distributions (see Section 3.2 for a formal presentation of these three types of distributions and their key properties). Therefore, the first main objective of our work is to give a complete landscape of the limit laws associated with critical composition schemes. We analyse the case 0 < λ H < 1 in this article, and the cases 1 < λ H < 2 and λ H > 2 in our companion article [START_REF] Banderier | Phase transitions of composition schemes: stable laws and Boltzmann distributions[END_REF].

Our second main objective is to explain the phase transitions observed for the number of H-components of a given size. This builds on the work of Panholzer and the second author [START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF], in which they started to unify the diversity of limit laws encountered in these phase transitions under the umbrella of mixed Poisson distributions, and relies on the study of a size-refined composition scheme that we detail in the next section.

2. New main results.

2.1.

Composition schemes analysed in this article. In this work we complete the analysis of composition schemes with exponent 0 < λ H < 1. We identify the corresponding limit laws as generalized Mittag-Leffler distributions and product distributions; this unifies and refines many previous studies. Motivated by models associated with directed lattice paths [START_REF] Banderier | Basic analytic combinatorics of directed lattice paths[END_REF][START_REF] Banderier | Some reflections on directed lattice paths[END_REF][START_REF] Banderier | Lattice paths with catastrophes[END_REF][START_REF] Banderier | Local time for lattice paths and the associated limit laws[END_REF][START_REF] Wallner | A half-normal distribution scheme for generating functions[END_REF] and triangular Pólya urns [START_REF] Flajolet | Some exactly solvable models of urn process theory[END_REF][START_REF] Janson | Limit theorems for triangular urn schemes[END_REF][START_REF] Janson | Moments of Gamma type and the Brownian supremum process area[END_REF], we also relax and extend the scheme by adding another component M (z), which allows us to apply our results to various examples in Section 6.

First, we consider the following extended composition scheme

F (z) = G H(z) • M (z), (1)
for some functions F /G/H/M analytic at the origin, with nonnegative coefficients. Such schemes are critical if ρ G = H(ρ H ) (like in Definition 1.1) and additionally satisfy ρ M ≥ ρ H , where ρ M is the radius of convergence of M (z) (the analysis is straightforward if the extended composition scheme is not critical). In order to enumerate the family F according to the occurrences of H-components, we consider

F (z, u) = G uH(z) • M (z), (2)
which from now on we will refer to as the extended composition scheme. Equivalently, [z n u k ]F (z, u) is the number of F -structures of size n having k H-components. The corresponding random variable X n has a probability mass function given by

P{X n = k} = [z n u k ]F (z, u) [z n ]F (z, 1) = g k [z n ]H(z) k • M (z) f n . (3)
Our first main result is Theorem 4.1, in which we give explicit expressions for the asymptotics of the factorial moments of X n , the limit distribution of X n (suitably normalized), and its density function. It appears that this limit distribution differs depending on some relationship between the singular exponents of G(z), H(z), and M (z), as summarized in Table 1 (where we write X n ∼ c n • D when c -1 n • X n → D in distribution for n → ∞). In addition to these convergences in distribution, we prove moment convergence for the continuous limit laws.

Singular λ M > λ G λ H λ M = λ G λ H λ M < λ G λ H exponent
(pure scheme) (confluent scheme) (degenerate scheme)

Limit law continuous linear combination discrete (gen. Mittag-Leffler ML)

(ML + B) (Boltzmann B) Example Xn ∼ Cn λ H ML Xn ∼ LinComb(n λ H ML, B) P{Xn = k} ∼ g k ρ k G G(ρ G )
TABLE 1 Three asymptotic behaviours according to the singular exponents λ G /λ H /λ M of G/H/M : the number Xn of H-components in F -structures of size n in the extended scheme (1) is given by three completely different types of limit laws, depending on whether the scheme is analytically pure/confluent/degenerate (Theorems 4. 1, 4.3, and 4.4).

Secondly, we consider a size-refined composition scheme which allows us to capture some threshold phenomenon via a bivariate generating function F (z, v) 2 :

F (z, v) = G H(z) -(1 -v)h j z j • M (z), j ∈ N. (4)
In this scheme, [z n v k ]F (z, v) is therefore the number of F -structures of size n having k H-components each of size j; this is combinatorially summarized by:

F = G H ̸ =j + vH =j × M.
Given n, j ∈ N, let X n,j denote the random variable counting the number of H-components of size j inside F -structures3 of size n. Note that the random variables X n,j naturally refine the distribution of the core size X n given in (3), since j∈N X n,j = X n .

Formula (4) implies that one has

P{X n,j = k} = [z n v k ]F (z, v) [z n ]F (z, 1) = h k j k! [z n-jk ]G (k) H(z) -h j z j M (z) f n . (5)
Our second main result is Theorem 5.1, in which we prove that the factorial moments of X n,j are asymptotically of mixed Poisson type, and establish a convergence in distribution, with convergence of all moments, towards explicit limit distributions.

We extend these two main results to the composition schemes involving a logarithmic singularity in Theorem 7.1 and 7.3, leading to Mittag-Leffler distributions.

Then, our third main result is Theorem 7.4 in which we give a multivariate generalization with arbitrary many variables, leading to Dirichlet product distributions.

Phase transitions.

Analogously to what can happen in physics or chemistry for some small change of temperature, pressure, or concentration, a phase transition in mathematics corresponds to a sudden non-smooth change of properties under smooth variation of the parameters. Such non-smooth changes are thus analytically reflected by a singularity of some function associated with these properties. For combinatorial structures, generating function methods were successfully used to analyse such phase transitions, e.g., for random graphs or planar maps [START_REF] Banderier | Random maps, coalescing saddles, singularity analysis, and Airy phenomena[END_REF][START_REF] Giménez | Asymptotic enumeration and limit laws of planar graphs[END_REF][START_REF] Giménez | Graph classes with given 3-connected components: asymptotic enumeration and random graphs[END_REF][START_REF] Janson | The birth of the giant component[END_REF], for satisfiability problems [START_REF] Daudé | Random 2-XORSAT phase transition[END_REF], and for many other problems [START_REF] Flajolet | Analytic Combinatorics[END_REF][START_REF] Fuchs | Subtree sizes in recursive trees and binary search trees: Berry-Esseen bounds and Poisson approximations[END_REF][START_REF] Fuchs | Limit theorems for subtree size profiles of increasing trees[END_REF][START_REF] Hwang | Asymptotics of Poisson approximation to random discrete distributions: an analytic approach[END_REF][START_REF] Hwang | Phase changes in random recursive structures and algorithms[END_REF][START_REF] Hwang | Uniform asymptotics of Poisson approximation to the Poisson-binomial distribution[END_REF].

Usually, the main problems are to locate the phase transition, to properly describe the phase transition via special functions in terms of the involved parameters, and to give intuitive explanations of the observed phenomena. Some cases even exhibit two successive phase transitions. In probability theory, such a double phase transition occurs for example with the binomial distribution X n d = B(n, p) when p depends on n, with p bounded away from one. It indeed leads to the following trichotomy involving a continuous to discrete to degenerate phase transition: First, when both E(X n ) = p(n) • n and the variance V(X n ) tend to infinity, then a central limit theorem for

X n follows. Second, if p(n) • n → λ > 0, then X n is asymptotically Poisson distributed. Third, if p(n) • n → 0, then X n
degenerates to a Dirac distribution with all mass at 0.

The analysis of the size-refined composition scheme unravels, though more subtly, such a double phase transition: the random variable X n,j given by Equation (5) has three phases, each leading to its own limit law, visualized in Table 2.

Scale j ≪ n λ H 1+λ H j = Θ n λ H 1+λ H j ≫ n λ H 1+λ H Limit law continuous discrete discrete (gen. Mittag-Leffler ML) (mixed Poisson MPo(ξ ML)) (Dirac) Example X n,j ∼ C h j ρ j H n λ H ML X n,j ∼ MPo(ξ ML) P{X n,j ≥ 1} ∼ 0 TABLE 2
Three consecutive régimes for the number X n,j of H-components of size j in F -structures of size n in the critical size-refined scheme (4), depending on the relation between j and n;

this double phase transition is proven in Theorem 5.1.

Table 2 also motivates the following important remark.

REMARK 2.1 (Ubiquity of the exponent 1/3). Generating functions often have a dominant singularity of the square-root type (i.e., λ H = 1/2). This phenomenon is explained by the Drmota-Lalley-Woods theorem: Whenever H(z) can be defined by a strongly connected set of polynomial equations with nonnegative coefficients, it has a singular exponent 1/2 (see, e.g., [START_REF] Banderier | Formulae and asymptotics for coefficients of algebraic functions[END_REF][START_REF] Drmota | Random Trees[END_REF][START_REF] Flajolet | Analytic Combinatorics[END_REF]). Accordingly, in conjunction with Table 2, this explains why one observes a threshold at j = n 1/3 in many phase transitions; see Section 6.

One pleasant consequence of our work is that it gives a unified explanation of phase transitions from continuous to discrete observed in many examples: descendants in increasing trees [START_REF] Kuba | Descendants in increasing trees[END_REF], node degrees in increasing trees [START_REF] Kuba | On the degree distribution of the nodes in increasing trees[END_REF], block sizes in k-Stirling permutations [START_REF] Kuba | Analysis of statistics for generalized Stirling permutations[END_REF], stopping times in urn models [START_REF] Kuba | Limiting distributions for a class of diminishing urn models[END_REF], death processes [START_REF] Kuba | On death processes and urn models[END_REF], inversions in labelled tree families [START_REF] Panholzer | Limiting distributions for the number of inversions in labelled tree families[END_REF], ancestors and descendants in evolving k-tree models [START_REF] Panholzer | Ancestors and descendants in evolving k-tree models[END_REF].

These case by case studies lacked a proper comprehensive and uniform description of the arising phase transitions. So, instead of treating these combinatorial structures individually, we directly study the size-refined composition scheme [START_REF] Banderier | Formulae and asymptotics for coefficients of algebraic functions[END_REF]. As summarized in Table 2, we show how the phase transitions for the random variable X n,j depend on the growth of j = j(n) with respect to the size n. We prove that the distribution of X n,j(n) is continuous for small values of j (a three-parameter generalization of the Mittag-Leffler distribution), or discrete for some threshold values of j (a Poisson distribution mixed with the previous Mittag-Leffler distribution), or a Dirac distribution for large values of j. We further exemplify these results on different processes, like the Chinese restaurant process, sign changes and returns to zero in random walks, and the branching structure of random trees.

2.3.

Plan of the paper. In Section 3 we collect results from analytic combinatorics. We also present our basic assumptions on the generalized composition scheme and collect properties of various distributions that appear later in our main results. Section 4 is devoted to our results on the random variable X n corresponding to the extended composition scheme (2) involving F (z, u). Section 5 contains the results for the random variable X n,j corresponding to the size-refined composition scheme (4) involving F (z, v) and exhibiting phase transitions. We also give the covariance and the correlation coefficient of X n,j1 , X n,j2 , observing again some phase transitions. In Section 6 we discuss various examples to which we apply our results. Finally, in Section 7, we analyse further extensions for a cycle scheme and for a multivariate critical composition scheme. We also present new examples for these two extensions.

3. Singularity analysis, stable laws, and mixed Poisson distributions. In this section, we first present a few important notions from analytic combinatorics [START_REF] Flajolet | Analytic Combinatorics[END_REF] which we use to identify the radius of convergence and the singular exponents in our composition schemes. Then, we present a few results on the family of moment-tilted stable laws, based on James [START_REF] Fitzgerald | Lamperti-type laws[END_REF][START_REF] Fitzgerald | Stick-breaking P G(α, ζ)-generalized gamma processes[END_REF][START_REF] Fitzgerald | Generalized Mittag Leffler distributions arising as limits in preferential attachment models[END_REF] and Janson [START_REF] Janson | Moments of Gamma type and the Brownian supremum process area[END_REF]. We also collect properties of mixed Poisson distributions and their factorial moments [START_REF] Grandell | Mixed Poisson Processes[END_REF][START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF]. All of this allows us to identify in Sections 4 and 5 the distribution of the H-components in our composition schemes.

Real-Tauberian

ρ

Darboux-Pólya Singularity analysis

FIG 1.
As visually summarized by Flajolet in [START_REF] Flajolet | Some exactly solvable models of urn process theory[END_REF], three fundamental methods of asymptotic analysis require information on the function in different parts (shown here in red) of the complex plane. Flajolet and Odlyzko's singularity analysis [START_REF] Flajolet | Singularity analysis of generating functions[END_REF] offers more powerful results, but requires analyticity in a ∆-domain (tastefully also sometimes called "camembert domain" or "Pac-Man" by Flajolet himself!). This is the domain inside the blue curve, defined by ∆ = {z ∈ C such that |z| < ρ + ε and arg(z -ρ) > θ}, for some ε > 0 and 0 < θ < π/2. This analyticity is agreeably typically granted for most combinatorial constructions (e.g., for the ones leading to meromorphic, algebraic-logarithmic, hypergeometric, or D-finite functions).

3.1. Singularity analysis and asymptotic expansions. Let F (z) = n≥0 f n z n be a function with nonnegative coefficients f n that is analytic in a ∆-domain (see Figure 1 for this notion) with a finite radius of convergence ρ and singular expansion

F (z) = P 1 - z ρ + c F • 1 - z ρ λF (1 + o(1)) , (6)
where λ F ∈ R\{0, 1, 2, . . . } is called the singular exponent (of F (z) at z = ρ), and where P (x) ∈ C[x] is a polynomial (of degree ≥ 1 for λ F > 1, of degree 0 for 0 < λ F < 1, and P = 0 for λ F < 0). Then, by standard singularity analysis [START_REF] Flajolet | Singularity analysis of generating functions[END_REF], if ρ is the unique singularity of F (z) in |z| ≤ ρ, the Taylor series coefficients of F (z) satisfy

[z n ]F (z) = f n = c F ρ n • n -λF -1 Γ(-λ F ) • (1 + o(1)) . (7)
As the f n 's are nonnegative, this implies the sign property

sgn(c F ) = sgn Γ(-λ F ) , (8)
i.e., due to the sign change of the gamma function at each negative integer we have c F < 0 for 0 < λ F < 1, c F > 0 for λ F < 0, and sgn(c

F ) = (-1) ⌈λF ⌉ for λ F > 1.
Note that it is in general easy to get more asymptotic terms in [START_REF] Banderier | Random maps, coalescing saddles, singularity analysis, and Airy phenomena[END_REF], and that singularity analysis directly translates them into more asymptotic terms in [START_REF] Banderier | Enumeration and asymptotics of restricted compositions having the same number of parts[END_REF]. What is more, if one has several dominant singularities, one just has to sum the contributions of the local expansions at each singularity to get the asymptotics of the coefficients (this standard process is well presented in [START_REF] Flajolet | Analytic Combinatorics[END_REF]Chapter IV.6]; see also the rotation law in [START_REF] Banderier | The kernel method for lattice paths below a rational slope[END_REF] for walks or trees, which involve multiple dominant singularities as soon as their offspring distribution has a periodic support).

Note that algebraic functions constitute one of the main sources of functions satisfying the conditions of the expansion (6); indeed, they admit a Puiseux expansion

F (z) = k≥k0 c k • (1 -z/ρ) k/r ,
for k 0 ∈ Z and an integer r ≥ 1. For example, in Section 6.1 we will encounter the Catalan generating function 1/2 -√ 1 -4z/2, for which one has ρ = 1/4, P (x) = 1/2, k 0 = 0, and r = 2. This Catalan example is pleasantly simple, and thus obviously not generic, as here the Puiseux expansion contains only two terms. In full generality it can involve an infinite number of terms whose sum is converging. Let us now present this general case. LEMMA 3.1 (Singular expansion). Let F (z) be a power series with nonnegative coefficients satisfying [START_REF] Banderier | Random maps, coalescing saddles, singularity analysis, and Airy phenomena[END_REF]. Then F (z) has the following singular expansion

F (z) =            c F 1 -z ρ λF (1 + o(1)) if λ F < 0, τ F + c F 1 -z ρ λF (1 + o(1)) if 0 < λ F < 1, τ F + ⌊λF ⌋ i=1 p i 1 -z ρ i + c F 1 -z ρ λF (1 + o(1)) if λ F > 1,
where

τ F = F (ρ) > 0 for λ F > 0 and p 1 = -ρF ′ (ρ) < 0 for λ F > 1.
PROOF. First, if λ F < 0, then the lowest order of the Puiseux expansion is λ F . Second, for λ F > 0, we rewrite (6) into

F (z) = k i=0 p i (1 -z/ρ) i + c F (1 -z/ρ) λF + . . . .
Then, we have p 0 = P (0) = F (ρ) which we define to be τ F . We get τ F > 0 as it is an infinite convergent sum of nonnegative not-all-zero terms. Next, observe that the nonnegative coefficients of F (z) imply that F ′ (ρ) > 0. Thus, taking the derivative in the expansion of F (z) we get lim z→ρ F ′ (z) = +∞ for 0 < λ F < 1 and p 1 = -ρF ′ (ρ) ̸ = 0 for λ F > 1.

We now consider the critical scheme F (z) = G(H(z))M (z) where we assume that each of the functions G(z), H(z), and M (z) has a finite radius of convergence with a unique dominant singularity (i.e., the one of smallest modulus). By Pringsheim's theorem [40, p. 240] applied to each of these functions, the nonnegativity of its coefficients implies that its dominant singularity lies on the positive real axis and corresponds therefore to its radius of convergence denoted by ρ G , ρ H , and ρ M , respectively.

Note that if M (z) has an infinite radius of convergence (denoted by ρ M = +∞) or if ρ M ̸ = ρ H , then the asymptotics are easily obtained via

[z n ]F (z) ∼ G(H(ρ M ))[z n ]M (z) if ρ M < ρ H , (9a)
M (ρ H )[z n ]G(H(z)) if ρ M > ρ H . ( 9b
)
Now, as in (6), we define for each function: • the singular exponents λ G , λ H , and λ M , • the constant terms τ G , τ H , and τ M , • and the singular coefficients c G , c H , and c M . Thus, thanks to Equations (9a) and (9b), we can now focus (without loss of generality, or rather "without loss of difficulty!") on the case ρ M = ρ H which is more involved as here G(z), H(z), and M (z) are all contributing to the asymptotics in a nontrivial way. Then, one gets different régimes (depending on λ H ) for the asymptotics of the coefficient of F (z). In this article we focus on the range 0 < λ H < 1, while we treat the other range λ H > 1 in a companion article [START_REF] Banderier | Phase transitions of composition schemes: stable laws and Boltzmann distributions[END_REF]. Note that, as our work extends to some interesting combinatorial cases where M (z) has a radius of convergence ρ M > ρ H , we will also encompass this case, for which it is then convenient to set λ M = +∞ (the singular exponent of M (z) at z = ρ H is infinite whenever M is analytical there). The case λ M = +∞ is thus archetypal of cases where M (z) only affects the asymptotics of f n by a multiplicative constant like in Equation (9b).

We can now express the singular exponent of F in terms of those of G/H/M . LEMMA 3.2. Let F (z) = G(H(z))M (z) be a critical composition scheme that is singular at ρ H . Then, the singular exponent λ H of H(z) satisfies λ H > 0. Moreover, for the range 0 < λ H < 1, the singular exponent λ F of F (z) satisfies

λ F = min(λ G λ H , λ H , λ M , λ G λ H + λ M ).
For λ H > 1, the singular exponent λ F of F (z) satisfies

λ F = min(λ G , λ H , λ M , λ G + λ M ).
PROOF. The claim λ H > 0 follows from H(ρ H ) = ρ G ∈ (0, ∞) as one would have H(ρ H ) = +∞ if H(z) had a pole (or any algebraic singularity of negative singular exponent) at z = ρ H . Now, we plug the singular expansions from Lemma 3.1 at z = ρ H for G(z), H(z), and M (z) into F (z) = G(H(z))M (z). When 0 < λ H < 1 we get the following expansions (in which we omit the terms not contributing to the first-order asymptotics):

F (z)=                                  c M c G -c H ρ G λ G 1 - z ρ H λ G λ H +λ M + . . . if λ G <0, λ M <0, (10a)
τ M c G -c H ρ G λ G 1 - z ρ H λ G λ H + . . . if λ G <0, λ M >0, (10b)
τ M c G -c H ρ G λ G 1 - z ρ H λ G λ H +c M τ G 1 - z ρ H λ M +. . . if 0<λ G <1, (10c)
G ′ (ρ G )τ M c H 1 - z ρ H λ H + c M τ G 1 - z ρ H λ M + . . . if λ G >1, (10d)
where, if

λ M = +∞, τ M = c M = M (ρ H ) and 1 -z ρH λM = 0.
When λ H > 1 the linear term of H(z) is nonzero and therefore λ G plays the rôle of λ G λ H in the above expansions. So, for λ H > 1, one gets (omitting the constants for simplicity):

F (z)=                                  C 1 1 - z ρ H λ G +λ M + . . . if λ G <0, λ M <0, C 2 1 - z ρ H λ G + . . . if λ G <0, λ M >0, C 3a 1 - z ρ H λ G + C 3b 1 - z ρ H λ M + . . . if 0<λ G <1, C 4a 1 - z ρ H λ G + C 4b 1 - z ρ H λ M + C 4c 1 - z ρ H λ H + . . . if λ G >1.
Finally, the singular exponent λ F is equal to the minimal exponent in each of these Puiseux expansions.

This lemma motivates the following definition.

pure degenerate λ G 0 1 min(λ G λ H , λ H ) -λ M confluent FIG 2.
The three different régimes (pure, confluent, degenerate) for extended or size-refined composition schemes:

The Puiseux expansions of G/H/M go into resonance (or not), thus leading to these three cases.

DEFINITION 3.3 (Pure/confluent/degenerate composition schemes). Consider an extended or size-refined composition scheme (1) or (4) with a unique dominant singularity ρ F = ρ H , and with 0 < λ H < 1. It is either analytically

• pure if λ G < 0 or 0 < λ G < 1 and λ M > min(λ G λ H , λ H ); • confluent if 0 < λ G < 1 and λ M = min(λ G λ H , λ H ); • degenerate if λ G > 1 or 0 < λ G < 1 and λ M < min(λ G λ H , λ H ).
REMARK 3.4 (Dominant asymptotics and degenerate/confluent/pure composition schemes). Consider an extended or size-refined composition scheme (1) or (4) with a unique dominant singularity ρ F = ρ H .

• In the subcases (10a), (10b), and (10c) with λ M > λ G λ H , one gets a pure composition scheme. This adjective stresses the fact that the singular exponent of F is expressible, as expected, in terms of λ G and λ H . • In the subcase (10c) with λ M = λ G λ H , one gets a confluent composition scheme. • In the subcases (10c) with λ M < λ G λ H and (10d) the first-order asymptotics are dictated by M (z) or H(z) only, leading to a degenerate composition scheme. This adjective stresses the fact that the singular exponent of F is not expressible in terms of λ G .

We characterize the distributions associated with critical composition schemes in the analytically pure/confluent/degenerate cases in Theorems 4. 1, 4.3, and 4.4. We now present some probabilistic results on the distributions which will appear in these theorems.

3.2. Probability distribution melting pot. First, we discuss properties of tilted probability distributions and study in particular positive stable distributions. Then, we collect properties of a family of discrete distributions called mixed Poisson distributions [START_REF] Grandell | Mixed Poisson Processes[END_REF][START_REF] Karlis | EM algorithm for mixed Poisson and other discrete distributions[END_REF][START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF][START_REF] Willmot | Mixed compound Poisson distributions[END_REF], and we will end this section with a brief introduction to Boltzmann distributions.

For a random variable X of density f (x), the tilt of f (x) by a nonnegative integrable function g(x) is the following density g(x) E(g(X))

• f (x).

An important class of tilted densities are the polynomially tilted densities, where one tilts by a polynomial g(x) = x c (with c being any real value such that E(X c ) is well defined). We then use the notation

tilt c (f (x)) = x c E(X c ) • f (x).
Such tilted densities occur in many places: in the degree distribution in preferential attachment trees [START_REF] Fitzgerald | Stick-breaking P G(α, ζ)-generalized gamma processes[END_REF][START_REF] Fitzgerald | Generalized Mittag Leffler distributions arising as limits in preferential attachment models[END_REF], in Lamberti-type laws [START_REF] Fitzgerald | Lamperti-type laws[END_REF], in triangular urn schemes [START_REF] Janson | Limit theorems for triangular urn schemes[END_REF][START_REF] Janson | Moments of Gamma type and the Brownian supremum process area[END_REF], in node-degrees in plane-oriented recursive trees [START_REF] Kuba | On the degree distribution of the nodes in increasing trees[END_REF], and in table sizes in the Chinese restaurant process [START_REF] Aldous | Exchangeability and related topics[END_REF][START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF][START_REF] Pitman | Exchangeable and partially exchangeable random partitions[END_REF][START_REF] Pitman | Combinatorial Stochastic Processes[END_REF]. Note that many classes of distributions like the beta distribution, generalized gamma distribution [START_REF] Banderier | Periodic Pólya urns, the density method and asymptotics of Young tableaux[END_REF], the F -distribution, the beta-prime distribution, and distributions with gamma-type moments [START_REF] Janson | Moments of Gamma type and the Brownian supremum process area[END_REF] are closed under the tilting operation.

The following lemma shows that the operator tilt c admits in fact several equivalent definitions using the density, the moments, or the Laplace transform (see also [START_REF] Janson | Moments of Gamma type and the Brownian supremum process area[END_REF]Remark 2.11]). LEMMA 3.5 (Polynomially tilted density functions and moment shifts). Consider a random variable X with moment sequence (µ s ) s≥0 and density f (x) with support [0, ∞). Now consider a random variable X c with c ∈ N, having a distribution uniquely determined by its moments. Then the following properties are equivalent:

1. Tilted density: X c is a random variable with density

f c (x) = x c µc • f (x). 2. Shifted moments: X c is a random variable with moments E(X s c ) = µs+c µc . 3. Differentiated moment generating function: X c is such that E(e tXc ) = 1 µ c d c dt c E(e tX ).
REMARK 3.6 (Tilt with c ∈ R). For the properties (1) and (2) of Lemma 3.5, it is possible to extend their equivalence to c ∈ R, assuming that the corresponding moments exist. More generally, the equivalence between properties (1) and (2) stays valid for any random variable with density f (x) such that only moments µ 1 , . . . , µ n exist up to a certain value n ≥ 1. REMARK 3.7 (Densities with support R). In Lemma 3.5, if c is even, then one can drop the restriction that the support of f (x) is in [0, ∞). REMARK 3.8 (The tilt operator for densities/moments/random variables). This lemma justifies a slight abuse of notation: Starting with the densities of X and X c linked by tilt c (f (x)) = f c (x), the operator tilt c is also used to denote the corresponding tilted random variable tilt c (X) := X c and the corresponding tilted moments tilt c (µ s ) := µs+c µc .

PROOF OF LEMMA 3.5. For (1) ⇒ (2), first observe that f c (x) is indeed a density: One has f c (x) ≥ 0 on [0, ∞) and

∞ 0 f c (x) dx = µc µc = 1. Then, one checks that E(X s c ) = ∞ 0 x s f c (x) dx = ∞ 0 x s+c f (x) µ c dx = µ s+c µ c .
The fact that X c is uniquely determined by its moments then implies (2) ⇒ (1).

For (2) ⇔ (3), observe that

d c dt c E(e tX ) = d c dt c s≥0 µ s t s s! = s≥c µ s t s-c (s -c)! = s≥0 µ s+c t s s! ,
and, on the other hand, (1) and ( 2) imply

E(e tXc ) = ∞ 0 f c (x)e tx dx = 1 µ c s≥0 t s s! ∞ 0 x s+c f (x) dx = 1 µ c s≥0 t s s! µ s+c ,
which proves the claim.

Let us now introduce positive stable laws (also called one-sided stable laws, as their density has support (0, +∞)). DEFINITION 3.9 (Positive stable laws and their negative powers). We say that a positive random variable S α follows a stable law of parameter α ∈ (0, 1) if its Laplace transform is E(e -tSα ) = e -t α (see [START_REF] Vasilevich | Chance and stability. Stable distributions and their applications[END_REF] for a general presentation involving skewness, scale, and location parameters; they are respectively always 0, 1, and 0 in our work). The density of S α is4

f Sα (x) = ∞ n=1 (-1) n n!Γ(-nα) x -nα-1 (12) = 1 π ∞ n=1 (-1) n+1 Γ(nα + 1) sin(πnα) n! x -nα-1 (13) = 1 π α 1 -α π 0 K(ϕ) x 1/(1-α) exp - K(ϕ) x α/(1-α) dϕ, (14)
where

K(ϕ) = sin(αϕ) sin(ϕ) 1/(1-α) sin((1 -α)ϕ) sin(αϕ) . (15)
Formula ( 12) was first obtained by Humbert [START_REF] Humbert | Nouvelles correspondances symboliques[END_REF], and then rigorously proven by Pollard [START_REF] Pollard | The representation of e -x λ as a Laplace integral[END_REF] (see also Feller [33, Chapter XVII.6, Lemma 1], with the parameter γ = -α therein). Formula ( 14) is due, up to a typo that we corrected here, to Ibragimov and Chernin [START_REF] Abdulovich | On the unimodality of stable laws[END_REF]. Now, let β > 0 and define S α,β = (S α ) -β . Since P{S α,β ≤ x} = 1 -P{S α < x -1/β }, we directly obtain from [START_REF] Banderier | Lattice paths with catastrophes[END_REF] the density of S α,β on its support (0, +∞):

f Sα,β (x) = 1 β ∞ n=1 (-1) n n!Γ(-nα) x nα/β-1 = x -1/β-1 β f Sα (x -1/β ). (16)
Its moments are given by (see, e.g., Janson's survey on moments of Gamma type [START_REF] Janson | Moments of Gamma type and the Brownian supremum process area[END_REF]):

E(S s α,β ) = Γ( sβ α + 1) Γ(sβ + 1)
, s > -α β .

We will encounter composition schemes leading to powers of stable laws, an important subcase of it being the Mittag-Leffler distribution. 5EXAMPLE 3.10 (Mittag-Leffler distribution). We say that a random variable M α follows a Mittag-Leffler distribution ML(α

) if M α d = S α,α . Its moment generating function E(e xMα ) is the Mittag-Leffler function E α (x) = k≥0 x k Γ(1+αk) . An important special case is M 1 2 , the half-normal distribution |N (0, σ 2 )| with σ = √ 2; see Example 3.17 hereafter. For c > -α β consider the moment-tilted random variable X c = tilt c (S α,β ). Then, E(X s c ) = Γ( (s+c)β α + 1) Γ((s + c)β + 1) Γ(βc + 1) Γ( cβ α + 1) = Γ( (s+c)β α ) Γ((s + c)β) Γ(βc) Γ( cβ α ) ; (17)
see James [START_REF] Fitzgerald | Lamperti-type laws[END_REF][START_REF] Fitzgerald | Stick-breaking P G(α, ζ)-generalized gamma processes[END_REF][START_REF] Fitzgerald | Generalized Mittag Leffler distributions arising as limits in preferential attachment models[END_REF]. By ( 16) and Lemma 3.5, the density of X c is given by

f Xc (x) = Γ(βc + 1) βΓ( cβ α + 1) ∞ n=1 (-1) n n!Γ(-nα) x nα/β+c-1 . E.g., tilt 1 (M 1 2
) is the Rayleigh distribution of parameter √ 2; see Example 3.18 hereafter. ■

Both of them are defined in terms of the function Eα(x) introduced in 1903 by Mittag-Leffler [START_REF] Mittag-Leffler | Une généralisation de l'intégrale de Laplace-Abel[END_REF][START_REF] Mittag-Leffler | Sur la nouvelle fonction E α (x)[END_REF]. The first distribution (which we use in this article) was popularized by Feller [START_REF] Feller | Fluctuation theory of recurrent events[END_REF] (with a slight change of variable) and by Darling and Kac [START_REF] Darling | On occupation times for Markoff processes[END_REF] for the study of the local time of Markov processes. It has an exponentially bounded tail. The second one, which has a heavy tail, was introduced by Pillai [START_REF] Raman | On Mittag-Leffler functions and related distributions[END_REF].

Some probabilistic processes are related to a two-parameter generalization of this Mittag-Leffler distribution (see, e.g., [START_REF] Goldschmidt | Stable graphs: distributions and line-breaking construction[END_REF][START_REF] Fitzgerald | Generalized Mittag Leffler distributions arising as limits in preferential attachment models[END_REF][START_REF] Pitman | Combinatorial Stochastic Processes[END_REF]). We will establish in the next section that some critical composition schemes lead to this generalized Mittag-Leffler distribution. DEFINITION 3.11 (Generalized Mittag-Leffler distribution). For α ∈ (0, 1) and θ > -α, we say that a random variable X follows a generalized Mittag-Leffler distribution ML(α, θ) if X d = (tilt -θ (S α )) -α ; see James [START_REF] Fitzgerald | Generalized Mittag Leffler distributions arising as limits in preferential attachment models[END_REF]. This distribution is uniquely defined by its moments

E(X s ) = Γ s + θ α + 1 Γ(θ + 1) Γ(αs + θ + 1)Γ θ α + 1 = Γ s + θ α Γ(θ) Γ(αs + θ)Γ θ α . ( 18
)
Comparing these moments with [START_REF] Bernstein | Some canonical sequences of integers[END_REF] we directly get that for β = α and c = θ/α we have

tilt θ/α (ML(α)) d = ML(α, θ), i.e., tilt θ/α (S -α α ) = (tilt -θ (S α )) -α . ( 19
)
In other words, the permutation of the tilt and the power creates a change of the tilt parameter.

Next, we discuss product distributions. First, we recall properties of the beta distribution. DEFINITION 3.12 (Beta distribution). A beta-distributed random variable B d = Beta(α, β) with parameters α, β > 0 has a probability density function defined on (0, 1) by

f (x) = Γ(α + β) Γ(α)Γ(β) x α-1 (1 -x) β-1 .
The moments of B are given by

E(B s ) = Γ(s + α)Γ(α + β) Γ(s + α + β)Γ(α) , s > 0, (20)
and the beta distribution is uniquely determined by the sequence of its moments. Furthermore, let the reader be convinced of the convenient convention Beta(α, 0)

d = 1.
We now have all the ingredients to present the main properties of the distribution which will play a key rôle in the next sections, namely, the beta-Mittag-Leffler distribution. DEFINITION 3.13 (Beta-Mittag-Leffler distribution). We define the beta-Mittag-Leffler distribution BML(α, θ, β) as the distribution of the following product of independent random variables

Z d = Y • B α where Y d = tilt θ/α (S α,α ) d = ML(α, θ) and B d = Beta(θ, β)
are respectively distributed like a Mittag-Leffler and a beta distribution, with 0 < α < 1, θ > 0, and β ≥ 0.

LEMMA 3.14. The beta-Mittag-Leffler distribution BML(α, θ, β) has the following moments of order s (for s > 0):

E(Z s ) = Γ s + θ α Γ (θ + β) Γ (αs + θ + β) Γ θ α . ( 21
)
One has the following identity

Z d = ML(α, θ) Beta(θ, β) α d = ML(α, θ + β) Beta( θ α , β α ). (22)
PROOF. Due to the independence of the random variables we have E(Z s ) = E(Y s ) • E(B αs ). Then, using ( 17) and ( 20), one gets [START_REF] Carleman | Sur les équations intégrales singulières à noyau réel et symétrique[END_REF]. With the relation [START_REF] Bodini | Boltzmann samplers for first-order differential specifications[END_REF], this gives [START_REF] Darling | On occupation times for Markoff processes[END_REF].

One thus has the relations BML(α, θ, 0) = ML(α, θ), BML(α, α, 1 -α) = ML(α), and more generally BML(α, αk, 1 -α) = ML(α, α(k -1)).

Note that the beta-Mittag-Leffler distribution is one important instance of a "distribution with moments of Gamma type", a class of distributions popularized by Janson in his nice thorough survey [START_REF] Janson | Moments of Gamma type and the Brownian supremum process area[END_REF]. Therein, amongst several properties of these distributions, he obtains the asymptotics of their tail; applied to BML(α, θ, β), this implies that its density f (x) satisfies

f (x) ∼ Cx d-1 exp(-cx 1 1-α ) for x ∼ +∞, with c = (1 -α)α α 1-α , d = θ/α -θ -β + 1/2 1 -α , and C = Γ(θ + β) Γ( θ α ) α 1-2β 2(1-α) 2π(1 -α) .
Let us end this short introduction on the beta-Mittag-Leffler distribution by mentioning that an explicit series representation of its density is given later (in Formula (25) in the next section); we establish it via an inverse Mellin transform.

Another important ingredient of this article will be the mixed Poisson distributions. These distributions were first introduced by Dubourdieu in 1938 for actuarial mathematics/insurance modelling [START_REF] Dubourdieu | Les principes fondamentaux du calcul des probabilités et la théorie de l'assurance-maladie[END_REF], and then also studied by Lundberg and others (sometimes under the name "compound Poisson processes", a term that has a different meaning nowadays); they were also used by Neyman for applications in bacteriology [START_REF] Neyman | On a new class of "contagious" distributions applicable in entomology and bacteriology[END_REF], or for the analysis of some point processes in [START_REF] Grandell | Mixed Poisson Processes[END_REF]. Their unimodality properties are studied in [START_REF] Massé | Neyman type A distribution revisited[END_REF], and their tail asymptotics are analysed in [START_REF] Willmot | Lundberg Approximations for Compound Distributions with Insurance Applications[END_REF]. DEFINITION 3.15 (Mixed Poisson distributions). Let X denote a nonnegative random variable with cumulative distribution function U . We say that the discrete random variable Y has a mixed Poisson distribution with mixing distribution U and scale parameter ξ ≥ 0, if its probability mass function is given for ℓ ≥ 0 by

P{Y = ℓ} = ξ ℓ ℓ! R + X ℓ e -ξX dU = ξ ℓ ℓ! E(X ℓ e -ξX ).

This is summarized by the notation

Y d = MPo(ξU ), or, indifferently, Y d = MPo(ξX).
Note that mixed Poisson distributions provide a common generalization of three major discrete laws ubiquitous in combinatorics (see [START_REF] Flajolet | Analytic Combinatorics[END_REF]Figure IX.5]), namely the Poisson, the geometric, and the negative binomial distributions, making them of great importance per se.

We emphasize that the factorial moments 6 of a mixed Poisson distribution are closely related to the classical raw moments of its mixing distribution:

E(Y s ) = ξ s E(X s ), s ≥ 1.
Additionally, like for any distribution, the factorial and raw moments of Y are related via the Stirling set partition numbers s k (also called Stirling numbers of the second kind):

E(Y s ) = r k=0 s k E(Y k ).
We note in passing that such general relations are called Stirling transforms [START_REF] Bernstein | Some canonical sequences of integers[END_REF]. We refer to [START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF] for more properties of the Stirling transformation and mixed Poisson distributions; therein, Panholzer and the second author also gave the following useful expression for the probability mass function of Y d = MPo(ξX) in terms of its factorial moments. 6 Throughout this work we denote by x n the n th falling factorial,

x n = x(x -1) • • • (x -n + 1), n ≥ 0, with x 0 = 1.
It will be used for E(X n ), the factorial moment of order n of a random variable X. PROPOSITION 3.16. Let X denote a random variable with moment sequence given by (µ s ) s∈N such that E(e zX ) exists in a neighbourhood of zero, including the value z = -ξ. If a random variable Y has factorial moments given by E(Y s ) = ξ s µ s , then

Y d = MPo(ξX).
What is more, the sequence of moments of Y is the Stirling transform of the moment sequence (µ s ) s∈N , and the probability mass function of Y is given by

P{Y = ℓ} = s≥ℓ (-1) s-ℓ s ℓ µ s ξ s s! , ℓ ≥ 0.
Let us give two short examples of mixed Poisson distributions, which, as we shall later see, correspond to ubiquitous cases in combinatorics.

f (x; σ) = √ 2 σ √ π e -x 2 2σ 2 , x > 0;
alternatively, it is fully characterized by its moment sequence

E(X s ) = σ s 2 s/2 Γ( s+1 2 ) Γ( 1 2 )
.

Thus, a discrete random variable Y with probability mass function

P{Y = ℓ} = ξ ℓ ℓ! • √ 2 σ √ π ∞ 0 x ℓ e -ξx-x 2 2σ 2 dx, ℓ ≥ 0,
f (x; σ) = x σ 2 e -x 2 2σ 2 , x ≥ 0;
alternatively, it is fully characterized by its moment sequence

E(X s ) = σ s 2 s/2 Γ s 2 + 1 .
Thus, a discrete random variable Y with probability mass function

P{Y = ℓ} = ξ ℓ ℓ!σ 2 ∞ 0 x ℓ+1 e -ξx-x 2 2σ 2 dx, ℓ ≥ 0, has a mixed Poisson distribution: Y d = MPo(ξX) with X d = Rayleigh(σ).
Another representation valid for all ξ > 0 can be stated in terms of the incomplete gamma function

Γ(s, x) = ∞ x t s-1 e -t dt: P{Y = ℓ} = (ξσ) ℓ ℓ! e (ξσ) 2 2 ℓ+1 i=0 ℓ + 1 i (-ξσ) ℓ+1-i 2 i-1 2 Γ i + 1 2 , (ξσ) 2 2 . ■
In our results, we shall also encounter another important family of discrete distributions: the Boltzmann distributions.

Boltzman distributions were introduced in combinatorics by Duchon, Flajolet, Louchard, and Schaeffer in order to perform sampling of combinatorial structures [START_REF] Duchon | Boltzmann samplers for the random generation of combinatorial structures[END_REF]. The starting point of these authors was the idea to give, like in statistical mechanics, a Gibbs measure/Boltzmann weight x n (for some fixed real number x) to each combinatorial object of size n. Note that objects of the same size then follow a uniform distribution. It was then a nice surprise that, if one deals with an assemblage of combinatorial objects, the corresponding Boltzmann weights are given by very simple probabilistic laws (similarly to the symbolic method [START_REF] Flajolet | Analytic Combinatorics[END_REF] which directly gives the generating functions of unions/products/cycles of objects). This led to an outstanding generic linear time sampling algorithm: Its astonishing efficiency is partially due to the fact that, thanks to these Boltzmann weights, the sampling of a product of two combinatorial structures is simply obtained by two independent recursive subsamplings. The sampling algorithm is thus essentially based on the following definition: DEFINITION 3.19 (Boltzmann distribution). For any generating function G(z) = n≥0 g n z n , and for any parameter x > 0 inside the radius of convergence of G, a random variable X follows a Boltzmann distribution (associated with G) of parameter x, denoted by B G (x), if

P{X = n} = g n x n G(x) , n ≥ 0.
Then, the key idea behind Boltzmann sampling (of objects of size n) is to choose x adequately to maximize P{X = n}. If the object generated is not of size n, one rejects it and restarts the sampling. This leads to a uniform sampling algorithm of optimal efficiency when x is the unique real root of the equation xG ′ (x) = nG(x). This equation is reminiscent of many probabilistic results with mean µ = xG ′ (x)/G(x) (e.g., when G encodes the offspring of a Galton-Watson process). This is no coincidence: By design, Boltzmann sampling "reverseengineers" these results [START_REF] Panagiotou | Exact-size sampling of combinatorial structures in linear time[END_REF][START_REF] Sportiello | Boltzmann sampling of irreducible context-free structures in linear time[END_REF].

As we shall see, these Boltzmann distributions also occur in our critical composition schemes. Retrospectively, it explains and puts in a unified framework earlier sporadic occurrences of such distributions for the limit law of the degree of a random node in simplygenerated trees, the root degree in simply-generated trees, as well as in subcritical composition schemes; see [40, pages 460, 629-633]. 4. Extended composition scheme. In the following we state and prove our main theorem on the extended critical composition scheme (2) for pure schemes. For the terms critical and pure, we refer to Definitions 1.1 and 3.3, respectively. This theorem shows the universality of the beta-Mittag-Leffler distribution (introduced in Definition 3.13). THEOREM 4.1 (Extended composition scheme: pure case). In a pure extended critical composition scheme F (z, u) = G uH(z) M (z), the core size X n , rescaled, converges in distribution and in moments 9 to a random variable X distributed like a beta-Mittag-Leffler distribution:

λ G = -10 λ G = -2 λ G = -1 λ G = 1 3 FIG 3. Boltzmann distributions B G (x)
X n κ • n λH d --→ m X, with X d = BML(α, θ, β), (23)
where

α = λ H , θ = -λ G λ H , β = -λ - M = -min(0, λ M ), and κ = τ H -c H .
What is more, one has a local limit theorem

P{X n = x • κn λH } ∼ 1 κn λH • f X (x), (24)
where f X (x) is the density of X:

f X (x) = Γ(θ + β) Γ(θ/α) j≥0 (-1) j j!Γ(β -jα) x θ/α+j-1 . ( 25
)
REMARK 4.2 (Two simplifications of the beta-Mittag-Leffler distribution). In the above theorem, if λ M ≥ 0 (which includes the critical scheme F (z, u) = G(uH(z)) as λ M = +∞), then β = 0 and the beta-Mittag-Leffler distribution [START_REF] Daudé | Random 2-XORSAT phase transition[END_REF] simplifies into a generalized Mittag-Leffler distribution:

X d = tilt -λG (ML(λ H )) d = ML(λ H , -λ G λ H ).
In particular, for λ H = 1 2 and λ G = -1 the random variable X follows a Rayleigh distribution of parameter σ = √ 2; see Example 3.10. Another noteworthy simplification occurs for λ M < 0, in the special case λ G = -1 and λ H -λ M = 1: We then obtain a Mittag-Leffler distribution of parameter λ H :

X d = ML(λ H ).
In particular, for λ H = 1 2 the random variable X follows a half-normal distribution of parameter σ = √ 2; see again Example 3.10. Note that the cases with λ G = -1 occur in a great many places in applied probability theory; they indeed correspond to a natural combinatorial framework where F -objects are essentially sequences of H-components (see [START_REF] Flajolet | Analytic Combinatorics[END_REF] and examples in our Section 6). PROOF OF THEOREM 4.1. The factorial moments satisfy 10 :

E(X s n ) = [z n ]∂ s u (F )(z, 1) [z n ]F (z, 1) = [z n ]H(z) s G (s) H(z) M (z) [z n ]G(H(z))M (z)
. 9 We write Xn d --→ m X to denote that Xn converges in distribution to X, with convergence of all moments, i.e., E(X s n ) → E(X s ) for all s. This notion of convergence in moments was, e.g., used in [START_REF] Fuchs | Profiles of random trees: Limit theorems for random recursive trees and binary search trees[END_REF][START_REF] Janson | Moment convergence of balanced Pólya processes[END_REF]. It is also indirectly used in [START_REF] Rösler | The contraction method for recursive algorithms[END_REF], which deals with more constrained models that offer a convergence in L p , for all p > 1. Note that in all our results involving convergence in distribution or in moments, we omit the speed of convergence, which is in fact easily obtained by considering the Puiseux expansions of order 2 of G/H/M . 10 Throughout this work, we denote by ∂u the differentiation operator with respect to the variable u. Accordingly, we use the shorthand notation ∂u(F )(z, 1) = ∂uF (z, u) u=1 .

In the following we use the notation ρ F , τ F , λ F , and c F from Section 3.1 for the singular expansions of F = G/H/M . As we are in a pure critical scheme (Definition 3.3), the unique singularity of F (z) = G(H(z))M (z) is at z = ρ H . Then, we unify the three cases (10a), (10b), and (10c) by using λ - M = min(0, λ M ) and choosing C M according to the specific case (C M is for λ M ̸ = λ G λ H either τ M or c M ; this quantity C M will anyway cancel in the end). Note that the case (10d) does not hold in a pure scheme. This gives

F (z) ∼ C M c G -c H ρ G λG 1 - z ρ H λGλH +λ - M .
Therefore, using singularity analysis we get

f n = [z n ]F (z, 1) ∼ C M c G (-c H /ρ G ) λG ρ n H • n -λGλH -λ - M -1 Γ(-λ G λ H -λ - M ) . ( 26
)
Using singular differentiation [START_REF] Allen Fill | Singularity analysis, Hadamard products, and tree recurrences[END_REF][START_REF] Flajolet | Analytic Combinatorics[END_REF] for G(z), we get the following singular expansion of the higher-order derivatives G (s) (z), for integer s ≥ 1:

G (s) (z) ∼ (-1) s c G ρ s G λ s G 1 - z ρ G λG-s
.

From this we get

G (s) (H(z)) ∼ (-1) s c G ρ -λG G λ s G (-c H ) λG-s 1 - z ρ H λGλH -sλH .
Next, from the singular expansion of H(z) we directly get

H(z) s ∼ τ s H -sτ s-1 H • (-c H ) 1 - z ρ H λH .
Combining these expansions with the one of M (z) gives the required expansion

H(z) s G (s) H(z) M (z) ∼ (-1) s τ s H C M c G ρ -λG G λ s G (-c H ) λG-s 1 - z ρ H λ - M +λH λG-sλH
.

Next we rewrite (-1) s λ s G using the gamma function:

(-1) s λ s G = Γ(s -λ G ) Γ(-λ G ) .
Hence, we obtain by extraction of coefficients and singularity analysis

[z n ]H(z) s G (s) H(z) M (z) ∼ τ H -c H s C M c G (-c H /ρ G ) λG ρ n H × Γ(s -λ G ) Γ(-λ G ) • n -λGλH -λ - M -1+sλH Γ(sλ H -λ G λ H -λ - M ) . (27)
Combining this expression with [START_REF] Dolera | A Berry-Esseen theorem for Pitman's α-diversity[END_REF] gives

E(X s n ) ∼ n sλH κ s • Γ(s -λ G )Γ(-λ G λ H -λ - M ) Γ(-λ G )Γ(sλ H -λ G λ H -λ - M )
,

with κ = τH -cH . What is more, one has E(X s n ) ∼ E(X s n
) since we can express the raw moments by using the factorial moments and the Stirling numbers of the second kind:

E(X s n ) = E s j=0 s j X j n = s j=0 s j E(X j n ). (28)
Consequently, we obtain the moment convergence to the stated moment sequence:

E(X s n ) n sλH κ s → µ s .
By Stirling's formula for the gamma function, one has

Γ(z) = z e z √ 2π √ z 1 + O 1 z ; (29)
this entails that for s → ∞ the moments satisfy

(µ s ) -1/(2s) ∼ e 1-λH λ λH H s H-1 2 1 + O 1 s .
As 0 < λ H < 1, the divergence in Carleman's criterion [21, pp. 189-220] is satisfied:

∞ s=0 µ -1/(2s) s = +∞; (30)
consequently, the moment sequence (µ s ) s∈N characterizes a unique distribution. Thus, by the Fréchet-Shohat theorem [START_REF] Fréchet | A proof of the generalized second limit theorem in the theory of probability[END_REF], we obtain the weak convergence of the normalized random variable Xn κ•n λ H to a random variable X with moment sequence (µ s ) s∈N . Concerning the local limit theorem, we have to analyse

P{X n = k} = gk[z n ]H(z) k M (z) fn , for k = x • κ • n λH ,
with x in a compact subinterval of (0, ∞). First, we note that by [START_REF] Banderier | Enumeration and asymptotics of restricted compositions having the same number of parts[END_REF] applied to G(z) we directly obtain

g k ∼ c G ρ k G • (κx) -λG-1 n -λH λG-λH Γ(-λ G ) . (31)
It remains to determine the asymptotics of [z n ]H(z) k M (z). We have

[z n ]H(z) k M (z) = 1 2πi H(z) k M (z) z n+1 dz.
Introducing the point A of coordinates ( 1 n , ρ H (1 + log 2 n n )), this Cauchy integral can be transformed into an integral over a larger contour in the Delta-domain (in blue in Figures 1 and4). Then, setting z = ρ H (1 + t/n) leads to an integral asymptotically concentrated on the Hankel contour C (which starts and ends at +∞): Then, expanding e -x(-t) λ H leads to

[z n ]H(z) k M (z) ∼ τ k H C M ρ n H n 1+λ - M • 1 2πi C (-t) λ - M e -t-x(-t) λ H dt.
[z n ]H(z) k M (z) ∼ τ k H C M ρ n H n 1+λ - M • 1 2πi C e -t j≥0 (-x) j j! (-t) jλH +λ - M dt, (32)
in which we recognize the following Hankel contour representation [START_REF] Taylor | A Course of Modern Analysis[END_REF]Section 12.22]:

1 Γ(-z) = 1 2iπ C (-t) z e -t dt.
Recall that we use, by analytic continuation, 1/Γ(m) = 0 for any integer m < 0; this allows us to avoid more cumbersome expressions which make use of Euler's reflection formula

1 Γ(-z) = -1
π sin(πz)Γ(1 + z). Now, combining the expansions of f n from ( 26), g k from ( 31), and term-wise integration in (32), we get the desired local limit theorem:

P{X n = k} = g k f n [z n ]H(z) k M (z) ∼ Γ(-λ G λ H -λ - M ) κn λH Γ(-λ G ) j≥0 (-1) j j!Γ(-jλ H -λ - M ) x j-λG-1 . ( 33
)
It remains to verify that f X (x) is the density function of X with moments

µ s = Γ(s -λ G )Γ(-λ G λ H -λ - M ) Γ(-λ G )Γ(sλ H -λ G λ H -λ - M )
.

Note that the gamma function is never zero and that its only singularities are simple poles at the negative integers. Therefore µ s (considered as a complex function of s) has simple poles at ρ 0 = λ G (if λ - M ̸ = 0, which entails λ G < 0 as we have a pure scheme) and at ρ j = λ G -j (for j ∈ {1, 2, . . . }). Then, as µ s-1 is the Mellin transform of the density f X (x) of X, an inverse Mellin computation 11 implies that f X (x) is expressible in terms of the residues of µ s :

f X (x) = 1 2πi σ+i∞ σ-i∞ µ s-1 x -s ds = σj = 1 + ρj pole of µs-1 σj<σ
Res s=σj (µ s-1 )x -σj , which is valid for x > 0, and where σ = 1 + ρ 0 + ϵ is in the fundamental strip of µ s-1 , and

Res s=σj (µ s-1 ) = Res s=ρj (µ s ) = (-1) j j! • Γ(-λ G λ H -λ - M ) Γ(-λ G )Γ(-jλ H -λ - M )
.

Summing for j ≥ 0 (we can include ρ 0 , even if it is not a pole, as the corresponding residue is then 0) gives the same density function f X (x) as in the local limit theorem [START_REF] Feller | An Introduction to Probability Theory and its Applications[END_REF]. Now, thanks to our probability distribution melting pot section 3.2, we identify the corresponding limit law by matching the parameters in the moments of the beta-Mittag-Leffler product given in Lemma 3.14.

Above, we have established the limit laws occurring for analytically pure composition schemes. Next we deal with the analytically confluent and degenerate cases of Definition 3.3. They lead either to discrete distributions, or, more interestingly, to a mixture of discrete and continuous distributions. This generalizes the phenomenon observed in [START_REF] Banderier | Random maps, coalescing saddles, singularity analysis, and Airy phenomena[END_REF], where for λ H = 3 2 the limit law also consists of a discrete part plus a continuous part, a map-Airy distribution (this phenomenon also occurs for variants of 3-connected graphs; see [START_REF] Giménez | Graph classes with given 3-connected components: asymptotic enumeration and random graphs[END_REF]). The following two theorems explain the connection between this discrete part and Boltzmann distributions (usually used for random sampling; see, e.g., [START_REF] Bodini | Boltzmann samplers for first-order differential specifications[END_REF][START_REF] Duchon | Boltzmann samplers for the random generation of combinatorial structures[END_REF]). THEOREM 4.3 (Extended composition scheme: degenerate case). Let a degenerate extended critical composition scheme F (z, u) = G uH(z) M (z) with 0 < λ H < 1 be given.

For 0 < λ G < 1 and λ M < λ G λ H , the core size X n converges for k ≥ 0 and n → ∞ to a Boltzmann distribution B G (ρ G ) (see Definition 3.19):

P{X n = k} → P{B G (ρ G ) = k} = g k ρ k G G(ρ G ) .
For λ G > 1, the core size X n has for k ≥ 0 and n → ∞ the following behaviour:

• For λ M < λ H , the random variable X n converges to a Boltzmann distribution B G (ρ G ):

P{X n = k} → g k ρ k G G(ρ G ) .
• For λ M = λ H , the random variable X n converges to a convex combination of two discrete Boltzmann distributions:

X n d -→ Be(p) • B G (ρ G ) + (1 -Be(p)) • B G ′ (ρ G ), with p = cM G(ρG)
cM G(ρG)+cH τM G ′ (ρG) and where the Bernouilli distribution Be(p) is independent of the other random variables. Therefore, we have

P{X n = k} → p • g k ρ k G G(ρ G ) + (1 -p) • kg k ρ k-1 G G ′ (ρ G ) .
• For λ M > λ H , the random variable X n converges to a Boltzmann distribution B G ′ (ρ G ):

P{X n = k} → kg k ρ k-1 G G ′ (ρ G ) .
PROOF. Let us start with the case 0 < λ G < 1. We study for arbitrary but fixed k ∈ N the probability P{X n = k}, as n tends to infinity. From (3) we obtain

P{X n = k} = g k [z n ]H(z) k • M (z) f n ∼ g k f n [z n ] τ k H + kτ k-1 H c H 1 - z ρ H λH τ M + c M (1 - z ρ H ) λM , (34)
where we apply the Singular-expansion Lemma 3.1 to the function H(z). From the expansion (10c), we get by standard singularity analysis [START_REF] Banderier | Enumeration and asymptotics of restricted compositions having the same number of parts[END_REF], for λ M < λ G λ H , that f n satisfies

f n ∼ τ G c M 1 ρ n H n -λM -1 Γ(-λ M ) ,
as the singular exponent λ M dominates the asymptotics. Thus we obtain from ( 34)

P{X n = k} ∼ g k τ k H τ G = g k ρ k G G(ρ G ) , because τ H = ρ G and τ G = G(ρ G
) by our assumptions on criticality.

Let us continue with the case λ G > 1. From expansion (10d) we see that we now have to distinguish three subcases:

λ M < λ H , λ M = λ H , or λ H < λ M .
The case λ M < λ H is exactly the same as before. The case λ M > λ H is different, yet the results are derived analogously: the asymptotics of f n depend now on λ H and are given by

f n ∼ G ′ (ρ G )τ M c H 1 ρ n H n -λH -1 Γ(-λ H ) .
Then, we obtain from [START_REF] Ferrari | Multivariate mixed Poisson distributions[END_REF] where now again the contribution from H(z) dominates

P{X n = k} ∼ kg k τ k-1 H G ′ (ρ G ) = kg k ρ k-1 G G ′ (ρ G ) .
Finally, in the case λ M = λ H the previous two contributions are mixed as the first two terms in (10d) contribute: The asymptotics of f n are given by

f n ∼ G ′ (ρ G )τ M c H + c M τ G 1 ρ n H n -λH -1 Γ(-λ H ) .
Note that the coefficient is not zero, as both terms are negative: G ′ (ρ g ), τ M , τ M > 0 and c H , c M < 0 due to the sign property (8) for 0 < λ H < 1. Then, we obtain from [START_REF] Ferrari | Multivariate mixed Poisson distributions[END_REF] where again both contributions have to be taken into account

P{X n = k} ∼ c M g k ρ k G + τ M c H kg k ρ k-1 G G ′ (ρ G )τ M c H + c M τ G = p • g k ρ k G G(ρ G ) + (1 -p) • kg k ρ k-1 G G ′ (ρ G )
,

where p = cM G(ρG)
cM G(ρG)+cH τM G ′ (ρG) ∈ (0, 1) by the sign property [START_REF] Banderier | Analysis of three graph parameters for random trees[END_REF]. We thus get a linear combination of two Boltzmann distributions, weighted by a Bernoulli random variable Be(p). THEOREM 4.4 (Extended composition scheme: confluent case). Let a confluent extended critical composition scheme F (z, u) = G uH(z) M (z) with 0 < λ H < 1 be given (i.e., 0 < λ G < 1 and λ M = λ G λ H ). Then the core size X n is a convex combination of a Boltzmann distribution B G (ρ G ) and an asymptotically continuous random variable Z n :

X n ∼ Be(p) • B G (ρ G ) + (1 -Be(p)) • Z n , Z n κ • n λH d -→ X,
where κ and the limit law

X d = ML(λ H , -λ G λ H )
= cM G(ρG) cM G(ρG)+τM cG(-cH /ρG) λ G ) is independent of B G (ρ G ), Z n ,
and X. PROOF. The proof follows the same lines as the one of Theorem 4.3. We start from expansion (10c). Due to λ M = λ G λ H both terms contribute, and we get

f n ∼ τ M c G (-c H /ρ G ) λG + τ G c M ρ n H n -λM -1 Γ(-λ M ) . (35)
Then we extract the asymptotics from [START_REF] Ferrari | Multivariate mixed Poisson distributions[END_REF], where now the contributions from M (z) dominate, as λ M < λ G :

P{X n = k} = c M g k τ k H τ M c G (-c H /ρ G ) λG + τ G c M + o(1) = p • g k ρ k G G(ρ G ) + o(1), (36)
Using the sign property (8) we get p ∈ (0, 1). So, for large n, X n behaves with probability p like B G (ρ G ), but what happens with probability 1 -p? Where is this missing mass in [START_REF] Flajolet | Some exactly solvable models of urn process theory[END_REF]? For sure, it is sneakily spread in k o(1): It turns out that more and more mass is distributed in the range k ∼ Θ(κn λH ), leading to an asymptotic continuous distribution therein.

In order to identify this distribution, we compute the factorial moments of X n like in the proof of Theorem 4.1. We again use the singular expansions of G (s) H(z) , H(z) s , and 27) holds verbatim with C M = τ M . The big difference lies now in the asymptotics of f n : It is given by [START_REF] Dolera | A Berry-Esseen theorem for Pitman's α-diversity[END_REF] in the pure case and by [START_REF] Allen Fill | Singularity analysis, Hadamard products, and tree recurrences[END_REF] in the confluent case. Thus, after rescaling (27) by f n , the factorial moments have the same shape, but with an additional prefactor 1 -p: E(X

H(z) s G (s) H(z) M (z). Formula (
s n ) ∼ (1 -p) • E(Z s n
), which proves the claim. REMARK 4.5 (Physical interpretation of the bimodal behaviour). As illustrated by Table 1 on page 4, the confluent case gives a bimodal distribution. The first mode is dictated by small values of k, where a Boltzmann distribution is dominant, while for larger values of k ≈ n λH , a second mode appears, which is associated with a continuous Mittag-Leffler distribution. This phenomenon explains the following behaviour: If one picks an atom at random in a structure of size n, then with probability p it lies in a large H-block of size Θ(n) and with probability 1 -p in a smaller H-block of size Θ(n 1-λH ).

In the next section, we will refine these considerations by having a closer look at the distribution of the H-blocks of any given size.

5. Size-refined composition scheme. In this section, we give the limit laws for the profile of combinatorial structures given by a critical composition scheme. We focus here on schemes which are analytically pure (see Definition 3.3), while we handle the confluent and degenerate cases in our companion article [START_REF] Banderier | Phase transitions of composition schemes: stable laws and Boltzmann distributions[END_REF] (as they require additional technical details and different families of limit laws, which also pop up for λ H > 1). The profile is captured by the size-refined composition scheme (4). As we see in the theorem below, we get three distinct asymptotic régimes, each leading to its own limit law. Two of these limit laws are expressible in terms of the generalized Mittag-Leffler distribution of Theorem 4.1 (see also Definition 3.11 and Definition 3.13). THEOREM 5.1 (Mixed Poisson limit behaviour for the size-refined scheme). Consider a size-refined pure critical composition scheme

F (z, v) = G H(z) -(1 -v)h j z j M (z), with j ∈ N. Let ξ n,j = ρ j H -cH h j n λH ,
and X be the beta-Mittag-Leffler distribution of Theorem 4.1:

X d = BML(α, θ, β),
where α = λ H , θ = -λ G λ H , and β = -min(0, λ M ).

Then, the random variable X n,j , which counts the number of H-components of size j possesses three distinct asymptotic régimes, with a phase transition at j = Θ(n λ H 1+λ H ):

(i) For j ≪ n λ H
1+λ H , we have ξ n,j → +∞ and convergence in distribution and in moments:

X n,j ξ n,j d --→ m X. (ii) For j ∼ r • n λ H 1+λ H , r ∈ (0, ∞), we have ξ n,j → ξ with ξ = r - λ H 1+λ H • 1 -Γ(-λH
) and convergence in distribution and in moments towards a mixed Poisson distribution:

X n,j d --→ m MPo(ξX).

(iii) For j ≫ n λ H

1+λ H , we have ξ n,j → 0, and X n,j converges to a Dirac distribution at 0. REMARK 5.2 (Phase transition I). The intuition behind the phase transition is as follows:

In the limit n → ∞, there are many small (j ≪ n λ H 1+λ H ), some giant (j ∼ rn λ H 1+λ H ), and no super-giant (j ≫ n λ H 1+λ H ) H-components of size j. It is interesting to compare this situation with the birth of the giant component in Erdős-Rényi random graphs; see [START_REF] Janson | The birth of the giant component[END_REF]. Note that the case j ∈ N fixed (i.e., independent of n as n tends to infinity) falls into the case j ≪ n λ H 1+λ H . REMARK 5.3 (Phase transition II). For the often observed case of a square-root singularity of H(z) (i.e., λ H = 1 2 ), we reobtain the critical range j = Θ(n 1/3 ), which was already observed in the mixed Poisson Rayleigh distributions in [START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF]. Furthermore, as one has

E(X s n,j ) = ξ s n,j • E(X s ) • (1 + o(1)
), this offers en passant a link between the ξ n,j 's and the rescaling factor κn λH in Theorem 4.1:

j≥1 ξ n,j = j≥1 ρ j H -c H h j n λH = H(ρ H ) -c H n λH = τ H -c H n λH = κn λH .
This link can be seen as an asymptotic avatar of the combinatorial relation j≥1 X n,j = X n , implying

j≥1 E(X n,j ) = E(X n ).
In order to prove Theorem 5.1, we first need the following lemma concerning the convergence of mixed Poisson distributions. LEMMA 5.4 (Factorial moments and limit laws of mixed Poisson type [START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF]). Let (X n ) n∈N denote a sequence of random variables, whose factorial moments are asymptotically of mixed Poisson type, i.e., they satisfy for n → ∞ the asymptotic expansion

E(X s n ) = ξ s n • µ s • (1 + o(1)
), s ≥ 1, with µ s ≥ 0 and ξ n > 0. Furthermore, assume that the moment sequence (µ s ) s∈N determines a unique distribution X satisfying Carleman's condition. Then, the following limit distribution results hold: PROOF OF THEOREM 5.1. The factorial moments E(X s n,j ) = k≥0 P{X n,j = k}k s of X n,j are obtained from F (z, v) by repeated differentiation and evaluation at s = 1:

(i) If ξ n → ∞,
E(X s n,j ) = [z n ]∂ s v (F )(z, 1) [z n ]F (z, 1) = h s j [z n-j•s ]G (s) H(z) M (z) f n .
We already know the asymptotics of f n from [START_REF] Dolera | A Berry-Esseen theorem for Pitman's α-diversity[END_REF]. For fixed j we can proceed by extraction of coefficients, while for j = j(n) tending to infinity, the asymptotic expansion of h j follows by singularity analysis applied to H(z) (see Equation ( 7)):

h j = c H ρ j H • j -λH -1 Γ(-λ H ) • (1 + o(1)) . (37)
What is more, one has

G (s) H(z) M (z) ∼ (-1) s c M c G ρ -λG G λ s G (-c H ) λG-s 1 - z ρ H λH λG-sλH +λ - M . (38)
This implies that X n,j has factorial moments of mixed Poisson type:

E(X s n,j ) ∼ h s j (-c H ) -s ρ js H • µ s • n sλH with µ s := Γ(s -λ G )Γ(-λ G λ H -λ - M ) Γ(-λ G )Γ(-λ H λ G + sλ H -λ - M )
.

We already observed that the moment sequence (µ s ) s∈N determines a unique distribution, by Carleman's criterion [START_REF] Dubourdieu | Les principes fondamentaux du calcul des probabilités et la théorie de l'assurance-maladie[END_REF]. Thus, the limit laws follow by using Lemma 5.4. Finally, the critical growth range is obtained via the closed-form expression for ξ n,j (in which one inserts the expansion (37)): Indeed, ξ n,j is of growth order n λ H j 1+λ H Γ(-λH ) and converges to a nonzero constant if and only if j(n) ∼ r • n λ H 1+λ H , therefore the critical growth range is Θ(n λ H 1+λ H ). Collecting all contributions from ξ n,j (for j = o(n) tending to infinity) gives the constant ξ. The Dirac case is now finally obtained by an additional analysis of the expected value in the remaining range j ≫ n λ H 1+λ H . There, we directly obtain E(X n,j ) → 0, which proves the stated result.

Next we turn to the dependence between the number of H-components of size j 1 and j 2 , determining the covariance and the correlation coefficient. THEOREM 5.6. In a size-refined pure critical composition scheme, the covariance of the random variables X n,j1 and X n,j2 , counting the number of H-components of size j 1 and j 2 (with

j 1 , j 2 = o(n)), satisfies Cov(X n,j1 , X n,j2 ) ∼ ξ n,j1 ξ n,j2 • V(X), (39)
where ξ n,j = ρ j H -cH h j n λH > 0 and X d = BML(α, θ, β) denotes the beta-Mittag-Leffler distribution from Theorem 4.1. Furthermore, the correlation coefficient between X n,j1 and X n,j2 satisfies

ρ(X n,j1 , X n,j2 ) ∼ 1 1 + E(X) ξ1 1 1 + E(X) ξ2 ,
where, for k = 1, 2, one has ξ k := lim n ξ n,jk and

1 + E(X) ξk ∼ 1 if j k ≪ n λ H 1+λ H .
REMARK 5.7. We observe that for small j 1 , j 2 (e.g., if j 1 , j 2 = O(1)), the random variables are asymptotically highly correlated: ρ(X n,j1 , X n,j2 ) ∼ 1.

PROOF OF THEOREM 5.6. The combinatorial scheme

F = G H ̸ =j1,j2 +v 1 H =j1 +v 2 H =j2 × M (where one takes distinct sizes 1 ≤ j 1 < j 2 ) directly translates into F (z; v 1 , v 2 ) = G H(z) -(1 -v 1 )h j1 z j1 -(1 -v 2 )h j2 z j2 M (z)
. Accordingly, the random variables X n,j1 and X n,j2 have the joint distribution

P{X n,j1 = k 1 , X n,j2 = k 2 } = [z n v k1 1 v k2 2 ]F (z; v 1 , v 2 ) [z n ]F (z; 1, 1)
.

We already know the asymptotics of f n = [z n ]F (z; 1, 1), given in [START_REF] Dolera | A Berry-Esseen theorem for Pitman's α-diversity[END_REF]. We get by differentiation and evaluation at

v 1 = v 2 = 1 E(X n,j1 X n,j2 ) = [z n ]∂ v1 ∂ v2 (F )(z; 1, 1) [z n ]F (z; 1, 1) = h j1 h j2 [z n-j1-j2 ]G ′′ (H(z))M (z) f n .
The asymptotics of h j1 and h j2 are given in [START_REF] Flajolet | Some exactly solvable models of urn process theory[END_REF]. The singular expansion of G ′′ (H(z))M (z) is a special case of (38), so we obtain for j 1 , j 2 = o(n):

E(X n,j1 X n,j2 ) ∼ h j1 h j2 c -2 H ρ j1+j2 H • E(X 2 )
• n 2λH . Hence, using the explicit form of E(X) in [START_REF] Carleman | Sur les équations intégrales singulières à noyau réel et symétrique[END_REF], we obtain for the covariance: 2 and the definition of ξ n,j , this implies [START_REF] Flajolet | Singularity analysis of generating functions[END_REF]. For the correlation coefficient, we observe that Theorem 5.1 together with Lemma 5.4 implies V(X n,j ) ∼ ξ 2 n,j V(X) + ξ n,j E(X). Collecting all contributions, this implies that

Cov(X n,j1 , X n,j2 ) = E(X n,j1 X n,j2 ) -E(X n,j1 )E(X n,j2 ) ∼ h j1 h j2 ρ j1+j2 H c 2 H E(X 2 ) -E(X) 2 n 2λH . By V(X) = E(X 2 ) -E(X)
ρ(X n,j1 , X n,j2 ) = Cov(X n,j1 , X n,j2 ) V(X n,j1 ) V(X n,j2 ) ∼ 1 1 + E(X) ξn,j 1 1 + E(X)
ξn,j 2

.

Taking the limit of this expression for n going to infinity concludes the proof.

Before to present in Section 7 a multivariate generalization of these results, we now list several applications to a variety of combinatorial structures.

Applications and examples.

Our main results, Theorems 4.1 and 5.1, can be readily applied to the problems considered by Drmota and Soria [START_REF] Drmota | Images and preimages in random mappings[END_REF], Flajolet and Sedgewick [START_REF] Flajolet | Analytic Combinatorics[END_REF], Dumas, Flajolet and Puyhaubert [START_REF] Flajolet | Some exactly solvable models of urn process theory[END_REF], Janson [START_REF] Janson | Random cutting and records in deterministic and random trees[END_REF], Meir and Moon [START_REF] Meir | Cutting down random trees[END_REF], Panholzer and Seitz [START_REF] Panholzer | Limiting distributions for the number of inversions in labelled tree families[END_REF] (see also [START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF] for many additional pointers to the literature) once the required singular expansions of the involved generating functions are established. This includes returns to record-subtrees in Cayley trees, edge-cutting in Cayley trees, returns to zero in Dyck paths, cyclic points and trees in graphs of random mappings, all leading to (mixed Poisson) Rayleigh laws, as well as block sizes in k-Stirling permutations.

In the following we discuss several new results for the distribution of different parameters such as returns to zero and sign changes in walks and bridges with arbitrary steps, the number of subtrees satisfying some constraint in different fundamental families of trees, as well as the table sizes in the Chinese restaurant process, and the evolution of the number of balls in triangular urn models. These examples illustrate that composition schemes F = G(H) × M lead universally to beta-Mittag-Leffler distributions (Definition 3.13), which simplify into generalized Mittag-Leffler distributions if M (z) has a nonnegative singular exponent.

6.1. Core size of supertrees and a confluent example. Let C denote the family of plane trees (i.e., trees with all arities allowed and embedded into the plane) defined by

C = Z × SEQ(C), which translates to C(z) = 1 - √ 1 -4z 2 .
Then, following [40, pp. 412-414, 714], we consider supertrees, or "trees of trees", defined by

K = C (Z red + Z blue ) × C , which translates to K(z) = C 2zC(z) .
Note that this is a critical scheme as one has ρ C = 1/4 and τ

C = C(ρ C ) = 1/2.
The tree family K corresponds to trees where onto each node we graft a red or blue tree; one can also draw them like done in Figure 5. By Lagrange inversion, these supertrees K are thus enumerated by a neat combinatorial sum:

K n = ⌊n/2⌋ k=1 2 k n -k 2k -2 k -1 2n -3k -1 n -k -1 ,
thus the sequence K n for n ≥ 2 starts like 2, 2, 8, 18, 64, 188, 656, 2154, . . . , constituting sequence A168506 in the OEIS 12 .

FIG [START_REF] Banderier | Basic analytic combinatorics of directed lattice paths[END_REF]. A bicoloured supertree is a "tree of trees": Each node (here drawn in a potato shape) of some initial plane tree is replaced by a red or a blue node to which one attaches another plane tree.

By the Laplace method or by singularity analysis, this directly leads to the asymptotic expansion

K n ∼ 4 n 8Γ(3/4)n 5/4 .
(See also DeVries [START_REF] Devries | Algorithms for Bivariate Singularity Analysis[END_REF][START_REF] Pemantle | Analytic Combinatorics in Several Variables[END_REF] for an alternative approach to this expansion via multivariate analysis.) This asymptotic behaviour is noteworthy, because one sees here an unusual occurrence of the exponent - 5 4 , while most tree-like structures in combinatorics usually involve the exponent -3 2 . In fact, one could similarly define super-supertrees, super-super-supertrees, and so on, by further iterations of the critical scheme: C k+1 = C k 2ZC k with C 0 = C. This gives an interesting family of combinatorial structures whose asymptotic enumeration involves a dyadic exponent -1 -1/2 k+1 ; see [START_REF] Banderier | Formulae and asymptotics for coefficients of algebraic functions[END_REF] for a complete characterization of the possible singular exponents for N-algebraic functions (i.e., generating functions of any structure which can be generated by a context-free language).

With respect to supertrees, the critical scheme is

K(z) = G(H(z)), G(z) = C(z), H(z) = 2zC(z),
where H(z) has the following Puiseux expansion at z ∼ 1 4 :

H(z) ∼ 1 4 - 1 4 √ 1 -4z.
Now, we can study the core size X n via the bivariate generating function

K(z, u) = C u • 2zC(z) ,
as well as the number of H-components of size j, as captured by

K(z, v) = C 2zC(z) + (v -1)2c j-1 z j .
We can then apply our main Theorems 4.1 and 5.1 (with λ G = λ H = 1/2, τ H = 1 4 , c H = -1 4 , and κ = 1). This directly gives the following corollaries. COROLLARY 6.1. The core size X n in supertrees of size n has factorial moments given by

E(X s n ) ∼ n s/2 • µ s , µ s = Γ(s -1 2 )Γ(-1 4 ) Γ(-1 2 )Γ( s 2 - 1 4 )
.

The scaled random variable X n /n 1/2 converges in distribution with convergence of all moments to a c = -1/2 moment-tilted stable distribution, more precisely, a generalized Mittag-Leffler distribution of index 1/2:

X n n 1/2 d --→ m X, with X d = tilt -1 2 ML 1 2 d = ML 1 2 , - 1 4 .
Moreover, we have the local limit theorem

P{X n = x • n 1/2 } ∼ 1 n 1/2 • f X (x),
with f X (x) denoting the density of the random variable X.

Note that by Legendre's duplication formula one has

µ s = 2 s • Γ( s 2 + 1 4 ) Γ( 1 4 )
, so the random variable X can also be seen as equal in law to the chi distribution χ( 1 2 ) of parameter 1 2 , which is itself a generalized gamma distribution [START_REF] Banderier | Periodic Pólya urns, the density method and asymptotics of Young tableaux[END_REF].

COROLLARY 6.2. The number of coloured trees of size j -1 in supertrees of size n has factorial moments of mixed Poisson type given by

E(X s n,j ) = ξ s n,j • µ s (1 + o(1)), with ξ n,j = 2 • ( 1 4 ) j-1 c j-1 • n 1/2 and mixing distribution X = χ( 1 2
) with E(X s ) = µ s . Furthermore, the random variable X n,j possesses the three distinct asymptotic régimes of Theorem 5.1, with a phase transition at j = Θ(n 1/3 ).

We note in passing that a very similar result holds for the family of binary supertrees S, occurring in Bousquet-Mélou's study [START_REF] Bousquet-Mélou | Limit laws for embedded trees: applications to the integrated superBrownian excursion[END_REF] of the integrated super-Brownian excursion. This family is defined in terms of the family of complete binary trees B:

S = B(Z × B), B = Z + Z × B × B,
with initial values of S n given by 1, 1, 3, 8, 25, 80, 267, 911, . . . , constituting sequence A101490 in the OEIS. These functional equations indeed lead to

B(z) = 1- √ 1-4z 2 2z
, and to the following Puiseux expansion for S(z) = S( √ z):

S(z) ∼ 1 - √ 2(1 -4z) 1/4 + (1 -4z) 1/2 + . . . ;
thus leading to limit laws similar to the ones of supertrees in Corollaries 6.1 and 6.2.

A confluent example. Next, we consider pairs of supertrees in which we mark the number of H-trees in the first part of the pair. This translates to the scheme F (z, u) = G(uH(z))K(z) (with G, H, and K defined as in (6.1)), which is confluent because we have λ H λ G = λ M = 1 4 . This case is interesting as P{X n = k} (the probability that a pair of supertrees has k H-trees in its first part) converges to the sum of a discrete law and a continuous law. More precisely, as given by Theorem 4.4, and illustrated in Figure 6, the limit is a linear combination of a Boltzmann distribution, namely B C ( 1 4 ) 13 , and a generalized Mittag-Leffler distribution, namely ML 1 2 , -1 4 .

FIG [START_REF] Banderier | Random maps, coalescing saddles, singularity analysis, and Airy phenomena[END_REF]. The empirical distribution P{Xn = k} (drawn with red dots), and its theoretical limiting distribution (in blue). This blue curve is a linear combination of a discrete and a continuous distribution (the middle and right curves drawn in black):

1 2 B C 1 4 + 1 2 √ n ML 1 2 , -1 4 .
The blue theoretical limit curve agrees almost perfectly with the red empirical distribution even for small values of n (here, n = 500). 13 C denoting the Catalan generating function, it could be natural to call this Boltzmann distribution B C ( 1 4 ) the Catalan distribution, but this name is already used for some other distributions popping up in random matrix theory, and having moments related to the Catalan numbers, like the Marchenko-Pastur distribution. 6.2. Root degree and branching structure in bilabelled increasing trees. Bilabelled increasing trees are a natural generalization of increasing trees [START_REF] Kuba | A combinatorial approach to the analysis of bucket recursive trees[END_REF] where every node is assigned two labels instead of just one. General families of bilabelled trees are in bijection with increasing diamonds, which are a natural type of directed acyclic graphs; see [START_REF] Kuba | On bucket increasing trees, clustered increasing trees and increasing diamonds[END_REF] for the general statement and Figure 7 for a concrete example. Increasing diamonds model partial orders and their linear extensions, as well as computational processes and their executions in parallel computing [START_REF] Bodini | Increasing diamonds[END_REF]. They possess nice combinatorial properties and are enumerated by variants of hook-length formulas [START_REF] Kuba | Bilabelled increasing trees and hook-length formulae[END_REF][START_REF] Kuba | Combinatorial families of multilabelled increasing trees and hook-length formulas[END_REF][START_REF] Kuba | On bucket increasing trees, clustered increasing trees and increasing diamonds[END_REF].

Given a degree-weight sequence (φ j ) j≥0 , the corresponding degree-weight generating function is defined as φ(t) = j≥0 φ j t j . The family T of bilabelled increasing trees can be described by the following symbolic equation:

T = Z □ * Z □ * φ T ,
where Z denotes single unilabelled nodes, A □ * B denotes the boxed product (i.e., the smallest label is constrained to lie in the A component) of the combinatorial classes A and B, and

φ(A) = φ 0 • {ϵ} + φ 1 • A + φ 2 • A 2 + .
. . denotes the class containing all weighted finite labelled sequences of objects of A (i.e., each sequence of length k is weighted by φ k ; ϵ denotes the neutral object of size 0); see [START_REF] Flajolet | Analytic Combinatorics[END_REF]. Note that increasing diamonds are associated with φ(t) = 1 1-t . For the exponential generating function

T (z) = n≥1 T n z n
n! where n counts the number of labels, the above construction translates into

T ′′ (z) = φ(T (z)), T (0) = 0, T ′ (0) = 0. ( 40
)
We now focus on the case of 3-bundled bilabelled increasing trees; see Figure 8. FIG 8. An increasing plane-oriented recursive tree and the bijectively equivalent 3-bundled bilabelled increasing plane-oriented recursive tree. (The grey ellipses are just here to sketch the bijective correspondences.)

The family of 3-bundled trees is defined by Equation ( 40) with the following degree-weight generating function

φ(t) = 1 (1 -t) 3 = k≥0 k + 2 2 t k . ( 41
)
In other words, each node may have any number k of children, and the binomial indicates two bars between these children, thus creating 3 (possibly empty) sequences (or bundles) of children. From [START_REF] Kuba | Combinatorial families of multilabelled increasing trees and hook-length formulas[END_REF] we get the remarkably simple closed form

T (z) = 1 -1 -z 2 = n≥1 (2n -1)!!(2n -3)!! z 2n (2n)! , ( 42
)
where the double factorials are defined as

(2n -1)!! = n k=1 (2k -1) = (2n)! n! • 2 n ,
with initial values of T 2n given by 1, 3, 45, 1575, 99225, 9823275, 1404728325, . . . , constituting sequence A079484 in the OEIS. We are interested in the random variable X n counting the root degree of these 3-bundled bilabelled increasing trees of size n, under the uniform random tree model. Note that by definition there are no bilabelled trees with an odd number of labels, so T 2n+1 = 0 and, consequently, X 2n+1 = 0. In the following we use the notation R n = X 2n for the root degree.

By ( 40) and ( 41), the generating function

T (z, u) = n≥1 T n E(u Xn ) z n n! satisfies ∂ 2 ∂z 2 T (z, u) = φ uT (z) = 1 1 -u(1 - √ 1 -z 2 ) 3 .
Therefore, the Taylor expansion of T (z, u) starts as follows

T (z, u) = 1 + 3u z 2 2! + (3u + 9u 2 ) z 4 4! + (36u + 135u 2 + 540u 3 ) z 6 6! + . . .
Since a derivative with respect to z is simply a shift in the coefficient sequence of exponential generating functions, we obtain

E(u Rn+1 ) = (2n)! T 2n+2 [z n ] 1 1 -u(1 - √ 1 -z) 3 .
Now we observe, by Stirling's formula for the gamma function ( 29) and singularity analysis [START_REF] Banderier | Enumeration and asymptotics of restricted compositions having the same number of parts[END_REF] applied to [START_REF] Flajolet | General combinatorial schemas: Gaussian limit distributions and exponential tails[END_REF], that

T 2n+2 (2n)! ∼ 2 √ n √ π ∼ [z n ] 1 (1 -z) 3/2 . ( 43
)
This implies that, except for the non-standard shift in the random variable, the problem is equivalent (for first-order asymptotics) to the composition scheme

(1 -u(1 - √ 1 -z)) -3 , i.e., ρ H = 1, λ H = 1
2 , λ G = -3, and λ M = +∞ (as M (z) = 1 is entire). We note in passing that in this special case, it is also possible to obtain a quite simple closed-form expression for the probability mass function, as well as the (factorial) moments, due to the explicit expressions for the involved generating functions:

E(X s 2n+2 ) = (2n)! T 2n+2 [z n ] ∂ s u 1 1 -u(1 - √ 1 -z) 3 u=1
.

However, here we use our general scheme and Theorem 4.1. We apply Legendre's duplication formula and obtain

Γ(s + 3)Γ( 3 2 ) Γ(3)Γ( s+3 2 ) = 2 s • Γ s + 4 2 .
This leads to the following result.

COROLLARY 6.3. The random variable R n , counting the root degree in a random strict bilabelled increasing three-bundled tree with 2n labels, with tree generating function given by φ(t) = (1 -t) -3 , satisfies

E(R s n ) ∼ n s/2 • 2 s • Γ s + 4 2 .
The random variable R n /n 1/2 converges (in distribution and in moments) to a multiple of a chi-distributed random variable X d = χ(4), with four degrees of freedom:

R n n 1/2 d --→ m √ 2 • X, X d = χ(4).
We can refine the root degree by looking at the branching structure. We denote by R n,j = X 2n,j the random variable counting the number of branches (i.e., subtrees) with 2j labels, 1 ≤ j ≤ n, attached to the root:

R n = j≥1 R n,j .
Such random variables naturally arise in the context of the Chinese restaurant process [START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF][START_REF] Pitman | Exchangeable and partially exchangeable random partitions[END_REF][START_REF] Pitman | Combinatorial Stochastic Processes[END_REF] and generalized plane-oriented recursive trees [START_REF] Kuba | Analysis of label-based parameters in increasing trees[END_REF]. See also Feng et al. [START_REF] Su | Uniform recursive trees: branching structure and simple random downward walk[END_REF] for the analysis of the branching structure of recursive trees. The generating function

T (z, v) = n≥1 T n E(v Xn,j ) z n n! satisfies ∂ 2 ∂z 2 T (z, v) = φ T (z) -(1 -v) T 2j (2j)! z 2j .
Consequently,

E(u Rn+1,j ) = (2n)! T 2n+2 [z 2n ] 1 1 -(1 - √ 1 -z 2 ) -(1 -v) T2j (2j)! z 2j 3 .
We can use the asymptotics [START_REF] Fréchet | A proof of the generalized second limit theorem in the theory of probability[END_REF] and apply Theorem 5.1 to obtain the following result.

COROLLARY 6.4. The random variable R n,j counting the number of size 2j branches attached to the root in a random strict bilabelled increasing three-bundled tree with 2n labels, with tree generating function given by φ(t) = (1 -t) -3 , has factorial moments of mixed Poisson type,

E(R s n,j ) = ξ s n,j • E(X s )(1 + o(1)), with ξ n,j = √ 2 • T2j (2j)! • n 1/2
and mixing distribution X = χ(4). Furthermore, the random variable R n,j possesses the three distinct asymptotic régimes of Theorem 5.1, with a phase transition at j = Θ(n 1/3 ).

Bridge = Seq(minimal bridges)

Walk never returning to 0

FIG 9.
A walk consists of an initial bridge which contains all returns to zero (red dots) and a final walk never returning to zero. The bridge is further decomposed into minimal bridges touching zero only twice.

6.3. Returns to zero: walks and bridges with drift zero. A lattice path of length n is a sequence (s 1 , . . . , s n ) of steps s i ∈ S for a fixed finite subset S ⊆ Z called step set. Geometrically, we fix the starting point 0 and consider the partial sums k i=1 s i which can be interpreted as appending the steps one after another. Each step s i gets a weight p i > 0 and the weight of a path is the product of all steps. The step polynomial P (u) = i p i u i connects the weights and the steps. We call a step set periodic if there exist integers b, p ∈ Z, p ≥ 2 such that P (u) = u b P (u p ); otherwise we call it aperiodic. Here and in the next section, we assume that the step set is aperiodic. Note that this is no major constraint as the asymptotics of walks with periodic steps can be deduced from the ones with aperiodic ones; see [START_REF] Banderier | The kernel method for lattice paths below a rational slope[END_REF]. We call a lattice path a walk if it is unconstrained, and a bridge if it ends at zero, i.e., n i=1 s i = 0. A return to zero is a point in the path such that k i=1 s i = 0; see Figure 9. Generalizing results from Feller [START_REF] Feller | An Introduction to Probability Theory and its Applications[END_REF] and Barton [START_REF] Skellam | Distributions associated with random walk and recurrent events[END_REF]Discussion p. 115], it was shown in [118, Section 3.2] that for drift P ′ (1) = 0 the law of the number of returns to zero follows a Rayleigh distribution for bridges, while it follows a half-normal distribution for walks. The proof utilized a general theorem on the singular structure of the generating functions. However, the situation is also amenable to Theorem 4.1 and we can give an alternative proof next.

Let w n,k be the number of walks of length n with k returns to zero. The bivariate generating function of walks W (z, u) = n,k≥0 w n,k z n u k is given by

W (z, u) = 1 1 -u 1 -1 B(z) W (z) B(z) , (44)
where B(z) and W (z) = 1/(1 -zP (1)) are the generating functions of bridges and walks, respectively; see [START_REF] Wallner | A half-normal distribution scheme for generating functions[END_REF]Equation (3.3)]. To explain [START_REF] Fuchs | Subtree sizes in recursive trees and binary search trees: Berry-Esseen bounds and Poisson approximations[END_REF], observe that every bridge is a sequence of minimal bridges, which are bridges that never return to the x-axis between the start-and endpoint; see Figure 9. Therefore, minimal bridges are enumerated by 1 -1/B(z). Hence, this is exactly the situation of the extended composition scheme [START_REF] Aldous | Exchangeability and related topics[END_REF] with

G(z) = 1/(1 -z), H(z) = 1 -1/B(z)
, and M (z) = W (z)/B(z). Now, for zero drift, [START_REF] Wallner | A half-normal distribution scheme for generating functions[END_REF]Equation (3.4)] shows that

B(z) = c B 1 -zP (1) + O(1) with c B = P (1) 2P ′′ (1) . ( 45
)
Hence, one has λ G = -1, λ H = 1 2 , and λ - M = -1 2 , therefore we are in the pure régime; see Definition 3.3. This is exactly the situation in Remark 4.2 and the number of returns to zero in walks thus follows a half-normal distribution with parameter σ = √ 2c B = P (1)/P ′′ (1). Now, the generating function for bridges is nearly the same as W (z, u) from [START_REF] Fuchs | Subtree sizes in recursive trees and binary search trees: Berry-Esseen bounds and Poisson approximations[END_REF] except that the last factor W (z)/B(z) is omitted. So, by Corollary 4.2, the number of returns to zero here follows a Rayleigh distribution with the same parameter σ = P (1)/P ′′ (1).

We can refine this result for the random variable X n,j counting the number of distance-jzeroes (which were introduced in [START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF]). These are the number of returns to zero which have a distance of exactly j steps to the previous zero contact. The union over j of distance-j-zeroes gives all returns to zero, and they therefore clearly represent a partition of all returns to zero. Using Theorem 5.1, we then get the following limit theorem. COROLLARY 6.5. Let X n,j be the number of distance-j-zeroes in lattice paths of length n. For walks (resp. bridges) with zero drift (i.e., P ′ (1) = 0), X n,j has factorial moments of mixed Poisson half-normal type (resp. mixed Poisson Rayleigh type)

E(X s n,j ) = ξ s n,j • E(X s ) (1 + o(1)) , (46)
with ξ n,j =

P (1) 2P ′′ (1)
hj

P (1) j • n 1/2
, where X is given by

X = HN(σ)
for walks, Rayleigh(σ) for bridges, σ = P (1) P ′′ (1) .

Furthermore, the random variable X n,j possesses the three distinct asymptotic régimes of Theorem 5.1, with a phase transition at j = Θ(n 1/3 ). REMARK 6.6 (Universality of the rescaling factor). Note that the rescaling factor ξ n,j in ( 46) is the same for walks and bridges, while the moment sequence µ s changes. This independence of ξ n,j can be explained: In Figure 9, the last factor of the walk is a walk not touching zero and is encoded by M (z) = W (z)/B(z). Then by Theorem 5.1 we know that ξ n,j is independent of this factor, and thus has the same value for walks and bridges.

Furthermore, this factor M (z) is also responsible for the often observed dichotomy between half-normal and Rayleigh distributions in the extended composition scheme which we will also observe in the next examples of initial returns and sign changes.

Note that Formula (46) offers a neat factorization for the moments: One can regroup in one factor the quantities with a probabilistic flavour (involving the variance P ′′ (1) of the allowed steps, and µ s ), while the remaining factor (h j , the number of minimal bridges of length j) corresponds to a quantity with a combinatorial flavour. This could also be explained using renewal theory. 6.4. Initial returns in coloured bridges. We generalize the previous model by introducing m-coloured bridges B m (see Figure 10): We append m non-empty bridges one after the other (and each one with a different colour): B m = (B -1) m . Then, we are interested in the number of returns to zero in the first bridge, i.e., the initial one that we uniformly coloured. We call such returns the initial returns. Reusing the combinatorial constructions of the previous section, this gives for the bivariate generating function B m (z, u) the following decomposition

B m (z, u) =   1 1 -u 1 -1 B(z) -1   (B(z) -1) m-1 . (47)
For m = 1 this is (44) except the factor W (z)/B(z) and the constraint to be non-empty. Asymptotically, and therefore for the law, the non-emptiness is negligible. The generating function W m (z, u) of m-coloured walks (m-tuples of bridges with a few more steps coloured in the same colour as the final bridge) is given by

W m (z, u) = (1 + B m (z, u)) W (z) B(z) .
Now, we can directly apply Theorem 4.1. From the reasoning above we see that

λ G = -1, λ H = 1 2
, and λ M = -m-1 2 for bridges and λ M = -m 2 for walks.

COROLLARY 6.7. The random variable X n counting the number of initial returns in a m-coloured walk (resp. bridge) of length n satisfies

E(X s n ) ∼ n s/2 σ √ 2 n µ s , σ = P (1) P ′′ (1) , µ s = Γ(s+1)Γ((m+1)/2)
Γ((m+s+1)/2) , for walks,

Γ(s+1)Γ(m/2) Γ((m+s)/2) , for bridges.
The scaled random variable X n /n 1/2 converges in distribution with convergence of all moments to a Rayleigh distribution and a scaled beta distribution (see Definition 3.12 and Example 3.18):

X n n 1/2 d --→ m X, X d = Rayleigh(σ) • B 1/2 ,
with independent Rayleigh(σ) and B = Beta 1 2 , m 2 , for walks, Beta 1 2 , m-1

2

, for bridges, where Beta(α, 0) = 1. Moreover, we have the local limit theorem

P{X n = x • n 1/2 } ∼ n -1/2 • f X (x),
with density f X (x) of the random variable X for bridges (for walks one replaces m by m + 1) given by

f X (x) = 2 πσ 2 Γ m 2 e -x 2 2σ 2 U m 2 -1, 1 2 , x 2 2σ 2 ,
where U (a, b, x) is the confluent hypergeometric function of the second kind which is the solution of

zy ′′ + (b -z)y ′ -ay = 0 such that U (a, b, x) ∼ z -a for z → ∞ and | arg(z)| < 3π/2; see [25, Section 13.2].
Observe the special cases U (-1/2, 1/2, x) = √ x and U (0, 1/2, x) = 1 which nicely give the density functions of a Rayleigh (see Example 3.18) and a half-normal distribution (see Example 3.17). Hence, for m = 1 we recover the results of the previous section and uncover a large family of connected probability distributions. It is interesting that this family also appears in the context of preferential attachments in graphs [104, Formula (1.1)].

It is also interesting to consider multicoloured bridges, where we allow any number of colours. We still mark by u the initial returns. The corresponding bivariate generating function B(z, u) is given by

B(z, u) = m≥1 B m (z, u) =   1 1 -u 1 -1 B(z) -1   1 2 -B(z)
.

The generating function for the number of multicoloured bridges is thus B(z, 1) = 1 2-B(z) -1. From [START_REF] Fuchs | Limit theorems for subtree size profiles of increasing trees[END_REF] we see that B(z) possesses a singularity of order -1/2 at z = 1/P (1), and hence B(z, 1) becomes singular at some z 0 > 0 which is the unique solution of B(z 0 ) = 2. Hence, the probability generating function reveals a geometric distribution of parameter 1/2:

[z n ]B(z, u) [z n ]B(z, 1) ∼ 1 1 -u 1 -1 B(z0) -1 = u/2 1 -u/2 .
As the truncated sum m0 m=1 B m (z, u) behaves asymptotically like B m0 (z, u), we see here a phase transition from a continuous law (for any finite m 0 ) to a discrete law (when m 0 goes to infinity). Note that this phenomenon holds verbatim for walks.

Finally, let us apply the size-refined scheme Theorem 5.1, counting initial returns in mcoloured bridges (or walks) which are a certain distance apart: COROLLARY 6.8. Let X n,j be the number of initial returns at distance j from the previous zero in m-coloured walks or bridges of length n. Then, X n,j has moments of mixed Poisson type

E(X s n,j ) = ξ s n,j • E(X s ) (1 + o(1)) , with ξ n,j = P (1) 2P ′′ (1)
hj

P (1) j • n 1/2 , h j = [z j ](1 -1/B(z))
, and where X is the product of a beta and a Rayleigh distribution given in Corollary 6.7. Furthermore, X n,j possesses the three distinct asymptotic régimes of Theorem 5.1, with a phase transition at j = Θ(n 1/3 ).

These results also hold for other variants of paths, in bijection with sequences that already appeared in the literature; see Table 3. 6.5. Sign changes in walks. Using the same notation as in Example 6.3, we now define the sign of the path (s 1 , . . . , s n ) after k steps as sgn( k i=1 s i ) ∈ {-1, 0, 1}. Thereby every lattice path is associated with a sequence of signs. A sign change is therein any subsequence (-1, 0 * , 1) or (1, 0 * , -1) where 0 * denotes a (possibly empty) sequence of 0s; see Figure 11.

In this section we consider Motzkin paths. They are composed of up steps +1, down steps -1, and horizontal steps 0; see again Figure 11. Their step polynomial is therefore given by P (u) = p-1 u + p 0 + p 1 u (with p -1 p 0 p 1 ̸ = 0). In the case of zero drift, let us show how to apply our results to get that the number of sign changes follows asymptotically a Rayleigh distribution for bridges and a half-normal distribution for walks, while for nonzero drift it follows a geometric distribution; see [START_REF] Wallner | A half-normal distribution scheme for generating functions[END_REF].

Steps GF Sequence OEIS {U, D} 8z 2 -2- √ 1-4z 2 16z 2 -3
1, 0, 2, 0, 10, 0, 52, 0, 274, 0, 1452, 0, 7716, . . . A075436

{U, D, H 1 } z+ √ 1-4z 2 1-5z 2
1, 1, 3, 5, 13, 25, 61, 125, 295, 625, 1447, . . .

A098615 {U, D, H 2 } z 2 + √ 1-4z 2 1-4z 2 -z 4
1, 0, 3, 0, 11, 0, 43, 0, 173, 0, 707, 0, 2917, . . .

A026671 TABLE 3

Multicoloured bridge models: they end at 0 and use up steps U = (1, 1), down steps D = (1, -1), and horizontal steps H i = (i, 0) allowed only at altitude 0. The limit laws of initial returns to zero in these models are all the same and special cases of Corollaries 6.7 and 6.8.

FIG 12.

A plane-oriented recursive tree of size 15, where the root node labelled zero has one size-one, one size-two, one size-three, and two size-four branches, as well as the corresponding table structure in the Chinese restaurant model.

6.6. Tables in the Chinese restaurant process. Following Aldous, Pitman, and Dubins (see [START_REF] Aldous | Exchangeability and related topics[END_REF][START_REF] Pitman | Exchangeable and partially exchangeable random partitions[END_REF][START_REF] Pitman | Combinatorial Stochastic Processes[END_REF]), we now consider the Chinese restaurant process. This a discrete-time stochastic process having as value at time n one of the B n partitions of the set [n] = {1, 2, . . . , n} (where B n denotes the Bell numbers found as sequence A000110 in the OEIS). One fancifully imagines a Chinese restaurant with an infinite number of tables, where each table has a possibly infinite number of seats. In the beginning the first customer takes a seat at the first table. At each discrete time step a new customer arrives and either joins one of the existing tables, or takes a seat at the next empty table. Each table corresponds to a block of a random partition. The process thus starts at time n = 1 with the partition {{1}}. Now, given a partition T = {t 1 , . . . , t k } of [n] with |T | = k parts t i , at time n + 1 the element n + 1 is either added to one of the existing parts t i ∈ T with probability

P{n + 1 ≺ t i } = |t i | -α n + θ , 1 ≤ i ≤ k,
where n + 1 ≺ t i denotes that n + 1 is a costumer sitting at table t i , or as a new singleton block with probability

P{n + 1 ≺ t |T |+1 } = θ + k • α n + θ .
This model (parametrized by the two parameters 0 < α < 1 and θ > -α) thus assigns a probability to any particular partition T of [n]. We are interested in the random variable C n , counting the total number of tables in the Chinese restaurant process, as well as the random variable C n,j , counting the number of parts of size j in a partition of [n]. As pointed out in [START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF], this process can be embedded into a variant of the growth process of generalized planeoriented recursive trees with two different connectivity parameters a > 0 and b > -1. This allows us to study properties of the Chinese restaurant process using analytic combinatorial tools. For the reader's convenience, we restate this embedding below.

We collect the results of [START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF] and complement them by extending the constraint b > 0 to the full range b > -1 as well as by providing the missing identification of the limit law as a (moment-shifted) stable law.

Combinatorially, we consider a family T a,b of generalized plane-oriented recursive trees, where the degree-weight generating function ψ(t) = 1

(1-t) b , b > 0, associated with the root of the tree, is different to the one for non-root nodes in the tree, φ(t) = 1

(1-t) a , a > 0. Then, the family T a,b is closely related to the corresponding family T of generalized plane-oriented recursive trees with degree-weight generating φ(t) = 1

(1-t) a , a > 0, via the following formal recursive equations (see Section 6.2 for the definition of the boxed product):

T a,b = Z □ * ψ(T ), T = Z □ * φ(T ).
The weight w(T ) of a tree T ∈ T a,b is then defined by

w(T ) = ψ d(root) v∈T \{root} φ d(v) ,
where d(v) denotes the outdegree of node v. Thus, the generating functions

T a,b (z) = n≥1 T a,b;n z n n! and T (z) = n≥1 T n z n n!
of the total weight of size-n trees in T a,b and T , respectively, satisfy the differential equations T ′ a,b (z) = ψ(T (z)), T ′ (z) = φ(T (z)). The ordinary tree evolution process to generate a random tree of arbitrary given size in the family T (see [START_REF] Panholzer | Level of nodes in increasing trees revisited[END_REF] for a detailed discussion) can be extended in the following way to generate a random tree in the family T a,b . The process, evolving in discrete time, starts with the root labelled zero. At step n + 1, with n ≥ 0, the node with label n + 1 is attached as a new child to any previous node v (this is denoted by n + 1 ≺ v) with probabilities

P{n + 1 ≺ v} = d(v)+b b+(a+1)n if v is the root, d(v)+a b+(a+1)n if v is not the root.
We recall the following result from [START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF]. PROPOSITION 6.10 (Chinese restaurant process and generalized plane-oriented recursive trees). A random partition of {1, . . . , n} generated by the Chinese restaurant process with parameters α > 0 and θ > 0 can be generated equivalently by the growth process of the family of generalized plane-oriented recursive trees T a,b when generating such a tree of size n + 1. The parameters α, θ and a, b > 0, respectively, are related via

α = 1 1 + a , θ = b 1 + a .
The random variable C n is distributed as the outdegree X n+1 of the root of a random generalized plane-oriented recursive trees of size n + 1 from the family T a,b :

C n d = X n+1 .
The random variable C n,j is distributed as X n+1,j , which corresponds to the number of branches of size j attached to the root of a random tree of size n + 1 from the family T a,b :

C n,j d = X n+1,j .
Note that in the above relation, θ cannot be negative, since b is assumed to be positive. As already observed in [START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF], the correspondence can be extended to the full range 0 < α < 1 and θ > -a, where one has a = 1 α -1 > 0 and b = θ α > -1. For -1 < b ≤ 0, we cannot directly use the degree-weight generating function ψ(t) = (1 -t) -b . Indeed, for -1 < b < 0 we would have ψ(t) = 1 + bt + . . . , involving a negative weight; while for b = 0 we would have ψ(t) = 1, a degenerate case. However, we can use a modified generating function, leading to a correct model of the Chinese restaurant process in the range -1 < b ≤ 0 (see [START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF][START_REF] Panholzer | Level of nodes in increasing trees revisited[END_REF] for more details on the growth process):

ψ(t) = 1 + t 0 1 (1 -x) 1+b dx = 1 + 1 b 1 (1 -t) b -1 = 1 + k≥1 b + k k -1 t k k ,
for -1 < b < 0, while for b = 0 one uses

ψ(t) = 1 -log(1 -t) = 1 + k≥1 t k k .
Thus, we have some generalized plane-oriented recursive trees attached to a root with a different tree-weight generating function ψ(t). Summarizing, we have

ψ(t) =        1 (1-t) b if b > 0, 1 -log(1 -t) if b = 0, 1 + 1 b 1 (1-t) b -1 if -1 < b < 0.
Here (except for the special case b = 0, which is handled by a slightly different approach, detailed later in Theorem 7.1), we can directly apply our results from Theorem 4.1 to

T ′ a,b (z, u) = n≥1 T n E(u Xn ) z n-1 (n -1)! = ψ(u • T (z)), T ′ (z) = φ(T (z)),
for the total number of tables, and from Theorem 5.1 to

R ′ a,b (z, v) = n≥1 T n E(v Xn,j ) z n-1 (n -1)! = ψ(T (z) -(1 -v)z j T j j! ),
for the number of tables of size j. This allows us to extend the corresponding result of [START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF] to the full range of b > -1, also providing the missing identification of the limit law as a (moment-tilted) stable law: THEOREM 6.11. Let a > 0, b > -1. The random variable X n,j counting the number of branches of size j in a random T a,b tree of size n (or, equivalently, the number of tables with j seated customers in a Chinese restaurant process of parameter α = 1/(1 + a) and θ = b/(1 + a), with a total of n -1 customers) possesses the three distinct asymptotic régimes of Theorem 5.1, with a phase transition at j = Θ(n 1/(a+2) ):

(i) For j ≪ n 1 a+2 we have ξ n,j = αn α j j-1-α j-1
→ ∞ and Xn,j ξn,j converges in distribution, with convergence of all moments, to a generalized Mittag-Leffler distribution:

X n,j ξ n,j d --→ m X with X d = ML(α, θ).
(ii) For j ∼ r • n 1 a+2 , r ∈ (0, ∞), we have ξ n,j → ξ, and the random variable X n,j converges in distribution, with convergence of all moments, to a mixed Poisson distribution:

X n,j d --→ m MPo(ξX).

(iii) For j ≫ n 1 a+2 we have ξ n,j → 0, so X n,j converges to a Dirac distribution at 0. REMARK 6.12. Our result above implies that there are only a few giant tables in the Chinese restaurant process (a mixed-Poisson number of tables with a number of customers proportional to n 1 a+2 ). In contrast, there are much more tables with a smaller number of customers, and an asymptotically negligible number of tables of size ≫ n 1 a+2 ; see Figure 13. REMARK 6.13. Closed formulas for the factorial moments E(X s n,j ), as well as a formula for the probability mass function of X n,j are readily obtained from the generating functions by extraction of coefficients. REMARK 6.14. Our results also allow recovering the limit theorem in [START_REF] Pitman | Combinatorial Stochastic Processes[END_REF] for the total number of tables C n in the Chinese restaurant process (via X n ), albeit with a totally different proof, as the normalized random variable X n /n α converges in distribution with convergence of all moments to a random variable X, with X given in the theorem before. For the reader's

E(X s ) = Γ(s + θ α )Γ(θ) Γ(θ + s • α)Γ( θ α )
.

PROOF OF THEOREM 6.11. We follow very closely [START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF] and sketch the remaining steps. We solve the differential equation T ′ (z) = φ(T (z)) and get

T (z) = 1 -(1 -(a + 1)z) 1 a+1 .
Thus, the probability generating function is given by

E(v Xn+1,j ) = n! T n+1 [z n ]ψ T (z) -(1 -v)z j T j j! ,
where the coefficient

T n+1 n! = [z n ]ψ T (z)
is computed by standard singularity analysis. Therefore, except for the non-standard shift, we can readily apply our scheme to the generating function

ψ T (z) -(1 -v)z j T j j! .
More precisely, since a > 0 and b > -1, we have λ H = 1 a+1 and λ G = -b for b ̸ = 0. Hence, the critical range is given by

j(n) = Θ n λ H 1+λ H = Θ n 1 a+2
.

Here, no additional factor M (z) is present, so λ - M = 0. In the case of b = 0 we apply the cycle scheme of Theorem 7.3. This gives

E(X s ) = µ s =    Γ(s+b)Γ( b a+1 ) Γ(b)Γ( b+s a+1 ) if β ̸ = 0, Γ(s+1) Γ( s a+1 +1) if β = 0.
Finally, we unify both expressions by simply using Γ(x + 1) = xΓ(x).

6.7.

Triangular urn models and beta-Mittag-Leffler distributions. Two-colour triangular urns are instances of generalized Pólya urn models [START_REF] Flajolet | Some exactly solvable models of urn process theory[END_REF][START_REF] Janson | Functional limit theorems for multitype branching processes and generalized Pólya urns[END_REF][START_REF] Hosam | Pólya Urn Models[END_REF]. At each time step n ≥ 1, a ball is drawn uniformly at random, reinserted, and depending on the observed colour, balls of both colours are added to the urn: If a white ball was drawn, we add α white and β black balls, whereas, if a black ball was drawn, we add γ white and δ black balls. The addition/replacement of balls can be described by the so-called ball replacement matrix

M = α β γ δ ,
where for balanced urn models it holds that α + β = γ + δ, such that the total number σ = α + β of added balls in each step is independent of the observed colour. The initial configuration of the urn consists of w 0 white balls and b 0 black balls, and the random variable W n counts the number of white balls in the urn after n draws. For balanced triangular urns with replacement matrix

M = α β 0 δ , α, β > 0, δ = σ = α + β,
it was shown by Flajolet, Dumas, and Puyhaubert [START_REF] Flajolet | Some exactly solvable models of urn process theory[END_REF] (and also by Janson [START_REF] Janson | Limit theorems for triangular urn schemes[END_REF][START_REF] Janson | Moments of Gamma type and the Brownian supremum process area[END_REF] via different analytic methods) that Wn αn α/σ d -→ W, for a random variable with moments

E(W s ) = Γ( b0+w0 σ ) Γ( w0 α ) • Γ(s + w0 α ) Γ(s • α σ + b0+w0 σ )
.

Figure 14 shows the evolution of the first three draws of a triangular urn.

In the special case (w 0 , b 0 ) = (α, β), and thus w 0 + b 0 = σ, this simplifies to a Mittag-Leffler distribution with parameter 0 < α/σ < 1. As it is often the case for urn models, the limit law strongly depends on the initial composition of the urn. For b 0 > 0 and either w 0 = 0 or w 0 = β, Janson observed a moment-tilted stable law, leaving the other cases open; see [START_REF] Janson | Limit theorems for triangular urn schemes[END_REF]Theorem 1.8 and Problem 1.15]. We can directly obtain the limit law W using our extended scheme. The key tool is the history generating function F (z, u) = n,k≥0 f n,k u k z n n! where f n,k is equal to the number of transitions (or histories) leading to a configuration with k white balls after n steps; see Figure 14. The closed form of this history generating function was derived in [START_REF] Flajolet | Some exactly solvable models of urn process theory[END_REF]Proposition 14]:

F (z, u) = u w0 (1 -σz) -b0/σ 1 -u α 1 -(1 -σz) α/σ -w0/α .
Putting aside the prefactor u w0 , and after a change of variable u α → u, this equation can be interpreted as an extended critical composition scheme

F (z, u) = M (z) • G uH(z) , (48)
involving the exponential generating functions with nonnegative integer coefficients M (z) = (1 -σz) -b0/σ , G(z) = (1 -αz) -w0/α , and

H(z) = (1 -(1 -σz) α/σ )/α. (49)
The fact that the singular exponents depend on b 0 and w 0 explains en passant why the limit distribution of W n differs according to the initial composition of the urn. Indeed, as the number of white balls at time n satisfies

P{W n = αk + w 0 } = n![z n u k ]F (z, u) n![z n ]F (z, 1) = g k k! [z n ]H(z) k M (z) [z n ]F (z, 1) ,
we can apply Theorem 4.1 and we then get the following limit distributions of W n for balanced triangular urn models, completing and extending earlier results [64, Theorem 1.8]: COROLLARY 6.15. Let W n be the random variable for the number of white balls in a balanced triangular Pólya urn with initially w 0 > 0 white and b 0 ≥ 0 black balls. Then, we have a convergence in distribution, with convergence of all moments, towards a beta-Mittag-Leffler distribution (see Definition 3.13)

W n αn α/σ d --→ m BML α σ , w 0 α , b 0 α .
REMARK 6.16 (Almost sure convergence and beyond). This limit was also recently identified by Goldschmidt et al. [START_REF] Goldschmidt | Stable graphs: distributions and line-breaking construction[END_REF] for urns with non-integer weights: A link with the Chinese restaurant model (for b 0 = 0) leads to a Mittag-Leffler distribution, then they show that the impact of b 0 > 0 on the process leads to a distribution with an additional beta law factor. It is interesting to stress that their approach implies an almost sure convergence. Note that the fluctuations around this almost sure limit are known: A second-order central limit theorem (that is, the random variable minus its almost sure limit converges, rescaled, to a Gaussian distribution), as well as a law of the iterated logarithm was obtained using discrete martingale [START_REF] Kuba | On martingale tail sums in affine two-color urn models with multiple drawings[END_REF]. What is more, following Gouet [START_REF] Gouet | Martingale Functional Central Limit Theorems for a Generalized Pólya Urn[END_REF], a continuous-time reparametrization leads to a functional second-order limit theorem for balanced urn models. There is currently no systematic way to obtain an almost sure convergence for all combinatorial models covered by our composition schemes; however, in a few cases (e.g., for walks, trees, and maps), some ad-hoc clever constructions entail this almost sure convergence [88-90, 93, 113].

Note that applying Theorem 5.1 to the size-refined version of the composition scheme (48), we get factorial moments of mixed Poisson type for a size-refined random variable X n,j and the corresponding limit laws. However, the combinatorial interpretation of the random variable(s) X n,j is more involved and will be given elsewhere.

We stress the fact that the methods and results presented in this article are thus holding both for ordinary generating functions (typically used for unlabelled structures) and for exponential generating functions (typically used for labelled structures); see [START_REF] Flajolet | Analytic Combinatorics[END_REF].

This concludes the list of applications for our results on the extended and size-refined composition schemes. We now give some extensions of our work to other schemes. THEOREM 7.4. In a multivariate pure extended critical composition scheme (52), the joint moments of the random vector X n = (X n,1 , . . . , X n,m ) satisfy E(X

= Γ -m ℓ=1 λ Gℓ λ Hℓ -λ M Γ m ℓ=1 s ℓ λ Hℓ -m ℓ=1 λ Gℓ λ Hℓ -λ M m ℓ=1 Γ(s ℓ -λ Gℓ ) Γ(-λ Gℓ ) . (53)
Consequently, one gets a convergence in distribution and in moments

X n,1 κ 1 n λ1 , . . . , X n,m κ m n λm d --→ m X, ( 54
)
where X is determined by its joint moment sequence µ s = µ s1,...,sm . Moreover, the random vector X has a scaled Dirichlet-stable product distribution,

X = (X 1 , . . . , X m ) d = (V 1 • W λH 1 1 , . . . , V m • W λH m m ), (55) where
W = (W 1 , . . . , W m , W m+1 ) follows a Dirichlet distribution W d = Dir(-λ G1 λ H1 , . . . , -λ Gm λ Hm , -λ M ), (56)
and where the V ℓ 's, for 1 ≤ ℓ ≤ m, are m independent generalized Mittag-Leffler distributions V ℓ d = ML(λ Hℓ , -λ Gℓ λ Hℓ ), also independent of W.

PROOF. We proceed similarly to the proof of the first part of Theorem 4.1. First, the mixed factorial moments of X n , which are obtained by differentiation and extraction of coefficients:

E X s1 n,1 • • • X sm n,m = [z n ]∂ s1 u1 . . . ∂ sm um (F )(z, 1, . . . , 1) [z n ]F (z, 1, . . . , 1)
.

The differentiation with respect to u ℓ only affects the factor G ℓ u ℓ H ℓ (z) , leading to a singular expansion covered in Section 3.1. Extraction of coefficients then gives an asymptotic expansion of the factorial moments. Converting all the factorial moments into moments using [START_REF] Drmota | Marking in combinatorial constructions: generating functions and limiting distributions[END_REF] gives the desired asymptotics and moments in [START_REF] Humbert | Nouvelles correspondances symboliques[END_REF]. It remains to identify the distribution. To this aim, note that a Dirichlet distributed random vector (W 1 , . . . , W m+1 ) d = Dir(a 1 , . . . , a m+1 ) with positive parameters a 1 , . . . , a m+1 has a density function supported on the m simplex {(x 1 , . . . , x m+1 ) ∈ R m+1 ≥0 | m+1 j=1 x j = 1}:

f (x 1 , . . . , x m+1 ) = Γ( m+1 j=1 a j m+1 j=1 Γ(a j ) m+1 j=1 x aj-1 j
.

Accordingly, the joint moments are given by

E W s1 1 • • • W sm+1 m+1 = Γ( m+1 j=1 a j ) Γ m+1 j=1 (s j + a j ) m+1 j=1 Γ(s j + a j ) m+1 j=1 Γ(a j ) . Now, consider a random vector (Z 1 , . . . , Z m ) satisfying (Z 1 , . . . , Z m ) d = (V 1 • W α1 1 , . . . , V m • W αm m ) with generalized Mittag-Leffler distributions V ℓ d = ML(α ℓ , a ℓ ) for ℓ = 1, . . . , m,
such that all random variables are mutually independent. Using the closed form [START_REF] Bodini | Increasing diamonds[END_REF] we get

E Z s1 1 • • • Z sm m = Γ( m+1 ℓ=1 a ℓ ) Γ a m+1 + m ℓ=1 (α ℓ s ℓ + a ℓ ) m ℓ=1 Γ(s ℓ + a ℓ /α ℓ ) Γ(a ℓ /α ℓ ) .
Comparing this expression with the moments (53), the claim follows.

REMARK 7.5. The marginals X ℓ of the random vector X are also covered by Theorem 4.1. The random vector X is closely related to Poisson-Dirichlet distributions PD(α, θ), [START_REF] Fitzgerald | Stick-breaking P G(α, ζ)-generalized gamma processes[END_REF][START_REF] Fitzgerald | Generalized Mittag Leffler distributions arising as limits in preferential attachment models[END_REF][START_REF] Pitman | The two-parameter Poisson-Dirichlet distribution derived from a stable subordinator[END_REF] and the joint limit law of node degrees in preferential attachment trees or generalized plane-oriented recursive trees [START_REF] Tamás | The maximum degree of the Barabási-Albert random tree[END_REF]; see also the subsequent example. Now, if one considers the multivariate size-refined scheme

F = M × m ℓ=1 G ℓ H ℓ,̸ =jℓ + v ℓ H ℓ,=jℓ ,
one gets the following multivariate version of Theorem 5.1. THEOREM 7.6 (Multivariate pure size-refined critical scheme). In a multivariate pure size-refined critical composition scheme

F (z, v 1 , . . . , v m ) = M (z) • m ℓ=1 G ℓ H ℓ (z) -(1 -v ℓ )h ℓ,jℓ z jℓ
the random variables X n,ℓ,jℓ , which count the number of H ℓ -components of size j ℓ , have joint factorial moments of mixed Poisson type:

E X s1 n,1,j1 . . . X sm n,m,jm = µ s1,...,sm • m ℓ=1 ξ sℓ n,ℓ,jℓ • (1 + o(1)), with ξ n,ℓ,jℓ = ρ j ℓ H ℓ
-cH ℓ h ℓ,jℓ n λH ℓ and joint mixing distribution X = (X 1 , . . . , X m ) as in Equation [START_REF] Hwang | Asymptotics of Poisson approximation to random discrete distributions: an analytic approach[END_REF].

Let X ℓ , for 1 ≤ ℓ ≤ m, denote the marginal distribution of the ℓth coordinate of X = (X 1 , . . . , X m ). For n → ∞, the limiting distributions of X n,ℓ,jℓ jointly undergo mixed Poisson type phase transitions with mixing distributions X ℓ . The phase transitions depend on the growth of j ℓ = j ℓ (n), with critical growth ranges given by j

ℓ = j ℓ (n) = Θ(n λ H ℓ 1+λ H ℓ ).
In particular, for j ℓ (n) ∼ ξ ℓ • n λ H ℓ 1+λ H ℓ , the random vector X n,j = (X n,1,j1 , . . . , X n,m,jm ) converges in distribution with convergence of all (factorial) moments to a multivariate distribution MPo(ξX).

For properties of multivariate mixed Poisson distributions we refer to [START_REF] Ferrari | Multivariate mixed Poisson distributions[END_REF] or [START_REF] Kuba | On moment sequences and mixed Poisson distributions[END_REF]. The key tool for proving Theorem 7.6 is the following multivariate extension of Lemma 5.4. LEMMA 7.7 (Joint factorial moments and limit laws of mixed Poisson type). Let (X n ) n∈N denote a sequence of m-dimensional random vectors, whose factorial moments are asymptotically of mixed Poisson type, i.e., they satisfy for n → ∞ the asymptotic expansion

E(X s n ) = E X s1 n,1 . . . X sm n,m = µ s1,...,sm • m ℓ=1 ξ sℓ n,ℓ • (1 + o(1)), s ≥ 1,
with µ s1,...,sm ≥ 0, and ξ n,ℓ > 0 for 1 ≤ ℓ ≤ m. Furthermore, assume that the sequence of joint moments (µ s ) s∈N m determines a unique distribution L = (L 1 , . . . , L m ). Then, for n → ∞, one has the following joint limit distributions:

(i) If ξ n,ℓ → ∞, the random variable Xn,ℓ ξn,ℓ converges in distribution, with convergence of all moments, to L ℓ . (ii) If ξ n,ℓ → ξ ∈ (0, ∞), the random variable X n,ℓ converges in distribution, with convergence of all moments, to Y

d = MPo(ξL ℓ ). (iii) If ξ n,ℓ → 0, X n,ℓ converges to a Dirac distribution: X n,ℓ d -→ 0.
PROOF. We follow the proof of the univariate case [83, Lemma 2]. Assume that the set {1, . . . , m} decomposes into two disjoint sets C and D such that for indices k ∈ C we have λ n,k → ∞, whereas for an index k ∈ D it holds λ n,k → ρ k for n → ∞. We observe that the mixture of joint raw moments and joint factorial moments converge:

E k∈C X sk n,k λ sk n,k • k∈D X sk n,k → µ s1,...,sm • k∈D ρ sk k .
The latter joint moment sequence, both raw and factorial moments, is exactly the joint moment sequence of a random vector Z = (Z 1 , . . . , Z m ), with

∀k ∈ C : Z k d = L k , ∀k ∈ D : Z k d = MPo(ρ k L k ),
such that

E k∈C Z sk k • k∈D Z sk k = µ s1,...,sm • k∈D ρ sk k .
PROOF OF THEOREM 7.6. The proof is very similar to proofs of Theorems 5.1 and 7.4, so we will be brief again. The mixed factorial moments of X n,jℓ are obtained by differentiation and extraction of coefficients:

E X s1 n,1,j1 • • • X sm n,m,jm = [z n ]∂ s1 v1 . . . ∂ sm vm (F )(z, 1, . . . , 1) [z n ]F (z, 1, . . . , 1)
.

The differentiation with respect to v ℓ only affects the factor

G ℓ H ℓ (z) -(1 -v ℓ )h ℓ,jℓ z jℓ ,
leading to singular expansions covered in Section 3.1 and additional factors h sℓ ℓ,jℓ , 1 ≤ ℓ ≤ m. The asymptotics of these factors for j ℓ → ∞ are all governed by [START_REF] Flajolet | Some exactly solvable models of urn process theory[END_REF]. The individual singular expansions are similar to [START_REF] Flajolet | Analytic urns[END_REF]. Thus, extracting the coefficient z

n-s1j1-•••-smjm from M (z) • m ℓ=1 G (sℓ) ℓ H ℓ (z)
leads to the asymptotic expansion of the joint factorial moments. Lemma 7.7 then yields the stated limit law.

We end this section with four examples covered by the multivariate scheme. EXAMPLE 7.8 (m-bundled plane-oriented recursive trees). We have previously encountered 3-bundled trees in Section 6.2 in the framework of bilabelled trees. In the following we study ordinary increasing trees [START_REF] Bergeron | Varieties of increasing trees[END_REF][START_REF] Drmota | Random Trees[END_REF][START_REF] Janson | Generalized Stirling permutations, families of increasing trees and urn models[END_REF][START_REF] Kuba | On the degree distribution of the nodes in increasing trees[END_REF][START_REF] Hosam | On the structure of random plane-oriented recursive trees and their branches[END_REF][START_REF] Panholzer | Level of nodes in increasing trees revisited[END_REF], where each node has only one label. As before, given a degree-weight sequence (φ j ) j≥0 , the corresponding degree-weight generating function is defined by φ(t) = j≥0 φ j t j . The associated family T of increasing trees can be described by the following symbolic equation using the boxed product (see Section 6.2): T = Z □ * φ T .

For the exponential generating function

T (z) = n≥1 T n z n n! ,
we thus have the following equation T ′ (z) = φ(T (z)), T (0) = 0.

We are interested in families of generalized plane-oriented recursive trees with degree-weight generating function φ(t) = 1/(1 -t) m , with m ∈ N.

For m = 1 we get the ordinary plane-oriented recursive tree, whereas for m > 1 we get the so-called m-bundled trees, with generating function T (z) = 1 -1 -(m + 1)z 1/(m+1) .

One may think of each node holding m -1 additional separation bars [START_REF] Janson | Generalized Stirling permutations, families of increasing trees and urn models[END_REF], which can be regarded as a special edge type. Naturally, this refines the root degree X n in a random tree of size n, since we may look at the number of nodes attached to the root in a specific cluster, induced by the m -1 bars: X n = m ℓ=1 X n,ℓ . By construction, the random variables X n,ℓ are exchangeable, but not independent. Writing u Xn = u 1

Xn,1 • • • u m Xn,m , the generating function

T (z, u) = n≥1
T n E(u Xn ) z n n! of the random vector X n = (X n,1 , . . . , X n,m ) is given by

∂ ∂z T (z, u) = m ℓ=1 1 1 -u ℓ T (z)
.

This refinement is covered by the multivariate scheme (52) (with a shift: X n+1,ℓ , instead of X n,ℓ ). Similarly, one may study the outdegree of a node labelled j, as well as multiple nodes, leading to closely related generating functions [START_REF] Fitzgerald | Generalized Mittag Leffler distributions arising as limits in preferential attachment models[END_REF][START_REF] Kuba | On the degree distribution of the nodes in increasing trees[END_REF][START_REF] Tamás | The maximum degree of the Barabási-Albert random tree[END_REF]. ■ EXAMPLE 7.9 (Bilabelled increasing trees and refined root degree). Continuing from Section 6.2, we can refine the root-degree X n in 3-bundled bilabelled increasing trees of size 2n, X n = X n,1 + X n,2 + X n,3 , where the X n,ℓ are exchangeable. The corresponding generating function T (z, u 1 , u 2 , u 3 ) = Except for the non-standard shift in the random variable, the problem corresponds directly to a multivariate pure critical composition scheme. ■ EXAMPLE 7.10 (Returns in coloured bridges and walks). Consider k-coloured bridges B k as defined in Section 6.4: A bridge is coloured in exactly k colours, where each colour consists of a non-empty bridge. Combinatorially, we append k non-empty bridges one after the other. We are interested in the individual number of returns to zero in each of the first k 1 bridges, followed by k 2 additional bridges, such that k 1 + k 2 = k. Reusing the combinatorial construction [START_REF] Giménez | Asymptotic enumeration and limit laws of planar graphs[END_REF], the multivariate generating function B k (z, u) satisfies

B k (z, u 1 , . . . , u k1 ) =   k1 j=1   1 1 -u j 1 -1 B(z) -1     (B(z) -1) k2 .
By our previous reasoning we see that the corresponding random variables are exchangeable and the multivariate scheme directly applies. Moreover, we can also consider walks and a refined generating function W k (z, u) of k-coloured walks with the tail coloured in the same colour as the final bridge, keeping track of the individual returns to zero of the first k 1 bridges,

W k (z, u) = (1 + B k (z, u)) W (z) B(z) .
Again, the corresponding random variables are exchangeable and the multivariate scheme directly applies again. ■ EXAMPLE 7.11 (Triangular urn models and node degrees in generalized recursive trees). We discuss a specific balanced triangular urn model with k + 1 colours, which generalizes the case with 2 colours from Section 6.7. Our urn model is specified by the following (k + 1) × (k + 1) balanced ball replacement matrix M with α r , β r ∈ N and α r + β r = σ, 1 ≤ r ≤ k + 1: We assume that there are initially A 0,r = a r balls of type r for 1 ≤ r ≤ k + 1 and that the random variables A n,r count the number of balls of type r after n draws. As before, we use the history-counting approach [START_REF] Flajolet | Some exactly solvable models of urn process theory[END_REF][START_REF] Flajolet | Analytic urns[END_REF] to analyse the urn model. For this purpose we will first derive the history generating function H(x 1 , . . . , x k , x k+1 ; z). Due to the balance condition, it suffices to study the function H(x 1 , . . . , x k , 1; z). In the associated differential system, the functions x ℓ (t) satisfy ẋr (z) = x αr+1 r x βr k+1 (z) (for 1 ≤ r ≤ k) and ẋk+1 (z) = x σ+1 k+1 (z), with initial conditions x ℓ (0) = x 0,ℓ . Now, by separation of variables, we directly obtain 14 x k+1 (z) = x 0,k+1 • 1 -σx σ 0,k+1 z -1/σ .

M =       
By integration, we readily get the other functions x r (z) for 1 ≤ r ≤ k:

x r (z) = x 0,r • 1 -x αr 0,r x -αr 0,k+1 1 -(1 -σx σ 0,k+1 z) αr/σ -1/αr .

Finally, using the basic isomorphism between differential systems and history generating functions [START_REF] Flajolet | Some exactly solvable models of urn process theory[END_REF], we obtain the desired closed form for the generating function of urn histories, which falls into the class of our multivariate pure critical composition scheme Theorem 7.4.

PROPOSITION 7.12. The history generating function H(x 1 , . . . , x k , 1; z), associated with the balanced triangular urn model with ball replacement matrix M in (57) and initial conditions A 0,r = a r , 1 ≤ r ≤ k + 1, is given by H(x 1 , . . . , x k , 1; z) = 1 -σz It is well known that node degrees in generalized plane-oriented recursive trees can be modelled by such urns [START_REF] Janson | Asymptotic degree distribution in random recursive trees[END_REF][START_REF] Tamás | The maximum degree of the Barabási-Albert random tree[END_REF]. As an application, we set α r = 1 and β r = k for 1 ≤ r ≤ k + 1 and a r = 1, 1 ≤ r ≤ k and a k+1 = 0 for the initial values. Then, the random variables A n,1 up to A n,k are exchangeable and count the refined root degree of k-bundled plane-oriented recursive trees, i.e., the number of children of the root in each bundle.

Moreover, the urn model M can be used to study the joint degree distribution of the nodes labelled 1, 2, . . . , k in generalized plane-oriented recursive trees; this will be discussed in a forthcoming paper [START_REF] Kuba | A note on central limit theorem analogues for node degree preferential attachment trees and Pitman's α-diversity[END_REF]. ■

14 There is a small sign typo in [37, page 36] where, in the corresponding result for the 2 × 2 model, y -σ 0 should be replaced by y σ 0 .

8. Conclusion. This work makes explicit the limit laws for combinatorial structures counted by schemes like G(H(z))M (z), where we additionally mark either the number of H-components, or the number of H-components of a given size. We focused on the technically more delicate and mathematically richer case where G, H, and M are simultaneously singular. We proved that when these functions have algebraic dominant singularities with exponents between 0 and 1, the limit laws of the schemes are moment-tilted stable distributions and product distributions involving Mittag-Leffler laws. In Table 4 we give an overview of the covered schemes and the corresponding limit laws, where we have convergence in distribution and convergence of all moments.

In the extended scheme, in which one follows the total number of H-components, we proved the appearance of three régimes for the corresponding limit law (continuous, Boltzmann, and linear combination of a continuous and a Boltzmann distribution), depending on the relation of the singular exponents of G, H, and M ; see Table 1. In the size-refined scheme, in which one follows the number of H-components of a given size, we proved the appearance of mixed Poisson distributions (see Definition 3.15) and a double phase transition for the limit law from continuous to discrete to degenerate, with explicit threshold sizes depending on the exponents; see Table 2.

We also presented several extensions (logarithmic singularities, multivariate cases) and a variety of applications to important probabilistic and combinatorial objects. This allowed us to obtain new results for the core size of supertrees, the number of subtrees in different increasing trees, the returns to zero and sign changes in walks and bridges, the table sizes in the Chinese restaurant process, and the number of balls in some urn models. The results of Table 4 hold for functions H(z) having a dominant singularity of algebraic type. Our methods can also deal with schemes having other types of singularities. For example, one could allow algebraic-logarithmic behaviours of the form

H(z) ∼ 1 - z ρ H λH 1 z log 1 - z ρ H ψH .
This would cover some instances of 3-colour balanced triangular urn models [START_REF] Flajolet | Some exactly solvable models of urn process theory[END_REF][START_REF] Puyhaubert | Modèles d'urnes et phénomènes de seuils en combinatoire analytique[END_REF], whose complete analysis, however, remains a challenge as other types of singularities appear, thus leading to new families of limit laws whose nature is unclear. Compare also with the related open problem by Janson [START_REF] Janson | Moments of Gamma type and the Brownian supremum process area[END_REF], asking for a more detailed description of the limit law of unbalanced 2-colour triangular urn models.

To keep this article readable, we eluded the question of the speed of convergence of X n , properly normalized, to its limit law X. Recently, a few articles addressed this for some preferential attachment rules [START_REF] Erol | Generalized gamma approximation with rates for urns, walks and trees[END_REF] and for the Chinese restaurant process [START_REF] Dolera | A Berry-Esseen theorem for Pitman's α-diversity[END_REF] (see Remark 6.16). In fact, following the Puiseux expansions at order 2 of G/H/M allow us to capture this speed of convergence, but this leads to many subcases. In a future work we plan to give uniform bounds on the moments, leading to Berry-Esseen-like inequalities for the speed of convergence of all the limit laws associated with critical composition schemes; see also [START_REF] Adell | On the uniform convergence of normalized Poisson mixtures to their mixing distribution[END_REF][START_REF] Hwang | Second phase changes in random m-ary search trees and generalized quicksort: convergence rates[END_REF].

Last but not least, in our companion article [START_REF] Banderier | Phase transitions of composition schemes: stable laws and Boltzmann distributions[END_REF], we further extend our analysis to schemes with algebraic singularities, where λ H > 1: We analyse a generalization of the composition scheme analysed in [START_REF] Banderier | Random maps, coalescing saddles, singularity analysis, and Airy phenomena[END_REF], leading to stable laws (e.g., the map-Airy distribution), Gaussian laws, as well as bimodal distributions. In extended composition schemes we observe a new behaviour, obtaining for example beta limit laws. In size-refined schemes we anticipate further continuous to discrete phase transitions. This will achieve the exhaustive exploration of the landscape of critical composition schemes with algebraic singularities and associated phase transitions.

EXAMPLE 3 . 17 (

 317 Mixed Poisson half-normal distribution). A half-normally distributed random variable X d = HN(σ) with parameter σ is the absolute value of a normally distributed random variable Y d = N (0, σ 2 ): X d = |Y |. Consequently, X has the probability density function

has a mixed Poisson distribution: Y d =

 d MPo(ξX) with X d = HN(σ). Note that we can readily expand the exponential function and obtain various series representations of P{Y = ℓ}. ■ EXAMPLE 3.18 (Mixed Poisson Rayleigh distribution). A Rayleigh distributed random variable X d = Rayleigh(σ) with parameter σ has the probability density function

  FIG 3. Boltzmann distributions B G (x) have different shapes 8 , related to the singular exponent λ G of the function G: It goes from an asymptotic Gaussian shape if λ G ≪ -1 (and for entire functions) to a spread shape if λ G ≫ -1. The above distributions are drawn with a value of the parameter x such that E(B G (x)) = xG ′ (x)/G(x) = n (with n = 200). It is interesting that they appear in the context of critical schemes, independently of uniform random generation motivations, but, as explained above, both aspects are in fact intimately related.

FIG 4 .

 4 FIG 4. The integration contour of the Cauchy integral (left) is transformed into a larger contour (in red in the middle image) inside a Delta-domain (in blue), becoming in the limit the Hankel contour C (right), allowing the identification of the gamma function.

=

  the random variable Xn ξn converges in distribution, with convergence of all moments, to X.(ii) If ξ n → ξ ∈ (0, ∞),the random variable X n converges in distribution, with convergence of all moments, to a mixed Poisson distributed random variable Y d MPo(ξX). (iii) If ξ n → 0, the random variable X n converges to a Dirac distribution: X n d -→ 0. REMARK 5.5. The second and third case could be grouped together, since for ξ = 0 we have Y d = MPo(0) d = 0. Furthermore, in the third case, for positive random variables (X n ) n∈N the assumptions can be relaxed to simply E(X n ) → 0, without requiring the specific structure of the moments. Note further that the discrete random variable Y d = MPo(ξX) converges, after scaling, to its mixing distribution X: One has Y ξ d -→ ξ→∞ X, with convergence of all moments.

11 FIG 7 .

 117 FIG 7. An increasing diamond and the bijectively equivalent bilabelled increasing plane-oriented recursive tree.

FIG 10 .

 10 FIG 10. A 4-coloured bridge, with all its initial returns to zero marked by red dots.

FIG 13 .

 13 FIG 13. Schematic representation of the tables in the Chinese restaurant model: When the number of customers goes to infinity, one observes, with probability one, many tables with a small number of customers (coloured teal), only a few larger tables of size Θ(n

FIG 14 .

 14 FIG 14. The evolution of a balanced triangular urn with replacement matrix M = 1 1 0 2 and initially w 0 = 2 white balls and b 0 = 1 black ball. Dashed arrows indicate that a white ball was drawn, solid arrows indicate that a black ball was drawn. The number below each node corresponds to the number of possible transitions (or histories) to reach this state.

∂z 2 T 3 ℓ=1 1 1 -

 2311 (z, u 1 , u 2 , u 3 ) = u ℓ (1 -√ 1 -z 2 ).

1 -

 1 x αr 1 -(1 -σz) αr/σ -ar αr .

G 6 TABLE 4

 64 , u) = M (z)G uH(z) , u) = -log 1 -uH(z) F (z, v) = M (z)G H(z) -z j h j (1 -v) v) = -log 1 -H(z) -(1 -v)h j z j /j! ℓ H ℓ (z) -z j ℓ h ℓ,j ℓ (1-v ℓ ) mv. mixed Poisson type phase transition 7.Overview of our results on critical composition schemes, where u = (u 1 , . . . , um) and v = (v 1 , . . . , vm).

  are the same as in Remark 4.2, and where Be(p) (with p

The word "core" comes from the theory of graphs and maps, where this composition scheme is natural and was, e.g., analysed in[START_REF] Banderier | Random maps, coalescing saddles, singularity analysis, and Airy phenomena[END_REF][START_REF] Janson | The birth of the giant component[END_REF].PHASE TRANSITIONS OF COMPOSITION SCHEMES

Note that we still use the letter F to denote the main function. The auxiliary variable is now v and no more u as we are marking a different parameter. This choice avoids more cumbersome notation and is also motivated by the fact that we always have F (z, 1) = F (z).

Flajolet and Sedgewick call the corresponding distribution the profile of the combinatorial object, by analogy with the profile of integer compositions; see[40, p. 169, p.

451, p. 632].

Throughout this article, we use that, by analytical continuation, 1/Γ(m) = 0 whenever m is an integer ≤ 0.

In the literature, there are unfortunately two distinct distributions which are called Mittag-Leffler distribution.

See https://lipn.fr/~cb/Papers/CriticalSchemes/ for several animations of the different limit laws occurring in this article.

See[START_REF] Flajolet | Analytic Combinatorics[END_REF] Appendix B.7] for more on this method; see also[START_REF] Janson | Moments of Gamma type and the Brownian supremum process area[END_REF] Theorem 5.4] with γ = γ ′ = 1 -λ H > 0,and [66, Equation (6.12)] for similar results on the class of functions with moments of Gamma type.

OEIS stands for the On-Line Encyclopedia of Integer Sequences, accessible via https://oeis.org.

Acknowledgments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Acknowledgments.

We are pleased to dedicate this article to our colleague/mentor Alois Panholzer, on the occasion of his 50th birthday. In fact, Alois played a central rôle in the birth of this article in several ways:

• As a byproduct of Alois' analyses of combinatorial structures and stochastic processes (such as permutations, parking functions, lattice paths, urns, and trees), several phase transitions involving the mixed Poisson distribution were uncovered [83,102]. This guided the second author to study a general framework for these phase transitions, and to join forces with the other authors to extend the results of [6,118] involving Rayleigh, half-normal, and stable distributions. This led us to the composition schemes analysed in this work. • Alois also connected the first and second authors via a French-Austrian international project (PHC Amadeus) in 2005/06, perpetuating a long tradition of exchanges in combinatorics between Austria and France. A part of this project led to the article [8], which established Gaussian limit laws for some composition schemes related to trees. • Alois was also a constant source of inspiration and guidance for the second author. Without him, we would have never started this work.

Funding. Michael Wallner was founded by the FWF Erwin Schrödinger Fellowship J 4162 and the FWF Stand-alone Project P 34142.

A Motzkin walk (i.e., step set S = {-1, 0, 1}) with 4 sign changes marked in red.

Combinatorially, we see that the bivariate generating function of bridges is

, where S(z) = 1 1 -p 0 z and H(z) = E(z) S(z) -1.

Here, S(z) counts sequences of horizontal steps, E(z) counts excursions (bridges constrained to be nonnegative; see [START_REF] Banderier | Basic analytic combinatorics of directed lattice paths[END_REF]), and H(z) counts excursions which start with an up or a down step (and not with a horizontal step). We now give the main corresponding Puiseux expansions. First one has

Then, as the radius of convergence 1/p 0 of S(z) is strictly larger than 1/P (1), which is the one of H(z), we see that the additive term S(z) is negligible for the limit law. Thus, we have a composition scheme [START_REF] Aldous | Exchangeability and related topics[END_REF] where M (z) = 2S(z)H(z) has the asymptotic expansion

.

Hence, we have λ - M = 0, which means that the factor M (z) is asymptotically negligible for the law. The asymptotic dominant part arises from P (1) . A similar reasoning (and Remark 4.2) allows us to prove that the number of sign changes in walks asymptotically follows a half-normal distribution with the same parameter σ. We now refine the analysis by counting sign changes which are j steps apart. Then we can apply Theorem 5.1 to get the following refined result which strongly depends on H(z) = j≥0 h j z j . Note that a statement analogous to Remark 6.6 also applies here. COROLLARY 6.9. For walks of length n of Motzkin paths, let the random variable X n,j be the number of sign changes at distance j from the previous sign change or the origin. For walks (resp. bridges) with zero drift (i.e., P ′ (1) = 0), X n,j has factorial moments of mixed Poisson half-normal type (resp. mixed Poisson Rayleigh type)

hj P (1) j • n 1/2 and mixing distributions

for walks, Rayleigh(σ) for bridges, σ = 1 2

.

Furthermore, the random variable X n,j (for walks and for bridges) possesses the three distinct asymptotic régimes of Theorem 5.1, with a phase transition at j = Θ(n 1/3 ).

7. Further extensions.

7.1. Critical cycle scheme. Many combinatorial structures are cycles of more basic building blocks (e.g., cyclic permutations or functional applications are cycles of Cayley trees). If one marks the number of such basic building blocks, this corresponds to

where G = CYC denotes the cycle operator. This scheme is analysed in Flajolet and Sedgewick's magnum opus [40, page 414] in the supercritical case, and we now extend this analysis to the critical case (i.e., by Definition 1.1 for H(ρ H ) = 1). Note that the previous sections were assuming Puiseux-like expansions for the generating function F (z, 1) at its dominant singularity z = ρ = ρ H . Now, for critical cycle schemes, F does not have a Puiseux expansion, so the previous results need to be adapted.

Let us begin with an example: For

which starts with 1, 2, 8, 48, 384, 3840, . . . , and constitutes the entry A000165 in the OEIS.

Here, the moments E(X s n ) are of order n s/2 , so the scaling with 1/ √ n leads directly to moment convergence. This is just one instance of the following more general result. THEOREM 7.1 (Critical schemes with a log). In a critical cycle composition scheme

if H(z) has a singular exponent 0 < λ H < 1, the core size X n (i.e., the number of Hcomponents in structures of size n) has factorial moments given by

.

The scaled random variable X n /(κn λH ) converges in distribution with convergence of all moments to a Mittag-Leffler distributed random variable

REMARK 7.2. Observe that this scheme leads to a distribution similar (except for a shift in the moments) to the one obtained for the scheme involving the sequence operator G = SEQ, i.e., G(z) = 1 1-z , for which one has c G = 1, λ G = -1, and ρ G = 1. Alternatively, we may think of this cycle scheme as the limit case of Theorem 4.1 when λ G → 0. Indeed, for λ - M = 0 the moments (21) of the extended composition scheme can be rewritten into

.

Thus, for λ G → 0 the moments of the random variable X converge to the moments of an ordinary Mittag-Leffler distribution; see Definition 3.9. Similarly, taking the limit λ G → 0 in Remark 4.2 gives λ G = 0 as the tilting parameter, resulting again in the ordinary Mittag-Leffler distribution.

Now, for j ∈ N, we can also look at the size-refined scheme

for which we get the following theorem.

THEOREM 7.3 (Size-refined critical schemes with a log). In the size-refined critical cycle composition scheme

if H(z) has a singular exponent 0 < λ H < 1, then the number X n,j of H-components of size j in structures of size n has factorial moments of mixed Poisson type,

with ξ n,j = ρ j H -cH h j n λH and Mittag-Leffler mixing distribution X d = M λH . The random variable X n,j converges to one of the three limit laws given in Theorem 5.1, depending on whether j = j(n) is smaller, equal, or larger than the critical growth range j = Θ(n

PROOFS OF THEOREM 7.1 AND 7.3. The proofs are analogous to those of Theorems 4.1 and 5.1, and we only point out the differences next. Let us start with the proof of Theorem 7.1. The factorial moments of order s satisfy

,

where F is defined by Equation ( 50). Hence, we first compute

Using the transfer theorems of [START_REF] Flajolet | Analytic Combinatorics[END_REF] we directly obtain

Note that the log function can be replaced by a quasi-inverse using

Thus, the sth factorial moment is obtained from the asymptotics in ( 27) computed for s -1 and with G(z) = 1 1-z . Then, we obtain the final result: the normalized moments converge to the moments of a Mittag-Leffler distribution.

For the proof of Theorem 7.

We measure the size of the G ℓ component by the variable u ℓ ; accordingly this gives

The random vector X n = (X n,1 , . . . , X n,m ) measures the sizes of the G ℓ -components,

. . u km m ]F (z, u 1 , . . . , u m ) [z n ]F (z, 1, . . . , 1)

.

We now introduce a suitable extension of the terms critical (Definition 1.1) and pure (Definition 3.3) for multivariate schemes. We call a multivariate scheme critical if all functions H ℓ (z) have the identical radius of convergence ρ Hℓ = ρ H such that τ ℓ := H ℓ (ρ H ) = ρ Gℓ and M (z) has radius of convergence ρ M ≥ ρ H . We call a multivariate scheme pure if

• H ℓ (z) has a singular exponent 0 < λ Hℓ < 1 for 1 ≤ ℓ ≤ m;

• G ℓ (z) has a singular exponent λ Gℓ < 0 for 1 ≤ ℓ ≤ m;

• M (z) has a singular exponent λ M ≤ 0 or M (z) is analytic at ρ H .

We can now state our multivariate result.