Periodic Pólya urns and an application to Young tableaux

Pólya urns are urns where at each unit of time a ball is drawn and is replaced with some other balls according to its colour. We introduce a more general model: The replacement rule depends on the colour of the drawn ball and the value of the time (mod p ). We discuss some intriguing properties of the diﬀerential operators associated to the generating functions encoding the evolution of these urns. The initial partial diﬀerential equation indeed leads to ordinary linear diﬀerential equations and we prove that the moment generating functions are D-ﬁnite. For a subclass, we exhibit a closed form for the corresponding generating functions (giving the exact state of the urns at time n ). When the time goes to inﬁnity, we show that these periodic Pólya urns follow a rich variety of behaviours: their asymptotic ﬂuctuations are described by a family of distributions, the generalized Gamma distributions, which can also be seen as powers of Gamma distributions. En passant, we establish some enumerative links with other combinatorial objects, and we give an application for a new result on the asymptotics of Young tableaux: This approach allows us to prove that the law of the lower right corner in a triangular Young tableau follows asymptotically a product of generalized Gamma distributions.


Periodic Pólya urns
Pólya urns were introduced in a simplified version by George Pólya and his PhD student Florian Eggenberger in [7,8,27], with applications to disease spreading and conflagrations.
They constitute a powerful model, still widely used: see e.g.Rivest's recent work on auditing elections [28], or the analysis of deanonymization in Bitcoin's peer-to-peer network [9].They are well-studied objects in combinatorial and probabilistic literature [2,11,22], and offer fascinatingly rich links with numerous objects like random recursive trees, m-ary search trees, branching random walks (see e.g.[3,6,15,16,30]).In this paper we introduce a variation which offers new links with another important combinatorial structure: Young tableaux.We solve the enumeration problem of this new model, derive the limit law for the evolution of the urn, and give some applications.
In the Pólya urn model, one starts with an urn with b 0 black balls and w 0 white balls at time 0. At every discrete time step one ball is drawn uniformly at random.After inspecting its colour it is returned to the urn.If the ball is black, a black balls and b white balls are added; if the ball is white, c black balls and d white balls are added (where a, b, c, d ∈ N are non-negative integers).This process can be described by the so-called replacement matrix: We call an urn and its associated replacement matrix balanced if K := a + b = c + d.In other words, in every step the same number K of balls is added to the urn.This results in a deterministic number of balls after n steps: b 0 + w 0 + Kn balls.Now, we introduce a more general model which has rich combinatorial, probabilistic, and analytic properties.Definition 1.A periodic Pólya urn of period p with replacement matrices M 1 , M 2 , . . ., M p is a variant of a Pólya urn in which the replacement matrix M k is used at steps np + k.Such a model is called balanced if each of its replacement matrices is balanced.
In this article, we illustrate the aforementioned rich properties on the following model (the results for other values of the parameters are similar to the case we now handle in detail).Definition 2. We call a Young-Pólya urn the periodic Pólya urn of period 2 with replacement matrices M 1 := 1 0 0 1 for every odd step, and M 2 := 1 1 0 2 for every even step.
Let us describe the state of the urn after n steps by pairs (number of black balls, number of white balls), starting with b 0 = 1 black ball and w 0 = 1 white ball shown in Figure 1.In the first step the matrix M 1 is used and gives the two states (2, 1), and (1, 2).In the second step, matrix M 2 is used, in the third step, matrix M 1 is used again, in the fourth step, matrix M 2 , etc. Thus, the possible states are (3, 2), (2, 3), and (1, 4), at time 2, and (4, 2), (3,3), (2,4), and (1, 5), at time 3.
In fact, each of these states may be reached in different ways, and such a sequence of transitions is called a history.Each history comes with weight one.Implicitly, they induce a probability measure on the states at step n.So, let B n and W n be random variables for the number of black and white balls after n steps, respectively.As our model is balanced, B n + W n is a deterministic process, reflecting the identity B n + W n = b 0 + w 0 + n + n 2 .So, from now on, we concentrate our analysis on B n .

