On Random Allocation Models in the Thermodynamic Limit

We discuss the phase transition and critical exponents in the random allocation model (urn model) for different statistical ensembles. We provide a unified presentation of the statistical properties of the model in the thermodynamic limit, uncover new relationships between the thermodynamic potentials and fill some lacunae in previous results on the singularities of these potentials at the critical point and behaviour in the thermodynamic limit. The presentation is intended to be self-contained, so we carefully derive all formulae step by step throughout. Additionally, we comment on a quasi-probabilistic normalisation of configuration weights which has been considered in some recent studies


I. INTRODUCTION
The random allocation model, also known as the balls-in-boxes model, the urn model or the backgammon model [1][2][3], is a simple statistical model describing weighted random partitions of particles between boxes.Despite its simplicity the model exhibits very rich critical behavior, including discontinuous and continuous phase transitions of different orders depending on the model parameters and the ensemble being considered.The phase transition in the model is related to a real-space condensation observed in many statistical problems including zero-range processes [4][5][6][7][8][9][10][11], mass-transport [12][13][14], random trees [15,16], quantum gravity [3,17,18].The model has been used to understand some aspects of non-equilibrium dynamics of condensate formation [8].The balls-in-boxes model has also been applied in studies of such diverse problems as wealth condensation [19] or diversity of Zipf's populations [20].The model can also be used to mimic phase separation [21,22], condensation in complex networks [23,24], fire-ball formation at the van Hove singularity [25], formation of the giant component/cluster in networks or percolation models and the statistics of the longest interval in tied-down renewal processes [26][27][28].
In the present paper we revisit the issue of the phase transition in the model, and analyse it carefully from the point of view of equilibrium statistical mechanics.We will describe in detail how the order of the transition depends on the parameters of the model in various ensembles and present some new results on the singularities of the thermodynamic potentials, the relationships between them and the (slow) approach to the limiting distribution in the condensed phase.We also discuss some subtleties in a quasi-probabilistic normalisation used in [20] which regularizes an otherwise divergent sum over the weights and elucidate its exact relation to the class of models here.
The paper is organized as follows.In the next section we will introduce a canonical ensemble and the basic quantities that describe the behaviour of the system in this canonical ensemble.We then apply the saddle point method to calculate the free energy density of the system for low particle densities that are smaller than some critical density.For densities larger than this critical density the system exhibits a real-space condensation.Whether there is a phase transition at a finite critical density or not depends on the asymptotic behaviour of the weights governing the particle distribution.We will divide the weight functions into three families for which the system has only the fluid phase, only the condensed phase, or has both phases with a phase transition between them at a finite critical density.We will concentrate on power-law weights to analyse the types of possible phase transitions.We will list possible scenarios that depend on the power in the weight function and determine the free energy and its behaviour at the critical point.We will then repeat these steps for the grand-canonical ensemble.As we will see, the canonical system has a continuous phase transition, while the grand-canonical system may additionally have a discontinuous phase transition.We then move on to discuss other ensembles and the quasi-probabilistic weights of [20].We conclude with a short summary.

