Structural properties of generalised Planck distributions

A family of generalised Planck (GP) laws is defined and its structural properties explored. Sometimes subject to parameter restrictions, a GP law is a randomly scaled gamma law; it arises as the equilibrium law of a perturbed version of the Feller mean reverting diffusion; the density functions can be decreasing, unimodal or bimodal; it is infinitely divisible. It is argued that the GP law is not a generalised gamma convolution. Characterisations are obtained in terms of invariance under random contraction of a weighted version of a related law. The GP law is a particular instance of equilibrium laws obtained from a recursion suggested by a genetic mutation-selection balance model. Some related infinitely divisible laws are exhibited.


Introduction
In 1900 Planck derived a formula describing the energy spectrum of black body radiation. The energy density at frequency ν of a body at temperature T (in degrees Kelvin) is where h is Planck's constant, k is Boltzmann's constant, and c is the velocity of light. Good accounts of Planck's reasoning are given by Cropper (1970, Chap. 1), Longair (1984, Chap. 10) and, highly recommended, Weinberg (2013). Planck's first derivation combined reasoning based on electromagnetic theory and statistical mechanics, as well as empirical curve fitting to then new experimental results. In the last connection see Pais (1982, figure on p. 367). The form (1) suggests defining a Planck family of distributions whose members have a probability density function (pdf ) of the form where λ > 1, a > 0 and K is the normalisation constant. This family was introduced to the world in Johnson and Kotz (1970, p. 273), but subsequently dismissed from service by omission from the second edition of this reference work. The family (2) appears indirectly in Kleiber and Kotz (2003) in that its reciprocal form was used by Davis (1941) to fit personal income data for the U.S.A. in 1918. Subsequent analysis by Champernowne has shown that his own three-parameter distribution yields a better fit. See Kleiber and Kotz (2003, p. 240). Some elementary computations with (2) are reported by Tomovski et al. (2012). Finally, cases of (2) where λ = d + 1, where d = 1, 2, . . . , arise in physics literature as describing black-body radiation in a d-dimensional universe. See Puertas-Centeno et al. (2017) and its references. Nadarajah and Kotz (2006) have proposed the generalized Planck family GP(λ; a, b, c, ) whose pdf is where a, b > 0 and −1 ≤ c ≤ 1 with λ > 0 if c < 1, and λ > 1 if c = 1. The constant K is a normalisation factor which is evaluated at (6) below. Their motivation seems to be that the case c = 0 yields the gamma family, and hence that the generalised family gives a more malleable family of pdf 's. We will extend the range of the parameter c to allow −∞ < c ≤ 1. If λ = 1 then choosing 0 < c < 1 (resp. −1 < c < 0) yields the Bose-Einstein distribution (resp. the Fermi-Dirac distribution); see Kittel (1969, Chapter 9) or Olver et al. (2010, #25.12(iii)). In addition, the GP( 1 2 ; a, b, c) family is implicitly present in a specific 'house of cards' model for mutation-selection balance used by Turelli (1984, §3) to study heritable variation. This term is used by Kingman (1978) to describe the situation where the rate of mutation at a locus is independent of the parental allele, thus collapsing the ancestral history of selection at that locus. Nadarajah and Kotz (2006) derive formulae for moments, the characteristic function and entropies, they mention shapes of the pdf, and they look briefly at maximum likelihood estimation.
Our objective in this paper is to investigate some structural properties of the GP laws and point out limits of this investigation. In §2 we present some fundamental identities, a representation as a scale mixture of gamma laws (Theorem 1), and limiting behaviour as λ → ∞ (Theorem 2) or b → ∞ (Theorem 3). In §3 we summarise for the readers convenience properties of the gamma law we aim to extend to the GP laws.
The gamma family occurs as the limiting/stationary law of the Feller diffusion process, i.e., the continuous-state branching process with immigration. In §3 we show that the GP laws arise in the same way through an additive non-linear adjustment of the drift term. Nadarajah and Kotz (2006) list with no proofs how the shape of g(x) := g(x; λ, a, b, c) depends on values of the parameters. Their assertions about multimodality are vague. Modal properties of g are described in detail in §5. Theorem 4 is an exact statement about when g is, or is not, strongly unimodal. Theorems 5-8 examine its modal properties for the respective cases a = b and c = 1, c = 1, 0 < c < 1, and c < 0. In summary, depending in a fairly complicated way on parameter combinations, g either is decreasing, unimodal, or bimodal.
In recent years the infinite divisibility of many common distributions has been proved by showing that they are generalised gamma convolutions (Nadarajah and Kotz (2006) or Steutel and van Harn (2004)). However, as we argue at the end of §6, that if c = 0, then the GP laws are not in general a generalised gamma convolution.
Proposition 1 gives sufficient conditions for self-decomposability of GP laws in the restricted case where g is decreasing. The proofs verify Bondesson's (1987) sufficient conditions. Note that since self-decomposable laws are unimodal, the parameter combinations in §5 for bimodality preclude self-decomposability.
Many laws can be characterised as the fixed point of a transformation defined by weighting the law and then applying a random rescaling. Often the weighting operation is length, or size, biasing. More specifically for this case, let r be real. For any random variable Y ≥ 0 with distribution function G and finite moment function M G (r) = E(Y r ), define the order-r length-bias operator L r which maps G into the distribution function G r (x) = x 0 y r dG(y)/M G (r). We extend our notational convention so that L r Y denotes a random variable whose distribution function is G r . If r > 0 then L r Y ≥ st Y . This property opens the possibility of characterising L(Y ) through a relation of the form Y d = V L r Y where 0 ≤ V ≤ 1 is a random variable and the factors in the product are independent. We obtain such characterisations of GP laws in §7.
In §8 we study what amounts to a formal generalisation of Kingman's (1978) 'house of cards' model for genetic mutation-selection balance. Equilibrium laws of this model generalise (3) and a representation (Theorem 17) generalises Theorem 9. Finally, in §9 we exhibit some discrete and continuous infdiv laws which are related in certain ways to the generalised Planck laws. Longer proofs of our various results comprise §10.
If a random variable X has the pdf (3) then we write X ∼ GP(λ; a, b, c), and similarly for other laws. The standard gamma law with shape parameter λ (i.e., a = 1 and c = 0 in (3)) is denoted by gamma(λ), and we use the notation γ (λ) for a random variable having this law. Note the tail equivalence P(X > x) ∼ P(γ (λ)/a > x) as x → ∞. Similarly, the beta law with parameters α and β is denoted by beta(α, β). Independence of random variables X and Y is denoted by X ⊥ Y . Finally, L(Y ) denotes the law of a random variable Y, and There is some duplication of notation between sections, but no confusion will result from this.