11:3
Figure 1 The evolution of a Young-Pólya urn with one initial black and one initial white ball.Black arrows mark that a black ball was drawn, dashed arrows mark that a white ball was drawn.Straight arrows indicate that the replacement matrix M1 was used, curly arrows show that the replacement matrix M2 was used.The number below each node is the number of possible transitions to reach such a state.In this article we give a formula for Hn (which encodes all the possible states of the urn at time n) and their asymptotic behaviour.
For the classical model of a single balanced Pólya urn, the limit law of the random variable B n is fully known: The possible limit laws include a rich variety of distributions.To name a few, let us mention the uniform distribution [10], the normal distribution [3], and the Beta and Mittag-Leffler distributions [15].Periodic Pólya urns (which include the classical model) lead to an even larger variety of distributions involving a product of generalized Gamma distributions [31].Definition 3. The generalized Gamma distribution GenGamma(α, β) with real parameters α, β > 0 is defined by the density function (having support (0, +∞)) where Γ is the classical Gamma function Γ(z) := ∞ 0 t z−1 exp(−t) dt.Remark.Let Γ(α) be the Gamma distribution 1 of parameter α > 0, given by its density Then, one has Γ(α) L = GenGamma(α, 1) and, for r > 0, the distribution of the r-th power of a random variable distributed according to Γ(α) is Γ(α) r L = GenGamma(α/r, 1/r).
Our main results are the enumeration result from Theorem 5, the application to Young tableaux in Theorem 7, and the following result (and its generalization in Theorem 6): Theorem 4. The normalized random variable 2 2/3 3 Bn n 2/3 of the number of black balls in a Young-Pólya urn converges in law to a generalized Gamma distribution: 1 Caveat: It is traditional to use the same letter for both the Γ function and the Γ distribution.Also, some authors add a second parameter to the distribution Γ, which is set to 1 here.

Periodic Pólya urns and an application to Young tableaux
We give a proof of this result in Section 3. Let us first mention some articles where this distribution has already appeared before: in Janson [17], for the analysis of the area of the supremum process of the Brownian motion, in Peköz, Röllin, and Ross [25], as distributions of processes on walks, trees, urns, and preferential attachments in graphs (they also consider what they call a Pólya urn with immigration, which is a special case of a periodic Pólya urn), in Khodabin and Ahmadabadi [19] following a tradition to generalize special functions by adding parameters in order to capture several probability distributions, such as e.g. the normal, Rayleigh, and half-normal distribution, as well as the MeijerG function (see also the addendum of [17], mentioning a dozen of other generalizations of special functions).
In the next section we translate the evolution process into the language of generating functions by encoding the dynamics of this process into partial differential equations.