II. CANONICAL ENSEMBLE
The balls-in-boxes model is defined by the partition function [1] Z S,N = (s1,...,s N ) w(s 1 ) . . .w(s N )δ S−(s1+...+s N ) (1) that describes weighted distributions of S particles in N boxes, where s i 's denote the occupations of boxes i = 1, . . ., N .The lowercase δ represents the Kronecker delta: δ n = 1 for n = 0 and δ n = 0 for all other integers n ̸ = 0.The delta in Eq. ( 1) selects configurations that have exactly S particles in N boxes.The statistical weight of a configuration (s 1 , s 2 , . . ., s N ), that describes the partition of particles between boxes, is the product of statistical weights w(s i ) of individual boxes, that depend only on the number of particles in the box.The numbers of particles in different boxes are almost independent of each other, but the constraint s 1 +. ..+sN = S makes them weakly dependent.In some conditions, the presence of this constraint leads to a phase transition, as will be seen later.The statistical weights w(s i ) are identical for all boxes.The model is thus symmetric under the permutation of box labels.The phase transition observed in the model is related to spontaneous permutation symmetry breaking, where one box behaves differently to the others and captures an extensive number of particles.
In general, the statistical weight s → w(s) is a non-negative real-valued function defined for s = 0, 1, 2, . ... However, we find it convenient to limit the range of s to positive integers s = 1, 2, . .., in which case each box must contain at least one particle and there must be at least S = N particles in a system with N boxes.The full range (including empty boxes) can be always recovered by introducing new weights w ′ (s) = w(s + 1) for s = 0, 1, 2, . . .and redefining the box occupation numbers s ′ i = s i − 1 and S ′ = S − N .In the sequel we assume that there are no empty boxes.This version of the model occurs naturally in many of the problems mentioned in the introduction.
Let us consider a system with q balls in the first box.Then the rest of the boxes contains S − q balls and are described by the partition function Z S−q,N −1 .This leads to the recurrence relation for S ≥ N ≥ 1, and for N = 0.The formula for Z S,0 ensures that as required.Consider a physical quantity O = O(s 1 , . . ., s N ).The ensemble average is defined as It is worth noting that the following transformation w(s) → e λ+σs w(s). ( leaves the ensemble averages (4) unchanged because it introduces the same factors e λN +σS in the numerator and denominator of (4) and they cancel out.We will use this invariance in the sequel many times.
As an example, consider the box-occupation probability which is the probability that a randomly selected box contains q particles: Because all the boxes are identical, the ensemble average of π(q) is the same as the ensemble average of the probability that the first box has q particles δ q−s1 .This probability can be calculated directly from the partition function by observing that if the first box contains q particles then the remaining N − 1 boxes contain S − q particles.This leads to the relation for S ≥ N ≥ 1 and q = 1, 2, . . ., S − N + 1. Relation (2) ensures that this distribution is normalised.
As an illustration consider the model with weights w(s) = 1 for s = 1, 2, . ... One finds that and for q = 1, 2, . ... Note that the invariance of the ensemble averages (6) under the transformation (5) leads to the somewhat counter-intuitive conclusion that exponentially increasing or decreasing weights w(s) = e λ+σs give exactly the same results as the constant weights w(s) = 1.In particular, the particle distribution for the exponentially increasing or decreasing weights will be given by (11).

III. THERMODYNAMIC LIMIT
We are interested in the behaviour of the system in the thermodynamic limit: where r ∈ (0, 1).The behaviour depends on the form of w(s) but also on r which is a free parameter of the model, being the reciprocal of the average particle density ρ: We will use r and ρ interchangeably.
To describe the behaviour of the system in the thermodynamic limit (12), it is convenient to introduce a thermodynamic potential where "lim" in this equation means the limit (12).With a mild misuse of terminology ϕ(r) may be called the free energy density (free energy per particle).The free energy density, ϕ(r), is the rate of the asymptotic growth of the partition function with the system size for given r Z S,N ∝ e Sϕ(r) (15) in the thermodynamic limit (12).The sub-exponential corrections are omitted.For trivial weights, that is for w(s) = 1 for s = 1, 2, . .., discussed in the previous section, the partition function (10) behaves in this limit asymptotically as so the free energy density is This can be seen by applying Stirling's formula to (10).It is also easy to see that the particle distribution (11) takes the following asymptotic form in the limit ( 12) for q = 1, 2, . ... The asymptotic expressions ( 16) and ( 18) hold for any r ∈ (0, 1).They break down for r = 0, (infinite density) that is when the number of particles S grows faster than linearly with the number of boxes N as N goes to infinity.They also break down for r = 1, that is when the difference S − N grows slower than linearly with N as N goes to infinity.In the general case, the thermodynamic potential ϕ(r) can be calculated using the saddle point method.Using an integral representation of the Kronecker delta we can rewrite the partition function (1) as follows The right hand side of this equation can be more concisely written as where K(α) is a cumulant generating function In the limit (12), the leading contribution to the integral is where and α(r) is a solution of the saddle point equation.
This result is obtained by deforming the integration contour so that it passes through a saddle point.The saddle point is located on the real axis and its position α(r) depends on the particle density.When the particle density ρ = 1/r increases, α(r) decreases.The minimal value that it may take is limited by the radius of convergence of the series ( 22) For α > α c the series in (22) is convergent.The saddle point solution α(r) holds only for r > r c For r → r + c , the saddle point α(r) → α c approaches the singularity of K(α) (22).In the thermodynamic limit (12), we can also derive an asymptotic form of the particle distribution (8) by replacing the partition functions in the numerator and denominator by their asymptotic forms (23): and where r = N/S and r = (N − 1)/(S − q).In the thermodynamic limit ( 12) Using this expansion we obtain log Z S−q,N −1 ∝ (S − q)α(r) The calculations are straightforward.In particular, they show that the O(1) terms cancel out on virtue of the saddle point equation (25).Dropping the O(S −1 ) terms we finally obtain in the thermodynamic limit ( 12) and ⟨π(q)⟩ S,N ∝ w(q)e −qα(r) e K(α(r)) = w(q)e −qα(r) ∞ s=1 w(s)e −sα(r) . ( We see that the distribution is suppressed by an exponential factor ∼ exp (−q(α(r) − α c )) for large q.
The first moment of the particle distribution should be exactly equal to the particle density Replacing ⟨π(q)⟩ S,N on the right hand side of this equation with (33) we find consistently with the saddle point equation (25).
Neglecting the O(S −1 ) terms we have assumed that q ≪ S.This is true for particle densities smaller than where this result is expected to hold.Above ρ c another condensed phase appears -the condensed phase where an extensive (proportional to S) number of particles condense in one box.In the next section we discuss in detail the phase transition between the two phases.

IV. PHASE STRUCTURE
As mentioned, the system can have two distinct phases: a symmetric one (also called the fluid phase) for small particle densities and a phase with a spontaneously broken permutation symmetry (also called the condensed phase) for large particle densities.The fluid phase is described by the solution of the saddle point equations (25).In the condensed phase, the equations are no longer valid.The critical density at which the transition from one phase to another occurs is given by (36) with r c defined by (27).Whether the system has a phase transition or not depends on the weight function w(s).Let us discuss the possible scenarios.
The set of possible weight functions can be divided into three families which are classified by the value of the parameter α c (26): the first family includes weight functions w(s) for which α c = −∞, the second includes weights for which α c = +∞, and the third weights for which α c is a finite number.
• Weight functions w(s) from the first family (α c = −∞) fall off asymptotically to zero faster than exponentially for s → ∞, for instance w(s) = 1/s!, w(s) = e −s 2 .Weights which vanish for all s above a certain value s 0 : w(s) = 0 for all s > s 0 also belong to this category.
• Weight functions w(s) from the second family (α c = +∞) grow faster than exponentially for s → ∞, for instance w(s) = s! or w(s) = e s 2 .
• Weight functions w(s) from the third family (−∞ < α c < +∞) neither increase nor decrease faster than exponentially for s → ∞, for instance w(s) = e σs s −β or w(s) = e ± √ s .One should note that power-like weights, w(s) = s −β , belong to this category, even if β is negative.
For the first family the system has only the fluid phase, for the second one it has only the condensed phase.The third family, in a way, interpolates between the first two cases and therefore it is the most interesting one.
From now on, we focus on weight functions from the third family, that is such that α c (26) is finite.Using the transformation (5), which does not affect the ensemble averages, we can transform the weights w(s) → e −αcs w(s) so that the critical value (26) after the transformation is α c = 0.So from now on, unless we specify otherwise, we set by default α c = 0.
The inverse particle density r is a free parameter that can change in the range r ∈ (0, 1).The saddle point equation ( 25) holds for r ∈ (r c , 1) where r c is given by ( 27) (with α c = 0).Two situations can be distinguished: −K ′ (0) is infinite or finite.In the former case the critical density is infinite, so there is no phase transition and the saddle point solution holds for the whole range of r ∈ (0, 1).In the latter case, the critical density is finite, so there is a phase transition at ρ c = −K ′ (0).The saddle point solution holds for r ∈ (r c , 1).As mentioned, in this case the particle distribution (33) falls off exponentially for large q as π(q) ∝ w(q)e −qα(r) with α(r) > 0. At r = r c the exponential factor disappears, α(r c ) = 0, and the particle distribution (33) tends in the thermodynamic limit (12) to its critical form whose mean is equal to the critical density To see what happens for ρ > ρ c , in the condensed phase, below we will perform a finite size analysis.More precisely, we will numerically determine the particle distribution ⟨π(q)⟩ S,N for finite N using formula (8) and the recursion (2).In fact, instead of (2) we have used the following recursion in the calculations.It is more efficient than (2), because it doubles N in one step, while (2) increases N by one.
As an example, we study the behaviour of the system for power-law weights with β = 4.0 and ρ = 1.4,which is greater than the critical density so the system is in the condensed phase.The results are presented in Fig. 1.As we can see from the first picture the distribution can be divided into a "bulk" part corresponding to the critical distribution denoted by the dashed line in the figure and a peak.The peak is located at Q ≈ N (ρ − ρ c ).The area under the peak tends to 1/N .This picture can be understood as a split of the system into a critical subsystem consisting of N − 1 boxes and a single box that captures the excess particles that did not fit in the critical subsystem.The critical part consists of S N −1 particles in N − 1 boxes.S N −1 can be approximated as a sum of N − 1 independent random numbers q 1 + . . .+ q N −1 distributed according to the critical distribution π c (q).By the generalised central limit theorem [29], such a sum grows on average as (N − 1)ρ c , and behaves in the limit N → ∞ like a normal random variable with the variance (N − 1)K ′′ (0) if K ′′ (0) < ∞, or as a one-sided α-stable random variable with α ∈ (1, 2), otherwise.This tells us that the number of particles in the remaining box is on average Q = S − S N −1 ≈ N (ρ − ρ c ), and its fluctuations about the mean are of order √ N if K ′′ (0) < ∞ or N 1/α otherwise.In the discussed example K ′′ (0) < ∞, so the peak can be approximated by a normalised Gaussian curve with mean Q and variance N K ′′ (0) with an additional 1/N factor that reflects the fact that the condensate is in one of the N boxes.The height of the peak is of order N −3/2 as follows from the last equation.
The approach to the limit distribution (42) is very slow and nonuniform as shown in Fig. 1.The peak which is present for any finite N does not disappear but moves away to infinity when N increases.The shape of the peak deviates from the Gauss-curve shape.The deviations are seen as long tails on the left side of parabolas in a semi-logarithmic plot.Slight deviations of tails on the right hand side can also been seen.The left tails are remnants of the power-law tail q −β of the critical distribution (37) [8].The contribution from the tails decreases as N grows, as discussed below.The excess probability accumulated in the tails decreases as N increases.The excess probability is calculated as the area between the curve ⟨π(q)⟩ S,N and a curve obtained as a weighted sum of the limiting distribution (42) and the Gaussian peak (43), to the left of the maximum of the peak.The weights (N − 1)/N and 1/N correspond to the fractions of boxes occupied by the critical subsystem and the condensate, respectively.For example, for the three graphs shown in the right plot in Fig. 2, which correspond to N = 2 k + 1 with k = 11, 12, 13, the excess probability accumulated in the tail as compared to the probability 1/N accumulated in the peak, P dev /(1/N ), is approximately 4.2%, 2.8% and 1.8%, respectively.This means that extra events associated with this part of the distribution are suppressed in the thermodynamic limit.The events in this part of the distribution correspond to the situation where also the number of particles in the second most occupied box is significantly larger than in the remaining boxes, or in other words that the excess of particles is shared by two boxes [8].
with β, λ, σ being arbitrary real parameters.These weights can be reduced to power-law weights by the transformation (5) w(s) → w(s) = e −λ−σs w(s).As follows from (6) the ensemble averages for the power-law weights (46) are the same as for the weights w(s) (45).We will therefore restrict ourselves to the weights (46), which are representative for the whole family (45).The generating function K(α) is where Li β (z) is the polylogarithm [30]: For non-integer β, the polylogarithm of e −α has the following series expansion at For integer β, the singular term contains a logarithmic singularity where H n = 1+1/2+. ..+1/n is the n-th harmonic number.The behaviour of K(α) for α → 0 + depends on β.Regarding the critical behaviour of the model, which is encoded in the singular behaviour of the free energy density ϕ(r) (14), we can distinguish four cases: The free energy density ϕ(r) can be calculated using the saddle point equations ( 24) and ( 25) which lead to a parametric representation of ϕ(r) with the parameter α which varies in the range (0, ∞).These equations hold for r > r c .The critical value is for (A) and (B) and r c = 0 for (C) and (D).For 0 < r ≤ r c the free energy grows linearly with r The behaviour of ϕ(r) is illustrated in Fig. 2 for β = 7/2, 5/2, 3/2, 1/2 which are representative for the four cases.The curves are obtained by the parametric equations (51).The linear part of the solution (53), which corresponds to the condensed phase, is shown in dashed line.The derivative ϕ ′ (r) is shown in the right chart in Fig. 2. It is calculated from Eq. ( 24) which gives ϕ ′ (r) = K(α(r)).In combination with (25) this leads to the following parametric equations for the derivative for r > r c , with the parameter α which varies in the range (0, ∞).For 0 < r ≤ r c the derivative is constant as follows from (51) and (53).For (A) and (B) the derivative ϕ ′ (r) is constant for 0 ≤ r ≤ r c .For r → 1 the derivative ϕ ′ (r) tends to −∞ .For r → 0 it tends to K(0) which is finite for (A-C) and infinite for (D).The main difference between (A) and (B) is that the second derivative ϕ ′′ (r) is continuous at r c for (B) and discontinuous for (A).

VI. SINGULARITIES OF FREE ENERGY
In this section we will discuss singularities of the free energy density ϕ(r) at the critical point r = r c , This is the point where the solid line changes to the dashed line for curves (A) and (B) in Fig. 2. The analysis will be performed on the basis of parametric equations (51), which allow us to extract the singularity of the free energy density.The main components of these equations are the function K(α) and its derivative K ′ (α), so to prepare the ground let's analyze the behavior of these functions for α close to α c = 0, which correspond to r close to r c .
For the cases (D) and (C), that is for β ∈ (−∞, 2] there is no phases transition as r c = 0.So we will concentrate on the range β ∈ (2, ∞).For non-integer we can write and where the symbol ⌊x⌋ is the largest integer less or equal x.The constants µ c and r c on the right-hand side of these equations stand for µ c = K(0) and r c = −1/K ′ (0), respectively.The coefficients a k 's and b k 's as well as A and B can be directly deduced from the expansion of K(α) (47).They depend on β, but the form of this dependence is inessential for the further analysis.The signs of the coefficients are chosen for convenience.Only terms up to α β−1 in the first equation, and up to α β−2 in the second are displayed.All others are included in the little-o symbols for α → 0. We will separately consider the ranges β ∈ (2, 3) and β ∈ (3, +∞), and the borderline case β = 3. Let's begin with β ∈ (2, 3).The parametric equations (54) for the derivative ϕ ′ (r) become with for α close to α c = 0, as follows from (56) and (57).From the second equation we find that α = ((r − r c )/B) 1/(β−2) + . . .when r tends to r c from above.Plugging this into the first equation, we get where C = a 1 /B 1/(β−2) .The second derivative behaves as 3), so this means that ϕ ′′ (r + c ) = 0. Also for r → r − c the second derivative is equal zero ϕ ′′ (r − c ) = 0 so the second derivative is continuous for β ∈ (2, 3).In the same way it can be checked that the third derivative is continuous for β ∈ (2, 5/2) and it is discontinuous for β ∈ [5/2, 3).So the transition is the third order for β ∈ [5/2, 3).More generally, the n-th derivative ϕ (n) (r) is discontinuous at r c for β ∈ [2 + 1/(n − 1), 3).Therefore the transition is of n-th order for β ∈ [2 + 1/(n − 1), 2 + 1/(n − 2)) for n = 3, 4, . ... The order of the transition increases to infinity when β approaches two.Eventually at β = 2, the transition disappears completely.Now consider the range β ∈ (3, +∞).We first focus on β ∈ (3,4).In this case Eq. ( 59) takes the form For r close to r c this gives Substituting this into (58) we get We see that the second derivative is discontinuous at r c , because lim Therefore the transition is of the second order.We note in passing that the third derivative has an infinite discontinuity for β ∈ (3, 4), coming from the term (r − r c ) β−2 in (64).The calculations can be repeated for β ∈ (n, n + 1], for n = 3, 4, . .., to see that in all these intervals the second derivative has a finite discontinuity and that the derivative ϕ (n) (r) has an infinite discontinuity at r c .Now consider the borderline case β = 3.The parametric equations (54) can be written as and as follows from (50).Inverting the second equation for α we find where Substituting this into the first equation we get Taking the derivative of both sides we find that the second derivative has a logarithmic singularity for r → r To summarize, the phase transition is of the second order for β ∈ (3, +∞), with a finite discontinuity of ϕ ′′ (r) at the critical value r = r c .For β = 3, the phase transition is also of the second order but ϕ ′′ (r) has a logarithmic discontinuity log(r − r c ), as r approaches r c from above.For β ∈ (2, 3), the phase transition is of the third or higher order.The order of the transition increases when β approaches 2 and the transition eventually disappears for β = 2.There is no phase transition for β ≤ 2.

VII. GRAND-CANONICAL ENSEMBLE
So far we have considered a system of S particles in N boxes and analysed its behaviour in the thermodynamic limit (12) as a function of the limiting particle density ρ = 1/r = S/N .Now we will consider a system with a variable number of boxes N , which is controlled by "chemical potential" µ, which is equal to the energy that can be absorbed or released in the system due to a change of the number of boxes by one.The corresponding partition function is [3] (72) The system described by the partition function (72) will be called grand canonical ensemble.The grand-canonical averages are defined as A general comment on notation: we will distinguish between canonical averages and grand-canonical ones by the subscripts of the brackets which will be S, N in the first case, and S, µ in the second.As before, let us discuss the particle distribution.The grand-canonical average is as directly follows from the definition (73).The contribution to the sum from N = 1 has to be treated carefully, because if there is only one box it must contain all particles.Therefore for N = 1 the expression Z S−q,N −1 on the right hand side (74) should be interpreted as Z S−q,0 = δ S−q (see Eq. ( 3)).Taking this into account we can rewrite the last equation as for q = 1, 2, . . ., S. We are now interested in the behaviour of grand-canonical averages in the thermodynamic limit, S → ∞.To describe properties of the system in this limit we define the following thermodynamic potential Again, with a slight misuse of terminology, ψ(µ) can be called grand potential or the Landau free energy density, in analogy to the free energy density ϕ(r) (14) which was defined for the canonical ensemble.The two thermodynamic potentials are related to each other by the Legendre-Fenchel transform Indeed, approximating the grand-canonical partition function (72) by an integral and extracting its leading exponential behaviour ( 23) leads to (77).The inverse transform is In most cases it reduces to with µ(r) being a solution of Comparing ( 24) and ( 80) we find the following consistency equations from which we can deduce that the grand-canonical thermodynamic potential ψ(µ) (76) is the inverse function of the cumulant generating function ( 22) ψ(µ) = K −1 (µ), or equivalently that It follows that K ′ (ψ(µ)) = 1/ψ ′ (µ), which shows that also Eqs. ( 25) and (81) are consistent.Using this exact relation we can calculate the average r = r(µ) in the thermodynamic limit S → ∞ in the grand-canonical ensemble This result holds for µ ≤ µ c , where µ c is the critical value of the chemical potential given by where α c is given by (26).For power-law weights (46) α c = 0.When the chemical potential exceeds the critical value, µ > µ c , the saddle point equation breaks down and the average r = r(µ) (85) drops to zero.In this case N grows sub-linearly with S for S → ∞.
The box-occupation probability in the fluid phase in the grand-canonical ensemble can be calculated in the thermodynamic limit as follows.For large S (S ≫ 1) the partition function Z S,µ (78) can be approximated by e Sψ(µ) (78).Substituting this into (75) we get ⟨π(q)⟩ S,µ = w(q)e −qψ(µ)−µ . ( Using the consistency relations ( 82) and ( 83) we see that the canonical and grand-canonical averages of the box-occupation probability (33) and (87) are identical ⟨π(q)⟩ S,µ = ⟨π(q)⟩ S,N for S → ∞ when the particle density in the canonical ensemble is related to the chemical potential in the grand-canonical ensemble as follows N/S = −K ′ (K −1 (µ)) (85).The equivalence holds in the saddle point regime (fluid phase), that is for µ < µ c (86).It breaks down for µ ≥ µ c .
To conclude this section let us draw curves representing ψ(µ) and ψ ′ (µ) for power-law weights (46) for different β.As for the free energy density (see Fig. 2) we will show in each plot four curves for (A) β = 7/2, (B) β = 5/2, (C) β = 3/2, and (D) β = 1/2.The curves will be generated as parametric plots.To plot ψ(µ) we use the fact that ψ is the inverse function of K (84) which is equivalent to the following parametric equation with the parameter α > α c , which for power-law weights (46) α c = 0.For ψ ′ (µ) we have the following parametric equations: which are a direct consequence of ψ being the inverse function of K.The results are shown in Fig. 3.For power-law weights (46) the critical value of the chemical potential is Note that the parametric equations for ψ ′ (µ) (90) are almost identical to those for ϕ ′ (r) (54).There are two differences.The first is that there is a minus sign in front of 1/K ′ (α) in the equations for ψ ′ (r) which is absent in the equations for ψ ′ (µ).The second is that the right-hand sides of these equations are swapped which means that the ordinate and of abscissa switch roles.This is not surprising because ψ(r) and ϕ(µ) are related to each other by the Legendre-Fenchel transformation.In other words the drawings of ψ ′ (µ) and ϕ ′ (r) (compare Fig. 3 and Fig. 2) can be obtained from each another by swapping the vertical and horizontal axis, and changing the direction of the vertical axis.So far we have discussed the case of power-law weights (46).The question is what the corresponding plots look like for the family of weights (45).Again this can answered using the transformation w(s) → e λ+σs w(s) (5).Under this transformation the free energy density ( 14) transforms as and grand potential transforms as as follows from (77).The derivatives transform as and respectively.This means that the curve ϕ ′ (r) (see Fig. 2) moves up by λ and the curve ψ ′ (µ) (see Fig.
The parameter σ has no effect on the critical values r c and µ c .

VIII. SINGULARITIES OF GRAND POTENTIAL
In this section will determine singularities of grand potential ψ(µ) for the power-law weights (46) at the critical chemical potential µ c = log ζ(β) for β > 1.There is no phase transition for β ≤ 1.
For β ∈ (2, +∞) we can use Eqs.( 56) and (57) to cast the parametric equations (90) into the form and From the second equation we can calculate α = (µ c − µ)/a 1 + o(µ c − µ) for µ → µ − c .Substituting this into the first equation we see that which is a negative number.On the other hand lim so the first derivative has a discontinuity at µ = µ c , so the phase transition is of the first order.In addition, the first derivative has a power-law singularity (µ c − µ) β−2 , (or (µ c − µ) β−2 log(µ c − µ) for integer β), for µ → µ − cr .This singularity makes the n-th and higher derivatives infinite for µ → µ − c , where n = ⌊β⌋.For β = 2 there is no discontinuity of the first derivative.In this case the first derivative behaves as as µ approaches µ c , so it vanishes in the limit µ → µ − c .Taking the derivative of both sides we see that the second derivative diverges as ψ ′′ (µ) ∼ (µ c − µ) −1 , as µ tends to µ c , so the transition is of the second order for β = 2.

IX. OTHER ENSEMBLES
One can consider a system with a variable number of particles.In this case the partition function is It is easy to see that the partition functions factorises so it describes N independent boxes, each contributing an energy −K(α) to the total energy of the system.The contributions from different boxes are entirely independent.The average particle density is ρ(α) = −K ′ (α).It becomes infinite when α is equal to or greater than α c .This means that each box absorbs infinitely many particles from the reservoir.
In principle, also an ensemble with varying numbers of particles and boxes can be defined by the following partition function which is convergent for µ > K(α) and α > α c , yielding This expression becomes infinite for µ → K(α), the sum (112) diverges.From this partition function one can also calculate the grand-canonical partition function Z S,µ (72): as an inverse Laplace transform.

X. PROBABILISTIC AND QUASI-PROBABILISTIC NORMALISATION OF WEIGHTS
In some problems it is convenient to normalise the weights w(s) in the partition function (1) so that the model after the normalisation can be interpreted in a probabilistic manner.This can be done only if the sum in the denominator on the right hand side in (115) is finite.The normalisation can be obtained by the transformation (5) with σ = 0 and λ = − log ∞ q=1 w(q) and therefore it has no effect on the statistical averages (6).There are however families of weights which cannot be normalised because the normalisation sum in (115) is infinite.For example, the constant weights w(s) = 1 for s = 1, 2, . . .cannot be normalised.As a way round this problem, some authors have proposed a quasi-probabilistic normalisation [20] by putting an upper limit on the weights (and hence sum) for s = 1, 2, . . ., S, and w(s, S) = 0 for s > S. The normalisation constant Ω(S) depends on S. It is finite for any finite S. For example, Ω(S) = S for the trivial weights w(s) = 1.
The model with the quasi-probabilistic normalisation can be obtained from the original weights by the transformation (5) with λ = − log Ω(S).This leads to the following relationship between the partition functions and where the partition functions on the right hand side of the equations above are defined by weights independent of S, see Eqs. ( 1) and (72).
We note that the model with the quasi-probabilistic normalisation (116) is no longer invariant with respect to the transformation (5), because Ω(S) changes under this transformation.It therefore makes a difference whether one considers purely power-law weights q −β /Ω(S) or exponentially damped ones q −β e −σq /Ω(S).We stress that the weights w(s) in the partition functions (1) and (72) depend only on the occupation of the box while in the quasi-probabilistic model the weights w(s, S) depend also on S. The quasi-probabilistic model therefore belongs to a different class.
In the rest of this section we will consider the quasi-probabilistic model with the power-law weights (46) in the grandcanonical ensemble for µ = 0, which corresponds to the model discussed in the main part of this paper (72) with a running chemical potential µ(S) = log Ω(S) = log for S → ∞, with some positive coefficients c 1 , c 2 and c 3 depending on β.The detailed results, also for the borderline cases of β = 1 and β = 2, can be found in [20].
Let us stress that the model with weights normalised by the pseudo-probabilistic condition (116) corresponds to the model discussed in the main part of this paper (72) with a running chemical potential µ = log Ω(S) (119).For β > 1, the effective chemical potential µ = log Ω(S) approaches the critical value µ c = log ζ(β) from below as S increases.For β ≤ 1, the effective chemical potential tends logarithmically to infinity, for S → ∞.In both cases the system stays in the fluid phase, as long as S is finite.The effective particle probability distribution (87) for the system with the effective chemical potential µ = log Ω(S) is ⟨π(q)⟩ S,log Ω(S) = w(q) Ω(S) e −qψ(log Ω(S)) . (124) We thus see that the quasi-probabilistic normalisation introduces weak exponential damping for large q into the effective particle distribution.The damping factor decays with S but disappears completely only in the limit S → ∞.

XI. SUMMARY
After some general discussion of the urn model we have focused in this paper on the zeta-urn model which is analytically solvable.We have used this to determine the thermodynamic potentials that control the exponential growth of the canonical and grand-canonical partition functions in the thermodynamic limit and elucidate the critical behaviour.
The second derivative of the free energy density with respect to the particle density describes density fluctuations in the canonical ensemble.For β > 3 the second derivative is discontinuous at the critical point, so the transition is of the second order.For β = 3 the second derivative has a logarithmic divergence.For β ∈ (2, 3] the order of the transition changes at the discrete values β n = 2 + 1/(n − 1), n = 2, 3, . .., where n is the order of the transition for β ∈ (β n+1 , β n ].There is no phase transition for β ≤ 2. For the grand-canonical ensemble the situation is different.The first derivative of the corresponding thermodynamic potential is discontinuous for β > 2, so the phase transition is of the first order in this case.For β ∈ (1, 2] the order of the transition changes at the discrete values β ′ n = 1 + 1/(n − 1), where n = 2, 3, . . . is the order of the transition for β ∈ (β ′ n+1 , β ′ n ].There is no phase transition for β ≤ 1.The thermodynamic potentials for the canonical and grand-canonical ensembles can be derived from each other.More specifically, the function: ϕ ′ (r) on the support r ∈ (r c , 1) is the inverse function of µ → −ψ ′ (µ) on the support µ ∈ (−∞, µ c ) which is a consequence of the Legendre-Fenchel transform which relates the two functions.Additionally, we have shown the grand-canonical potential ψ(µ) on the support µ ∈ (−∞, µ c ) is the inverse function of the cumulant generating function K(α) on the support α ∈ (α c , ∞).
In [32] we use some of the results presented in this paper to study the Rényi entropy and diversity measures [20,31] for the zeta-urn model and discuss their singular behaviour and other properties in the thermodynamic limit (12) in the canonical ensemble.