Fundamental facts
The Lerch zeta function is defined by the series where θ > 0, and the series converges for all positive λ if |c| < 1, for λ > 1 if |c| = 1 and it is conditionally convergent if c = −1 and λ = 1. The function (1, λ, θ) is usually called the Hurwitz zeta function, and the Riemann zeta function ζ(λ) is the still more special case c = θ = 1. The function (θ + n) −λ is a completely monotone function of θ and in fact it is the Laplace transform of v λ−1 e −nv / (λ). Substituting the resulting integral representation of (θ + n) −λ into (4) and summing gives the known integral representation which clearly is well defined if −∞ < c < 1 and λ > 0, or if c = 1 and λ > 1, and θ > 0 in either case. We will take (5) as the definition of the Lerch function (or Lerch transcendent as it is called by some). See Gradshteyn and Ryzhik (1980, 9.550,6) for these formulae and Olver et al. (2010, Chapter 25) for graphs and references. There is a considerable literature devoted to representations and identities of the Lerch zeta as a function of λ, i.e., the integral (5) is regarded as a Mellin transform. More relevant for our purposes is the fact that (5) is a Laplace transform in the variable θ with c and λ understood as parameters.
Integrating the pdf (3) and setting α = a/b shows that the normalisation constant is It follows that the Laplace-Stieltjes transform of X ∼ GP(λ; a, b, c) is