A functional equation for periodic Pólya urns
Let h n,k, be the number of histories of a periodic Pólya urn after n steps with k black balls and white balls, with an initial state of b 0 black balls and w 0 white balls, and with replacement matrices M 1 for the odd steps and M 2 for the even steps.We define the polynomials Note that these are indeed polynomials as there are just a finite number of histories after n steps.We collect all these histories in the trivariate exponential generating function In particular, we get for the first 3 terms of H(x, y, z) the expansion (compare Figure 1) Observe that the polynomials H n (x, y) are homogeneous, as we have a balanced urn model.Now it is our goal to derive a partial differential equation describing the evolution of the periodic Pólya urn model.For a comprehensive introduction to the method we refer to [10].
In order to capture the periodic behaviour we split the generating function H(x, y, z) into odd and even steps.We define Next, we associate to the replacement matrices M 1 and M 2 from Definition 2 the differential operators D 1 and D 2 , respectively.We get where ∂ x and ∂ y are defined as the partial derivatives ∂ ∂x and ∂ ∂y , respectively.These model the evolution of the urn.For example, in the term x 2 ∂ x , the derivative ∂ x represents drawing a black ball and the multiplication by x 2 returning the black ball and an additional black ball into the urn.The other terms have analogous interpretations.
With these operators we are able to link odd and even steps with the following system Note that the derivative ∂ z models the evolution in time.This system of partial differential equations naturally corresponds to recurrences at the level of coefficients h n,k, , and vice versa.This philosophy is well explained in the symbolic method part of [12] and see also [10].
As a next step we want to eliminate the y variable in these equations.This is possible as the number of balls in each round and the number of black and white balls are connected due to the fact that we are dealing with balanced urns.First, as observed previously, one has number of balls after Therefore, for any x k y z n appearing in H(x, y, z) with b 0 = w 0 = 1 we have This translates directly into ( Finally, combining (1) and ( 3), we eliminate ∂ y H e and ∂ y H o .After that it is legitimate to insert y = 1 as there appears no differentiation with respect to y anymore.As the urns are balanced, the exponents of y and x in each monomial are bound (see Equation ( 2)), so we are losing no information on the trivariate generating functions by setting y = 1.Hence, from now on we use the notation H(x, z), H e (x, z), and H o (x, z) instead of H(x, 1, z), H e (x, 1, z), and H o (x, 1, z), respectively.All of this leads to our first main enumeration theorem: Theorem 5 (Linear differential equations and hypergeometric expressions for histories).The generating functions describing the 2-periodic Young-Pólya urn at even and odd time satisfy the following system of differential equations: if n is even, Alternatively, this sequence satisfies h(n + 2) = 2  3 h(n + 1) + 1 4 (9 n 2 + 21 n + 12)h(n).This sequence is not found in the OEIS 2 , we added it there, it is now A293653, and it starts like this: 1, 2, 6, 30, 180, 1440, 12960, 142560, 1710720, 23950080, 359251200, . . . 2 On-Line Encyclopedia of Integer Sequences, https://oeis.org.

Periodic Pólya urns and an application to Young tableaux
In the next section we will use Equations ( 4) to iteratively derive the moments of the distribution of black balls after n steps.

Moments of periodic Pólya urns
In this section, we give a proof via the method of moments of Theorem 4 stated in the introduction.Let m r (n) be the r-th factorial moment of the distribution of black balls after n steps, i.e.
Expressing them in terms of the generating function H(x, z), it holds that Splitting them into odd and even moments, we have access to these quantities via the partial differential equation ( 4).As a first step we compute h n := [z n ]H(1, z), the total number of histories after n steps.We substitute x = 1, which makes the equation independent of the derivative with respect to x.Then, the idea is to transform (4) into two independent differential equations for H e (1, z) and H o (1, z).This is achieved by differentiating the equations with respect to z and substituting the other one to eliminate H e (1, z) or H o (1, z), respectively.This decouples the system, but increases the degree of differentiation by 1.We get In this case it is easy to extract the underlying recurrence relations and solve them explicitly.This also leads to the closed forms (5) for h n , from which it is easy to compute the asymptotic number of histories for n → ∞.Interestingly, the first two terms in the asymptotic expansion are the same for odd and even number of steps, only the third ones differ.We get As a next step we compute the mean.Therefore, we differentiate (4) once with respect to x, substitute x = 1, decouple the system, derive the recurrence relations of the coefficients, and solve them.Note again that the factor (x − 1) prevents higher derivatives from appearing and is therefore crucial for this method.After normalization by h n we get if n is even, if n is odd.
For the asymptotic mean we discover again the same phenomenon that the first two terms in the asymptotic expansion are equal for odd and even n.Differentiating (4) to higher orders allows to derive higher moments in a mechanical way (this however requires further details, which will be included in the expanded version of this article).In general we get the closed form for the r-th factorial moment

11:7
Therefore we see that the moments E (B * n r ) of the rescaled random variable B * n := 2 2/3 3 Bn n 2/3 converge for n to infinity to the limit Note that one has m 1)) for large r, so the following sum diverges: Therefore, a result by Carleman (see [5, pp. 189-220] or [33, p. 330]) 3 implies that there exists a unique distribution (let us call it D) with such moments m r .Furthermore, by the asymptotic result from Equation ( 6) there exist an n 0 > 0 and constants a r and b r independent of n such that a r < m r (n) < b r , for all n ≥ n 0 .Thus, by the limit theorem of Fréchet and Shohat [13] 4 there exists a limit distribution (which therefore has to be D) to which a subsequence of our rescaled random variables B * n converge to.And as we know via Carleman's criterion above that D is uniquely determined by its moments, it is in fact the full sequence of B * n which converges to D. Now it is easy to check that if X ∼ GenGamma(d, p) is a generalized Gamma distributed random variable (as defined in Definition 3), then it is a distribution determined by its moments, which are given by E(X r ) = Γ d+r p /Γ d p .In conclusion, the structure of m r in Formula ( 7) implies that the normalized random variable B * n of the number of black balls in a Young-Pólya urn converges to GenGamma (1,3) .This completes the proof of Theorem 4.
The same approach allows us to study the distribution of black balls for the urn with 1 0 0 1 and M p = 1 0 1 + .We call this model the Young-Pólya urn of period p and parameter .
Theorem 6.The renormalized distribution of black balls in the Young-Pólya urn of period p and parameter is asymptotically a distribution, which we call ProdGenGamma(p, , b 0 , w 0 ), defined as the following product of independent distributions: with δ = p/(p + ), and where Beta(b 0 , w 0 ) is as usual the law with support [0, 1] and density Sketch.This follows from the following r-th (factorial) moment computation:

11:8
Periodic Pólya urns and an application to Young tableaux which in turn characterizes the ProdGenGamma distribution.Indeed, if for some independent random variables X, Y, Z, one has E(X r ) = E(Y r )E(Z r ) (and if Y and Z are determined by their moments), then This is consistent with our results on the Young-Pólya urn introduced in Section 1. Indeed, there one has w 0 = b 0 = 1, p = 2, = 1, and therefore the renormalized distribution of black balls p δ p+ B n /n δ is asymptotically Unif(0, 1) • GenGamma(4, 3) = GenGamma (1, 3).We will now see what are the implications of this result on an apparently unrelated topic: Young tableaux.

Urns, trees, and Young tableaux
As predicted by Anatoly Vershik in [32], the 21st century should see a lot of challenges and advances on the links of probability theory with (algebraic) combinatorics.A key role is played here by Young tableaux 5 , because of their ubiquity in representation theory.Many results on their asymptotic shape have been collected, but very few results are known on their asymptotic content when the shape is fixed (see e.g. the works by Pittel and Romik, Angel et al., Marchal [1,24,26,29], who have studied the distribution of the values of the cells in random rectangular or staircase Young tableaux, while the case of Young tableaux with a more general shape seems to be very intricate).It is therefore pleasant that our work on periodic Pólya urns allows us to get advances on the case of a triangular shape, with any slope.
For any fixed integers n, , p ≥ 1, we introduce the quantity N := p n(n + 1)/2.We define a triangular Young tableau of slope − /p and of size N as a classical Young tableau with N cells with length n and height np such that the first p rows (from the bottom) have length n , the next p lines have length (n − 1) and so on (see Figure 2).We now study what is the typical value of its lower right corner (with the French convention for drawing Young tableaux, see [21] but take however care that on page 2 therein, Macdonald advises readers preferring the French convention to "read this book upside down in a mirror"!).
It could be expected (e.g. via the Greene-Nijenhuis-Wilf hook walk algorithm for generating Young tableaux, see [14]) that the entries near the hypotenuse should be N − o(N ).Can we expect a more precise description of these o(N ) fluctuations?Our result on periodic urns enables us to exhibit the right critical exponent, and the limit law in the corner: (Recall that ProdGenGamma is defined by Formula 9.) 5 A Young tableau of size n is an array with columns of (weakly) decreasing height, in which each cell is labelled, and where the labels run from 1 to n and are strictly increasing along rows from left to right and columns from bottom to top, see Figure 2. We refer to [21] for a thorough discussion on these objects.

Remark:
The simplest case ( = 1, p = 2) relates to the Young-Pólya urn model which we analysed in the previous sections.