Remark 1
The parameter b can be understood as a scaling parameter, and for most of our purposes it can be set to unity. In cases where this simplification is made, we denote the family by GP(λ; a, c) and the pdf by g(x; λ, a, c).
The quotient of gamma functions in (8) equals the moment function of γ (λ), and so one conjectures that the quotient of Lerch zetas also is a moment function. The series definition shows this is the case if c > 0, thus giving one proof of the following result which extends the mixture representation of the family (2) observed by Johnson and Kotz (1970). and Proof The pdf (3) can be written as a mixture of scaled gamma pdf 's, that is, X = γ (λ)/(a + n) with probability p(n). This is the essence of (9). The moment function (10) arises by multiplying (9) by (a + n) −t , summing, and appealing to (4).
Let N L denote a random variable having the Lerch law meaning that for n = 0, 1, . . . , P(N L = n) = p(n) as defined in (9). See Johnson et al. (2005, pp. 526-530) and the Pakes Journal of Statistical Distributions and Applications (2021) 8:12 Page 5 of 33 references there for information about the Lerch law and its special cases. Our designation 'Lerch law' for L(N L ), where N L = N − a, is adopted from Gupta et al. (2008). It is called a discrete Pareto law by others, e.g., Steutel and van Harn (2004, p. 68) for the case α = c = 1. It is shown in these references that L(N L ) is compound Poisson. The Lerch law has been shown to provide a better fit to certain biological count data than the generalised Poisson law and it is discussed as the limiting law of certain positive-recurrent birth and death processes by Klar et al. (2010). With this definition we see that the random variable N in Theorem 1 has a shifted Lerch law, N = N L + a.
Suppose that X i (i = 1, 2) has the GP(λ i ; a i , c i ) law with 0 ≤ c i ≤ 1, then it follows It follows that W ⊥ D, W has a type-2 beta law, and D has a discrete law which assigns positive mass to points comprising a countable dense subset of [ 0, ∞). This outcome was noted by Johnson and Kotz (1970) for the family (2). Theorem 1 has another interesting consequence. Let V have a standard normal law and V ⊥ γ (λ). Then V √ 2γ (λ) has the normal mixture law (Steutel and Van Harn (2004, p. 394)) whose characteristic function is (1 where γ 1 (λ) ⊥ γ 2 (λ). The law of L(λ) has been called a generalised Laplace law because setting λ = 1 yields the Laplace law; see Kotz et al. (2001, p. 180). Now let N have the shifted Lerch law (9) and N ⊥ L(λ). It follows that if X ∼ GP(λ; a, c) with the parameter restrictions in Theorem 1, then has the characteristic function (c, λ, a + t 2 )/ (c, λ, a).
The following result shows that if λ is large then L(X) is approximately normal. For this purpose we write X(λ) = X to emphasise the dependence of its law on λ. A simple proof follows from Theorem 1 in the case that c ≥ 0, but the proof we give in §10 for the general case follows from a probabilistic interpretation of the representation (5). Note that the outcome is independent of b and c.
the standard normal law.
The following result is an obvious consequence of (3).

Some properties of the gamma law
The pdf of the general gamma random variable γ (λ)/a is

Unimodality
The gamma law is unimodal with mode In addition, it is strongly unimodal if and only if λ ≥ 1; see Dharmadhikari and Joag-Dev (1998, p. 23).

Infinite divisibility
The gamma law is infinitely divisible, in fact self-decomposable for all values of the shape parameter, λ > 0. See Steutel and van Harn (2004, p. 225).

Diffusion limiting law
Let (Z t : t ≥ 0) be a diffusion process with Z 0 ≥ 0, drift function m(x), and diffusion function v(x). Assume that these functions are such that the weak limit Z ∞ of (Z t ) exists with pdf ζ(x). See Horsthemke and Lefever (1984, p. 110) for these matters. In this case the integrated stationary forward Kolmogorov equation is If ν > 0 and v(x) = 2x ν , then the drift function is determined by ζ(x) as In particular, if ζ(x) is the gamma pdf g(x; λ, a), then Feller (1951) solves the associated forward Kolmogorov equation in the case ν = 1, thus justifying calling (Z t ) the Feller diffusion. This process arises in a demographic context as a continuous-time and continuous sample path analogue of the Galton-Watson branching process with immigration. In fact, under suitable restrictions, it emerges as a diffusion limit of such processes as the initial population size grows unboundedly and the individual reproduction mean tends to unity. See Li (2011, Chapter 3).
Incorporating mean reversion is accepted as a fundamental requirement for models of the short/spot rate of interest rates (Brigo and Mercurio (2001), p. 46). The Feller diffusion arises in quantitative finance theory, where it is named the Cox-Ingersoll-Ross (CIR) model, as a popular way of incorporating linear mean reversion. See Cox et al. (1985), and Lamberton and Lapeyre (1996, p. 131) for a monograph treatment. Its time varying one-dimensional distributions are derived in this reference. The fact that the gamma law is the limiting/stationary distribution seems not to be clearly documented, but see Wong (1984), and it is an obvious consequence of formulae on page 131 in Lamberton and Lapeyre (1996). The case ν = 2 is a diffusion limit of the logistic model of population growth, for which see Horsthemke and Lefever (1984, §6.4) and references in Dennis and Patil (1984).

Characterisation
The fundamental characterisation property of the gamma law (Lukacs (1965)) asserts that if X and Y are independent positive random variables such that B = X/(X + Y ) and S = X + Y are independent, then X and Y have gamma laws with possibly different shape parameters, λ and μ, say. In this case B ∼ beta(λ, μ) and S ∼ gamma(λ + μ).
Many other characterisations are closely related to the Lukacs version. For example, we observe that length-biasing a gamma law yields another gamma law. More exactly (and recalling the length-bias operator L r ), if r > 0, then Hence one direction of the Lukacs characterisation can be expressed as a fixed-point . This is a characteristic property in the following sense (Pakes (1997, Theorem 4.1)).
Lemma 1 Let λ, r > 0. Any two of the following implies the third:

Planck laws as diffusion stationary laws
Let ζ(x) be the pdf of the gamma diffusion limit law in §3. Following Dennis and Patil (1984), let w(x) be a positive weight function such that is that of the limiting/stationary law for the diffusion with v(x) = 2x ν and infinitesimal mean function The In particular, if ν = 2 then the logistic per-capita growth rate r(x) is perturbed by a localised 'bump' which alters the limiting/stationary law in a quite significant manner. In general, the perturbation term lessens the equilibrium level, x e,w < x e (defined in §3.4), though much less so if x e is large. The rate of mean reversion is less for x < x e,w , but if c = 1 then the rate near the origin is close to λ + ν − 1 . The classification of {0} in Table 1 is unchanged if c < 1, but if c = 1 then the table entries are valid if λ is replaced by λ − 1 (which is positive in this case).
One aim of Dennis and Patil (1984) is exhibiting perturbations of the Feller/CIR diffusion whose limit laws are multi-modal, at least for certain parameter combinations. Their motivation is constructing stochastic versions of one-dimensional dynamical systems which have several equilibria. We examine modal behaviour of the GP laws in the next section. Nadarajah and Kotz (2006) list parameter conditions under which a Planck pdf is decreasing, unimodal (i.e., a unique positive mode), or 'maybe' multi-modal. The modal behaviour of GP laws also is of interest in statistical physics, e.g., Valluri et al. (2009) and references. We give here more precise statements, beginning with the more easily answered question of strong unimodality.

Strong unimodality
Recall that a pdf g(x) is strongly unimodal if and only if g is log-concave p. 17)). We work with the GP pdf g(x) = g(x; λ, a, c) and since this is twice continuously differentiable, then this condition is equivalent to where, for convenience, we write = λ − 1. Thus > 0 if c = 1 and > −1 if c < 1. The following result is specified in terms of the shape parameter λ. We write g ∈ SU to denote that g is strongly unimodal.

Unimodality
Strong unimodality implies unimodality, but the two notions differ and the algebraic work for establishing the latter is not as clean cut as for the former. In this subsection we work with the full family GP(λ; a, b, c), recalling that α = a/b. Denote the pdf by The following preliminary result is a simple deduction from (3). Proofs of the following results are in §10.

Lemma 2
Values of g at the origin are Setting y = bx, we see that the sign of g (x) is determined by the sign in (0, ∞) of We begin with the case of most physical significance, α = c = 1.
In this case the mode is (1)).
Remark 2 Spectroscopists define the normalised emissive power of a radiation distribution as the spectral density function divided by its modal value. Under the conditions of Theorem 5, this has the explicit form .

See Stewart (2012) (and its references) for the Planck distribution.
The next result relaxes the condition on α.

Example 1 Nadarajah and Kotz
is either unimodal or bi-modal.
Example 2 To see that bi-modality really is possible, we use parameter values based on those cited with no definite conclusion by Nadarajah and Kotz (2006): λ = 3/2, a = 1, b = 10 and c = 0.9.  (3) is decreasing for all α > 0 and the condition c ≥ −(1 − λ)e 2−λ is necessary for this outcome, and (iiib) If 0 < λ < λ c , then there is a critical value α c such that g(x) is decreasing if α > α c , and if 0 < α < α c , then there is an anti-mode M A and a mode M O such that 0 < M A < M O . (2006)

Infinite divisibility
We begin by proving a mixture representation for the case c > 0 which extends Theorem 1 and which can be regarded as an example of the general principle which can Table 2 Modes are M O1 = x 1 and M O2 = x 3 and the anti-mode M A = x 2 be inferred from Stuart (1962) and which is explicit in Steutel and van Harn (2004, Proposition 4.2, p. 344). This principle asserts that if 0 < λ < ν, then a scale mixture of the gamma(λ) law can be expressed as a scale mixture of the gamma(ν) law. We write the fundamental characterisation as γ (λ) The following mixture representation follows from Theorem 1.
and N has the shifted Lerch law as in Theorem 1. The pdf of V is The evaluation (19) results by differentiating the identity P substituting the algebraic form of the beta pdf and observing that The nature of f V (v) changes quite drastically as ν decreases to λ. If ν − λ > 1, this pdf is continuous; if ν − λ = 1 then it is decreasing with jump discontinuities everywhere on the support of N; if 0 < ν − λ < 1 then this pdf has singularities on the support of N.
If 0 < c < 1 and X ∼ GP(λ; a, c), then the representation assertion of Theorem 9 can be expressed as X d = BW where W = γ (ν)/N whose pdf is This defines a family of laws which includes the GP(λ; a, c) family as a special case; just observe that (c, 0, a) = 1/(1 − c).
In what follows we use the abbreviation infdiv to mean infinitely divisible. It follows from Theorem 9 that the GP(λ; a, b, c) law is a mixture of any gamma(ν) law with ν ≥ λ. Mixtures of gamma(2) laws are infdiv (Steutel and van Harn (2004, p. 346)), implying Assertion (i) in the following theorem. Similarly, any parameter configuration allowing the choice ν = 1 in Theorem 9 gives L(X) in the form of a mixture of exponential distributions and hence, as is obvious from (3), the pdf g(x) is completely monotone and hence convex. In fact, inspection of (14) shows that g(x) is log-convex if λ ≤ 1 and 0 < c ≤ 1. Log-convexity of a pdf is a sufficient condition that it be infdiv (e.g. Sato (2013, Theorem 51.4)) (and hence not strongly unimodal). Assertion (ii), stated for the case b = 1, extends a little the ambit of Assertion (i). Its proof is in §10.
Theorem 10 (i) If 0 ≤ c ≤ 1 and λ ≤ 2, then the GP(λ; a, b, c) law is infdiv.  (20) is equivalent to 0 < ν < 1 in which case L(X) is infdiv if It follows from Theorem 10 that the Bose-Einstein laws (λ = 1) are infdiv, but the status of the Planck law (λ = 4) is unknown. In passing we mention Varrò (2007) whose Section 5 is titled 'The infinite divisibility of the Planck variable' . His discussion distills to the assertion of the well-known fact that the geometric law is infdiv -he ascribes the title 'Planck-Bose' to this law.
Suppose that 0 < c ≤ 1, λ ≤ 2 and that b = 1, so X ∼ GP(λ; a, c) is infdiv. Note that the corresponding Lévy process is an infinite-activity (Type 1) subordinator, i.e., its sample paths are non-decreasing. It seems to be difficult to determine the nature of the Lévy measure of L(X). We can obtain some insight as follows. It is clear from its definition that the distribution family {GP(λ; a, c) : a > 0} is the natural exponential family generated by the measure whose density is where C(θ) = − log (c, λ, θ) is a Bernstein function, i.e., C (θ) is completely monotone (Schilling et al. (2012)). Thus there is a measure which, on the basis of the following discussion, we conjecture is absolutely continuous with density denoted by (x), such that It follows that the Laplace exponent of L(X) is C(a + θ) − C(a) and its Lévy measure is e −ax (x)dx. Next, we have Proceed formally by writing C (θ) = ∞ n=0 c n k n (θ) and equate the coefficients of c n in the identity obtained from (22) by multiplying through by (c, λ, θ) and using the representation (4). This yields n j=0 k j (θ) (n − j + θ) λ = λ (n + θ) λ+1 , allowing the step-by-step determination of the k n (θ). Thus k 0 (θ) = λ/θ, the Laplace transform of the measure λdx. This yields the Lévy measure (λ/x)e −ax dx for the case c = 0, i.e., for the gamma laws. Continuing, we find that Denoting the Kummer confluent hypergeometric function (Olver et al. (2010, p. 322)) by M (ψ, ρ, x), it follows easily from an integral representation that the Laplace transform of xM(1 + a, 2, −νx) is θ 1−a (ν + θ) −1−a . Thus k i (θ) (i = 1, 2) can formally be inverted to give the expansion However, the forms of k n (θ) for n ≥ 3 become more complicated and the corresponding inverse transforms involve convolutions of algebraic and Kummer functions with more and more factors.
We remind the reader that X has a self-decomposable (SD) law if for all ρ ∈ (0, 1] it has the auto-regression representation X d = ρX + Y ρ , where Y ρ ⊥ X. Self-decomposable laws are infdiv and they arise precisely as the possible weak limits of normed sums of independent random variables. We have seen that Planck laws can be bimodal under certain interacting parameter combinations, and since self-decomposable laws are unimodal, we may expect similar interacting parameter combinations to appear in the following partial result on self-decomposability of Planck laws. Its proof is in §10. Recall the constant c from Remark 4.

Remark 6
The bound for y c can be improved. If 0 < c < 1, it follows from (23) and y c < 2 that (y c − 2)e y c = −(2 + y c )c > −4c. Hence This lower bound subsists if c < 0 and λ c = 0 if c = c .   Table 3.

Example 4 Calculation yields the representative values of σ c and λ c displayed in
In the self-decomposable case of Proposition 1, the random variable X has a stochastic integral representation which we express as where (L t : t ≥ 0) is a Lévy process, which, in this case has non-decreasing sample paths (i.e., is a subordinator), and called the background driving Lévy process (conventionally abbreviated to BDLP); see Wolfe (1982). The Laplace exponent of the BDLP is C b (θ) = θC (θ). It follows from (22) that and hence from (43) below that In particular, C b (∞−) = λ, implying that the BDLP is a compound Poisson process. We infer too from Proposition 1 that the quotient (24) is a Bernstein function at least for the range of c delineated in the assertion. We have the following consequence, implying that X is the quotient of independent infdiv random variables.

The following result is a statement about what is sometimes called multiplicative infinite divisibility. Its proof is in §10.
Lemma 3 If c = 1, λ > 1 and b = a or b = 2a, then log X is infdiv.
What if c < 1? The answer is not known except when c = 0, the log-gamma law, := − log γ (λ). It follows from Euler's product representation of the gamma function that n j=1 γ j (1) This exhibits as the limit of a centred sum of independent gamma distributed random variables. The set of all such limit laws comprise the class of extended generalised gamma convolutions and they are SD laws. See Bondesson (1992, p. 112) and Steutel and van Harn (2004, p. 322). We remark in passing that a direct corollary of the derivation of this result is that has an extended generalised gamma convolution law with This shows in particular that the sequence {n! λ : n = 0, 1, . . . } is a Stieltjes moment sequence. Berg (2005, §2) has discussed this moment sequence in detail, showing in particular that it is Stieltjes-determinate if and only if λ ≤ 2. (Note that Berg denotes our parameter λ by c.) Observing that log(n + 1) − log n → 0 (n → ∞), we have the random series representation has the pdf e λ (x) identified in Lemma 2.1 of Berg (2005). Clearly B λ ×B ν d = B λ+ν , thus providing another proof that the family of pdf 's ( λ : λ > 0) comprises a product convolution semigroup of pdf 's. We note that Berg's apparently different analytical proofs of his results derive in part from a product representation of the gamma function.
We now recall that the weak limits of sums of independent gamma random variables (i.e., no centring) are called generalised gamma convolutions (GGC's) and they inherit the self-decomposability of gamma laws. Thus, the GGC class is closed under independent addition of random variables and, remarkably, it also is closed under independent multiplication (Bondesson (2015)).
The discrete-law version of GGC's are the generalised negative binomial convolution's (GNBC's), defined as the weak limits of sequences of sums of independent negative binomial random variables. The GNBC's coincide with the mixtures of Poisson laws where the mixing law is a GGC. Hence GNBC laws are infdiv. The next result records the known fact that a Lerch law is a GNBC; see Bondesson (1992, p. 135) or Steutel and van Harn (2004, p. 420). The proofs in these references are distributed over widely separate portions of their expositions, so it is worthwhile to record a consolidated proof. This proof is in §10 together with necessary concepts and facts.
Remark 7 A referee of an earlier version of this paper observed that the GNBC property may not hold if 0 < α, λ < 1, and indeed it does not if α = 0.5 and λ = 0.25.
It follows from Theorems 1 and 12 that X is equal in law to γ (λ) divided by an independent GNBC random variable. However this is not sufficient to conclude that L(X) is Pakes Journal of Statistical Distributions and Applications (2021) 8:12 Page 16 of 33 infdiv. Recent success in proving that many common laws are infdiv follows from showing that they are GGC's. An important proper subset of GGC's comprise absolutely continuous laws whose pdf 's are hyperbolically completely monotone (HCM). A function H(x) is called HCM if, for each u > 0, the product H(uv)H(u/v) is a completely monotone func- Bondesson (1992) or Steutel and van Harn (2004, §5 in Chap. VI). Any gamma pdf is HCM and many other pdf 's have this property. However, we have the following counter-result.

Theorem 13 The pdf g(x; λ, a, b, c) is HCM if and only if c = 0
Proof An HCM function extends as a holomorphic function to the complex plane cut along the negative reals (Bondesson (1992, p. 68)). But clearly, if c = 0 and ρ := − log |c|, the denominator of the pdf (3) has simple poles at

The assertion follows.
This outcome leaves open the possibility that a GP law is a GGC, but the following argument suggests not. Suppose that 0 < |c| ≤ 1. It follows from (7) that the moment generating function of the GP law is proportional to m(s) = (c, λ, α − s) and which, by virtue of (4), is analytic in the complex plane except for singularities at α, α + 1, . . . . If the GP law is a GGC, then for s > 0 its Thorin measure U(ds) is computed as U((0, s)) = π −1 A(s) where A(s) := arg(m(s)); see Bondesson (1992, p. 31). In particular A(s) is nondecreasing in (0, ∞). Clearly A(s) ≡ 0 in (0, α). Using the evaluation −1 = e −iπ , where i = √ −1, we see for α < s < α + 1 that where σ (s) = n≥1 c n (α +n−s) −λ is real and positive-valued in this interval if 0 < c ≤ 1, and negative-valued if −1 ≤ c < 0. Hence This is identically zero if λ is a natural number, so assume this is not the case. Letting s ↓ α shows that U(ds) assigns mass λ to s = α, exactly as for a gamma law, the case c = 0. However |σ (s)| ↑ ∞ as s ↑ α + 1 and hence tan A(s) → 0. Consequently A(s) cannot be non-decreasing in (α, α + 1), and this contradiction implies that the GP(λ; a, b, c) law is a GGC if and only if c = 0.

Characterisations by weighting and random scaling
In this section, factors in products of random variables are assumed to be independent. As mentioned in §1, many laws L(Y ) have been characterised in terms of a fixed point relation expressed as Y d = V L r Y , e.g., Pakes (1997). We begin by seeking to extend Lemma 1 to the Planck laws. We approach this by observing that if Z ≥ 0 is independent of Y and Hence characterising a product law such as GP(λ; a, b, c) can be achieved by attending separately to its factors. The gamma factor in Theorem 1 is addressed by Lemma 1, so we need only to consider the shifted Lerch zeta law divisor N. Pakes Journal of Statistical Distributions and Applications (2021) 8:12 Page 17 of 33 We will find it useful to highlight the dependence of L(N) on λ by writing N(λ) for N. It follows from (10) that the moment function of L −r N(λ) is According to (10) the right-hand side equals EN −t (λ + r), and hence a property which is reciprocal to (13). Let P denote the set of prime numbers and let {N p } be a sequence of independent random variables having the 'stretched' geometric laws Clearly Next, define independent random variables U p,r (p ∈ P) such that Finally, we state a useful general property of the length bias operator which extends Lemma 2.1 (c) in Pakes (1997).

Lemma 4 If Y is a positive random variable and r > 0, then L r Y s d = (L rs Y ) s for any real s such that E(Y rs ) < ∞.
Proof The definition of L r yields The following result is effectively a characterisation of L(N p ) in terms of an inverse power-law weighting.
Proof It follows from the evaluation (29) and using (30) to compute This identity implies that where the last step follows from Lemma 7.2 with s = −1.
Now assume that α = c = 1 in which case the moment formula (10) becomes ν(t) = ζ(λ + t)/ζ (λ). The Euler product formula and (28) yield the product representation Lemma 1, (31) and (32) together imply the following characterisation result for a subset of the Planck laws. The comments near the end of the Proof of Lemma 3 imply a similar result for the case α = 1 2 .
Theorem 14 does generalise Lemma 1, but at the expense of a complicated random scaling law and restrictions on the parameters. Another approach is to transfer complexity from the scaling to the weight function. Choose w(x) ≥ 0 for x ≥ 0 such that 0 < m w := E(w(X)) < ∞. Define the weighted distribution function L w G = x 0 w(y)dG(y)/m w and, if X has the distribution function G, then denote by L w X a random variable having the distribution function L w G. It is quite easy to show that if q ≥ λ > 0, V ∼ beta(q, 1) and Pakes and Navarro (2007) for this and similar results. The following result (proved in §10) extends this to a larger family of Planck laws. The restrictions on α and c ensure that the weight function is non-negative.
Theorems 1 and 9 can be interpreted in terms of weighting and scaling. The pdf of Y := γ (λ)/a can be regarded as arising from weighting the GP(λ; a, b, c) pdf; Y d = L w X where w(x) = 1 − ce −bx . So if W = α/N then the representation in Theorem 1 can be expressed as Since P(W ≤ 1) = 1 it follows that X is stochastically smaller than Y, an observation which also is a consequence of the fact that the weight function is strictly increasing. Similarly, if ν > λ then the representation of L(X) in Theorem 9 can be expressed as (33) where w(x) = (1 − ce −bx )x ν−λ and W = a/V .

A fixed-point generalisation
We extend the weighting operation defined in the previous section by allowing X to have an arbitrary distribution function G(x) supported in the real line and w(x) ≥ 0 to be a weight function satisfying 0 < m w := E[ w(X)] < ∞. We write L w X w to denote a random variable which has the weighted distribution function L w G( Finally let H(x) be an arbitrary distribution function and let 0 < β < 1 be a constant.
Here we are interested in fixed points of the transformation where c = (1 − β)/m w , thus ensuring that G(∞) = 1. A formal solution of (34) is where X(0) has a specified distribution function, I(n) ∼ Bern(β), the distribution function of V (n) is H, and I(n), L w X w (n) and V (n) are independent. The distribution function version of this recursion was analysed by Kingman (1978) in the particular case where all distribution functions are supported in [ 0, 1] and w(x) = x. This case is his 'house of cards' model for genetic mutation-selection balance where X(n) represents the fitness of an individual chosen randomly from the nth generation of an evolving population, the weight function represents the skewing effect of greater reproductive success of fitter individuals, β is a probability of mutation, and H expresses the fitness distribution of mutants.
The following result is a uniqueness assertion implying that G is determined by (34) for a given H. Proof Expand the denominator in (35) to obtain and interpret the summands in terms of quantities defined in the assertion.
The transformation Y = w(X) reduces the distribution function version of (36) to the recursion studied by Kingman (1978). The case c < 1 corresponds to his concept of 'democracy' , and then the distribution functions G(x; n) := P(X(n) ≤ x) converge in total variation norm to G(x) as given by (35). Moreover this convergence occurs geometrically fast. The case c = 1 is more problematic. The fact that G is a distribution function translates into the equality case of Kingman's inequality (4.1). Next, the discrete probability distribution (37) is non-increasing and hence lim sup n→∞ p n+1 /p n ≤ 1. This condition is sufficient to ensure that the renewal sequence {u n } generated by the discrete law {p n } satisfies the strong ratio limit property lim n→∞ u n+1 /u n = 1; see Garsia et al. (1962). Thus Kingman's meritocracy condition is fulfilled and hence we again have convergence to G in total variation norm. As a final remark, observe that the sequence {p n } is a Hausdorff moment sequence, and hence so is the renewal sequence {u n }. Kingman (1978) proves this, and so too does Horn (1970). It follows that {u n } is a Kaluza sequence (Kingman (1972, §1.5)), that is, {u n+1 /u n } is non-decreasing.
It is evident now that the GP(λ; a, b, c) law arises as the unique solution of (34) with the exponential weight in Theorem 17(ii), V (1) = γ (λ)/a and m w = βα λ (c, λ, α). We see too that the random scaling in Theorem 1 arises from the fact that members of the natural exponential family generated by a gamma law are simply scale changes of the generator. Finally, we observe that the case λ = 1 2 is implicit in Turelli (,§3) where V (1) has a centred normal law and the weight function has a Gaussian profile. The square of Turelli's limiting fitness random variable has a generalized Planck law with λ = 1 2 ; see his equation (3.6).

Related infinitely divisible laws
In this section we present some discrete and continuous infdiv laws which are related to the GP law in various ways. 1. Our first relation can be regarded as an alternative approach to the discussion of discrete Pareto laws in Steutel and van Harn (2004, p. 420). Recall that N L has the Lerch law, P(N L = n) = p(n) where {p(n) : n = 0, 1, . . . } is specified by the right-hand side of (9). Also, let (N (·)) denote a unit rate Poisson process. The following result adds detail to Theorem 11. Its proof is in §10.

Theorem 18
The Lerch law L(N L ) is a Poisson mixture, where T is independent of N and its density function is the mixture of exponential densities , where ε has a standard exponential law and it is independent of the gamma random variable γ (λ). Hence L(T) is infdiv.

2.
Recall that N (·) denotes the unit-rate Poisson process and that a positive random variable X has an infdiv law if and only if := N (ξ X) is infdiv (i.e., compound Poisson) for all ξ > 0 (Steutel and van Harn (2004, p. 369)). In view of the incompletely resolved infdiv status of the GP laws, it is of interest to see what emerges if X ∼ GP(λ; α, c). In this case It follows from Theorem 9 above that if 0 < c ≤ 1, then for any ν ≥ λ this is a mixture of negative-binomial laws with common shape parameter ν. In addition, it follows from Theorem 10 above and Theorem 12 in Steutel and van Harn (2004, p. 369) that L( ) is compound Poisson if λ ≤ 2, and that {p(j; ξ)} is log-convex if and only if λ ≤ 1; see Theorem 6.13 in Steutel and van Harn (2004, p.374).
It is easy to see that the sequence {ξ j (λ + j)/j! } is log-convex if and only if λ ≤ 1, so the following result is a little unexpected.

Pakes Journal of Statistical Distributions and Applications
(2021) 8:12 Page 22 of 33 We thus obtain a compound Poisson law by normalising, i.e., define the discrete law We remark in passing that the PGF Q(s) of this law can be expressed as a kind of quasimixture with respect to X ∼ GP(λ; α, c): , is the Mittag-Leffler function, i.e., a sub-family of confluent hypergeometric functions which generalises the exponential function; See Erdelyi (1955, p. 210) or Olver et al. (2010, p. 261).

3.
Let Y be a non-negative random variable with finite first order moment m Y , and in this subsection only we denote the ordinary size-biased version of Y by Y := L 1 Y . Then the law of Y is that of the stationary total lifetime in the renewal process generated by L(Y ).
It is known that the total lifetime can be represented as the sum of Y and an independent increment if and only where Z ≥ 0 and Y ⊥ Z. This is shown with the renewal context by van Harn and Steutel (1995), and independently by Pakes et al. (1996) emphasising the length-bias connection. Suppose that (40) holds and let ψ(θ) denote the cumulant function of L(Y ), where l Y is the infimum of the support of L(Y ) and (dx) is its Lévy measure, that is, a measure on I R + satisfying The following two examples explore the cases where one of L(Y ) or L(Z) has a GP law.
0 ≤ c ≤ 1 and 0 < λ ≤ 2, then it follows from Theorem 10(i) and (40) that where, from (22), the Laplace-Stieltjes transform of L(Z(λ)) is This emphasises that the question of whether the GP(λ; a, b, c) law is infdiv is equivalent to deciding for which values of λ is ξ(θ; λ) a completely monotone function of θ.
In preparation for our second example we recall the following facts. Pakes et al. (1996) observe that if Z = γ (λ)/α in (40) then L(Y ) has the Hougaard law H(η, δ, α), an infdiv law whose cumulant function is where η = 1 − λ < 1 and m Y = δα −λ . If 0 < η < 1 then the Hougaard laws for α > 0 comprise the natural exponential family generated by the positive stable law with index η; see Seshadri (1993). The case η = 0 is just the gamma law Y = γ (δ)/α, and if η < 0 then L(Y ) is a compound Poisson law which has a gamma jump law. Suppose in (40) that the distribution function of Z has the mixture form where {p(n)} is an arbitrary discrete law and G n is a distribution function satisfying G n (0) = 0. Theorem 2.2 in Pakes et al. (1996) asserts that Y in (40) can be expressed as where the summands are independent and L(Y n ) is infdiv with cumulant function Example 6 Choosing Z ∼ GP(λ; α, c) in (40), then Theorem 1 implies the form (41) where the p(n) are given by (9) and G n (x) is the distribution function of γ (λ)/(α + n). It follows that (42) holds with Y n ∼ H(η, δ n , α + n) where

and the cumulant function of Y can be shown to be
It follows from the form of the Lévy density in this representation that L(Y ) is a compound Poisson law if 0 < c < 1 and λ > 1, or c = 1 and λ > 2, and in either case the jump law is the GP(λ − 1; α, c) law; if 0 < c ≤ 1 and 1 < λ ≤ 2, then L(Y ) belongs to the Bondesson class of laws (i.e., the smallest class containing exponential mixtures and closed under convolution and weak limits); and if 0 < λ < 1 and 0 ≤ c ≤ 1, or if λ = 1 and 0 ≤ c < 1, then L(Y ) is a GGC law.

Proofs
Proof of Theorem 1 The change of variable y = αv in (5) yields It follows that the Laplace-Stieltjes transform of (aX( The exponents in the quotient of expectations converge to infinity in probability and hence this quotient converges to unity. The first factor on the right-hand side converges to  (14) holds. If λ < 1, i.e., < 0, then C(x) > 0 and Assertion (i) follows.
Suppose that c < 0. Then C(x) < 0 if λ ≥ 1. If, instead, 0 < λ < 1, then < 0 and then (14) can be recast as − This clearly cannot hold because, as x ↑ ∞, the right-hand side decreases to zero much faster than the left-hand side. Assertion (ii) follows. If c > 0 and λ ≥ 1, i.e. ≥ 0, then C(x) can be factorised as The first factor is positive, so (14) holds if and only if the second factor is negative-valued, i.e., The right-hand side is convex increasing from 1 − c > 0 at y = 0. The parameter decreases through (0, ∞) as increases from zero. Thus (46) holds if is sufficiently large. On the other hand, equality in (46) holds for two distinct values of y if is small. It follows that there is a critical value ∈ (0, ∞) such that the equality has a unique positive solution ϒ for y. This configuration is equivalent to the conditions 2 y = e y − ce −y and 2 = e y + ce −y ; the second condition is a tangent condition. Eliminating yields y = r(y) := e 2y − c e 2y + c .
Observing that r(0) = ε ∈ (0, 1), that r(y) ↑ 1 as y ↑ ∞, and that r(y) is concave increasing, it follows that the unique solution ϒ ∈ (ε, 1). Hence the critical value of r is Observe that the right-hand side lies in (−e −1 , 0) and recall that the functional equation we w = −z, where −e −1 < z < 0, has two real-valued solutions, the principal and secondary Lambert functions of −z. The solution y relevant for the present context tends to zero as → 0, and hence −(y − ) = W (− e − ), where W(·) denotes the principal Lambert function. Assertion (ii) now follows because M O = ay. Assertion (iii) follows from the expansion W (x) = x + o(x) as x → 0.

Proof of Theorem 6
The relation (17) Clearly R(0) = 1 and R(y) decreases convexly to zero as y ↑ ∞. Assertion (i) follows from graphical considerations. Assertion (ii) follows similarly by considering intersections of the graphs of each side of (49). In particular, for a fixed the critical value α c corresponds to a solution of (49) subject to the tangent condition −α = R (y). Eliminating α yields the equation z −1 sinh z = − 1 2 , where z = 1 2 y. Denoting the solution by z, (49) yields the critical value as Observing that R (0) = − 1 2 , it follows that α c < 1 2 . The reasoning is similar if = 1 and the assertion follows because the two graphs intersect at y = 0.
If 2α < , then S(·) has a unique inflection point y i = ( /α) − 2 ∈ (0, ∞) and it is convex in [ 0, y i ] and concave in [ y i , ∞). Now So if > , then S(·) is convex increasing in [ 0, y i ]. It follows that S(·) has a unique global maximum, at y M say, and a unique zero y 0 with y i < y M < y 0 < / . Hence g(x) is unimodal. If < , then S(·) is convex decreasing near y = 0. The above assertions about y M and y 0 hold (though y 0 may not be the only zero), and now S(·) has a local minimum y m ∈ (0, y i ). If S(y m ) < 0, then there are a further two zeros y i0 such that 0 < y 10 < y m < y 20 < y M and g(x) is bimodal in such a case. c.f. (23). It follows, as above, that y c > 2 and that