Sketch of proof.
We first establish a link between Young tableaux and linear extensions of trees.Then we will be able to conclude via a link between these trees and periodic Pólya urns.Let us start with Figure 2, which describes the main characters of this proof.The bottom part of Figure 2 presents two trees (the "big" tree T , which contains the "small" tree S).More precisely, we define the rooted planar tree S as follows: The left-most branch of S has n + 1 vertices, which we call v 1 , v 2 , . . ., v n +1 , where v 1 is the root and v n +1 is the left-most leaf of the tree.For 2 ≤ k ≤ n − 1, the vertex v k has p + 1 children.
The vertex v n has p − 1 children.
All other vertices v j (for j < n , j = k ) have exactly one child.Now, define T as the "big" tree obtained from the "small" tree S by adding a vertex v 0 as the father of v 1 and adding N + 1 − n(p + ) children to v 0 (see Figure 2).Remark that the number of vertices of T is equal to 1 + the number of cells of Y.Moreover, the hook length of each cell in the first row (from the bottom) of Y is equal to the hook length of the corresponding vertex in the left-most branch of S.
Let us now introduce a linear extension E T of T , i.e. a bijection from the set of vertices of T to {0, 1, . . ., N } such that E T (u) < E T (u ) whenever u is an ancestor of u .A key result, which will be proved in the expanded version of this abstract, is the following: if E T is a uniformly random linear extension of T , then X n (the entry of the lower right corner in a uniformly random Young tableau with shape Y) has the same law as E T (v n ): What is more, recall that T was obtained from S by adding a root and some children to this root.Therefore, one can obtain a linear extension of the "big" tree T from a linear extension of the "small" tree S by a simple insertion procedure.This allows us to construct a uniformly random linear extension E T of T and a uniformly random linear extension E S of S such that So, to summarize, we have now The last step (which we just state here, see our long version for its full proof) is that Indeed, more precisely N − E S (v n ) has the same law as the number of black balls in a periodic urn after (n − 1)p steps (an urn with period p, with adding parameter , and with initial conditions w 0 = and b 0 = p).Thus, our results on periodic urns from Section 3 and the conjunction of Equations ( 10), (11), and (12) gives the convergence in law for X n which we wanted to prove.
A o f A 2 0 1 8  2 In this section, we see that there is a relation between Young tableaux with a given periodic shape, some trees, and the periodic Young-Pólya urns.The lower right corner of these Young tableaux is thus following the same generalized Gamma distribution we proved for urns.

Conclusion and further work
In this article, we introduced Pólya urns with periodic replacements, and showed that they can be exactly solved with generating function techniques, and that the initial partial differential equation encoding their dynamics leads to linear (D-finite) moment generating functions, which we identify as a product of generalized Gamma distributions.Note that [20,23] involve the asymptotics of a related process (by grouping p units of time at once of our periodic Pólya urns).This related process is therefore "smoothing" the irregularities created by our periodic model, and allows us to connect with the usual famous key quantities for urns, such as the quotient of eigenvalues of the substitution matrix, etc.Our approach has the advantage to describe each unit of time (and not just what happens after "averaging" p units of time at once), giving more asymptotic terms, and also exact enumeration.
In the full version of this work we will consider arbitrary periodic balanced urn models, and their relationship with Young tableaux.It remains a challenge to understand the asymptotic landscape of Young tableaux, even if it could be globally expected that they behave like a Gaussian free field, like for many other random surfaces [18].As a first step, understanding the fluctuations and the universality of the critical exponent at the corner could help to get a more global picture.Note that our results on the lower right corner directly imply similar results on the upper right corner: just use our formulae by exchanging and p, i.e. for a slope corresponding to the complementary angle to 90 o .Thus the critical exponent for the upper right corner is 2 − δ.In fact, it is a nice surprise that there is even more structure: there is a duality between the limit laws X and X of these two corners and we get the factorization as independent random variables (up to renormalization and slight modifications of the boundary conditions) L = Γ(b 0 ).Similar factorizations of the exponential law, which is a particular case of the Gamma distribution, have appeared recently in relation with functionals of Lévy processes, following [4].

Theorem 7 .
Choose a uniform random triangular Young tableau Y of slope − /p and size N = p n(n + 1)/2 and put δ = p/(p + ).Let X n be the entry of the lower right.Then (N − X n )/n 1+δ converges in law to the same limiting distribution as the number of black balls in the periodic Young-Pólya urn with initial conditions w 0 = , b 0 = p and with replacement matrices M 1 = • • • = M p−1 e. we have the convergence in law, as n goes to infinity: