The Role of the Central Limit Theorem in the Heterogeneous Ensemble of Brownian Particles Approach

Abstract

The central limit theorem (CLT) and its generalization to stable distributions have been widely described in literature. However, many variations of the theorem have been defined and often their applicability in practical situations is not straightforward. In particular, the applicability of the CLT is essential for a derivation of heterogeneous ensemble of Brownian particles (HEBP). Here, we analyze the role of the CLT within the HEBP approach in more detail and derive the conditions under which the existing theorems are valid.

1. Introduction

The heterogeneous ensemble of Brownian particle (HEBP) [1,2] describes a large class of anomalous diffusion phenomena, observed in many physical and biological systems [3,4,5,6]. The HEPB approach is based on the Langevin stochastic equation of diffusion of a free particle, i.e., the mesoscopic description of Brownian motion (Bm). The relaxation time ($\tau$) and the diffusivity ($\nu$) of the particle constitute two important scales of Bm process. In the classic Langevin approach, $\tau$ and $\nu$ are constant parameters. Instead, in the HEPB approach, it is assumed that $\tau$ and $\nu$ are time-independent random variables. In HEPB the single-particle trajectory (SPT) follows a classic Langevin dynamics and it is characterized by a stochastic realization of the parameters (${\tau }_i$ and ${\nu }_i$ for particle i). The random nature of the scales $\tau$ and $\nu$ of the SPTs mimics the heterogeneity of the environment and/or the heterogeneity of an ensemble of particles diffusing in the environment. In fact, the anomalous diffusion behavior described by HEPB is generated by the heterogeneity of $\tau$ and $\nu$ values in different SPT realizations.

Long time and space correlation, characteristics of many anomalous diffusion processes [7,8,9], are often introduced by modifying the laws of the dynamics, by including memory kernels and/or integral operators [10,11] in the equations, for example, the fractional derivatives [12]. These changes in the dynamics introduce often non-Markovianity and/or non-locality in the processes.

The HEPB approach maintains the Markovianity of the process because the fundamental process remains the classical Brownian motion (Bm), but the heterogeneity of the scales involved in the system permits to describe a process with stationary features different from the Bm [13]. We will refer to this heterogeneity as to a population of scales, because the values of these scales follow a given probability distribution. Furthermore, the model structure permits to keep the standard dynamical laws, with integer time derivative of physical variables like the velocity (V) and the spatial coordinate (X) of the particle, and to avoid the introduction of fractional time derivatives.

One of the random scales contributing to the anomalous behavior the Langevin description of HEPB [1] is the time scale $\tau$. When the distribution of $\tau$ is properly chosen and $\nu$ is kept constant, the HEPB describes the same one-time one-point probability density function (PDF) of the fractional Brownian motion (fBm), i.e., a normal distribution with variance (the mean squared displacement of the process, MSD) scaling as a power law of time in the long time limit:

$${\sigma }_x^2\left(t\right)=\langle {(x(t+t_0)-x\left(t_0\right))}^2\rangle =D_{\alpha }{t}^{\alpha },$$

(1)

where $0<\alpha \le 2$ and $D_{\alpha }$ is the constant playing the role of diffusion coefficient. Depending on the value of the exponent $\alpha$, it is possible to distinguish what is called super-diffusion and sub-diffusion, associated respectively to super-linear and sub-linear values of the parameter.

The convergence of the PDF to a normal distribution depends on the applicability of the classical central limit theorem (CLT). We will demonstrate later that by choosing properly the population of the time scales according to certain PDFs, both the Gaussian shape of the PDF and the anomalous scaling of the variance can be guaranteed.

The CLT represents a cornerstone in probability theory. It states that when a large amount of one -or multi-dimensional, real-valued and independent (or weakly dependent [14]) random variables are summed, the probability distribution of their sum will tend to the Gaussian distribution $\mathcal{G}$, defined by its characteristic function:

$${\widehat{g}}_{\mathcal{G}}\left(k\right)=exp(-i\mu k-\frac{k^2\sigma }{2}).$$

(2)

This result has been generalized by the generalized CLT to a larger class of stable distributions described by the following characteristic function [15]:

$${\widehat{g}}_{\alpha }\left(k\right)={exp(-i\mu k-C|k|}^{\alpha }[1+i\beta \left(sign\left(k\right)\right)\omega (k,\alpha )]),$$

(3)

where $\alpha ,\beta ,\mu ,C\in \mathbb{R}$, $\omega (k,\alpha )=tan(\alpha \pi /2)$ if $\alpha \ne 1$, else $\omega (k,\alpha )=2/\pi ln\left(\right|k\left|\right)$. The Gaussian distribution can be found to be a special yet fundamental case when $\alpha =2$.

The generalized CLT [15] describes the convergence of the sum of stable variables with also infinite variance, for example, the symmetric Levy stable distribution. The stability property of the symmetric Levy stable distribution is fundamental to obtain a random walk with infinite large displacements as the well known Lévy–Feller diffusion process [8,16,17]. The PDF of this process converges in fact to the non-Gaussian but symmetric Levy stable distributions.

In the following sections, we briefly review the CLT formulation, then we introduce the problem of the convergence in the distribution of a mixture of Gaussian random components with random variances when the variance distribution is particularly extreme. Thereafter we recall the HEBP model formulation and define the sufficient conditions over $\tau$ to obtain a fBm-like process.

2. The Classical CLT Formulation

For completeness, we summarize here the most famous versions of the CLT and introduce some useful notation and definitions.

For parameters $\mu \in \mathbb{R}$ and $\sigma \in {\mathbb{R}}_{+}$, a normal (or Gaussian) distribution $\mathcal{N}(\mu ,{\sigma }^2)$ is a continuous probability distribution defined by its density function

$$f(x\mid \mu ,{\sigma }^2)=\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{-\frac{{(x-\mu )}^2}{2{\sigma }^2}},$$

(4)

where $\mu$ and ${\sigma }^2$ are the expectation and variance of the distribution, respectively. For $\mu =0$ and $\sigma =1$, we obtain what is usually called the standard normal distribution $\mathcal{N}(0,1)$.

For the sequence of random variables ${\left(X_n\right)}_{n\ge 1}$, we define random variables ${\left(S_n\right)}_{n\ge 1}$ as partial sums $S_n=X_1+X_2+\cdots +X_n$. The theory of central limit theorem derives conditions for which there exist sequences of constants ${\left(a_n\right)}_{n\ge 1}$, $a_n>0$, and ${\left(b_n\right)}_{n\ge 1}$ such that the sequence ${\left(\frac{S_n-b_n}{a_n}\right)}_{n\ge 1}$ converges in distribution to a non-degenerate random variable. In particular, CLT describes the convergence to standard normal distribution with constants defined as $a_n^2={\sum }_{k=1}^{n}Var\left[X_k\right]$ and $b_n={\sum }_{k=1}^n\mathbb{E}\left[X_k\right]$.

Different constraints on variables $X_1,X_2,\dots$ lead to different versions of the CLT. We will briefly review the most prominent results of the theory of central limit theorems. For a more pedagogical and/or historical perspective, see [18,19,20,21,22,23,24].

We start with the case when variables $X_1,X_2,\dots$ are independent and identically distributed. With additional requirements of finite mean $\mu$ and positive and finite variance ${\sigma }^2$, we obtain:

$$\frac{S_n-n\mu }{\sigma \sqrt{n}}\stackrel{d}{\to }\mathcal{N}(0,1)\mathrm{as}n\to \infty .$$

(5)

Dealing with independent, but not necessary identically distributed, random variables $X_1,X_2,\dots$ with finite variance, we define ${\mu }_k=\mathbb{E}{X}_k$, ${\sigma }_k^2=Var{X}_k$ and $s_n^2={\sum }_{k=1}^n{\sigma }_k^2$ for every $k\ge 1$. To obtain the main result, we need two Lindeberg’s conditions:

$$L_1\left(n\right)=\underset{1\le k\le n}{max}\frac{{\sigma }_k^2}{s_n^2}\to 0\mathrm{as}n\to \infty ,$$

(6)

and

$$L_2\left(n\right)=\frac{1}{s_n^2}\sum _{k=1}^n\mathbb{E}{\left|X_k-m_k\right|}^{2}I\left\{\left|X_k-m_k\right|>ϵs_n\right\}\to 0\mathrm{as}n\to \infty \left(\mathrm{for}\mathrm{every}ϵ>0\right).$$

(7)

The Lindeberg–Lévy–Feller theorem provides sufficient and necessary conditions for the following result:

$$\frac{S_n-\mathbb{E}{S}_n}{s_n}\stackrel{d}{\to }\mathcal{N}(0,1)\mathrm{as}n\to \infty .$$

(8)

Lindeberg and Lévy proved (using different techniques) that if (7) holds, so do (6) and (8). Feller proved that if both (6) and (8) are satisfied, then so is (7).

Since Lindeberg’s condition (7) can be hard to verify, we can instead use the Lyapounov’s condition which assumes that for some $\delta >0$, $\mathbb{E}{\left|X_k\right|}^{2+\delta }<\infty$ (for all $k\ge 1$) and

$$\frac{1}{s_n^{2+\delta }}\sum _{k=1}^n\mathbb{E}{\left|X_k-{\mu }_k\right|}^{2+\delta }\to 0\mathrm{as}n\to \infty .$$

(9)

If for independent random variables $X_1,X_2,\dots$ the Lyapounov’s condition is satisfied, then the central limit theorem (8) holds. Since the Lyapounov’s condition implies the Linderberg’s second condition this result follows directly from the Lindeberg–Lévy theorem.

In all versions of the CLT mentioned so far, the assumption of finite variance was crucial. To extend our observations to the case when variance does not exist (or is infinite), we introduce the notion of domains of attraction. We are observing a sequence $X,X_1,X_2,\dots$ of independent, identically distributed random variables. We say that X, or, equivalently, its distribution function $F_X$, belongs to the domain of attraction of the (non-degenerate) distribution G if there exist normalizing sequences ${\left(a_n\right)}_{n\ge 1}$, $a_n>0$, and ${\left(b_n\right)}_{n\ge 1}$, such that

$$\frac{S_n-b_n}{a_n}\stackrel{d}{\to }G\mathrm{as}n\to \infty .$$

(10)

Another important concept is the one of stable distribution. Retaining the same notion, the distribution X is stable if there exist constants ${\left(c_n\right)}_{n\ge 1}$, $c_n>0$, and ${\left(d_n\right)}_{n\ge 1}$ such that $S_n\stackrel{d}{=}{c}_{n}X+d_n$ (for all $n\ge 1$).

It can be shown that only stable distributions possess a domain of attraction [18]. The most notable stable distribution is Gaussian and by the classical CLT we know that all distributions X with finite variance belong to the domain of attraction of the Gauss Law. However, there are also some distributions with infinite variance that belong to it. More precisely, it can be shown [25] that random variable X with the distribution function $F_X$ belongs to the domain of the attraction of the Gauss law if and only if

$$\underset{x\to +\infty }{lim}\frac{x^2\left[1-F_X\left(x\right)+F_X(-x)\right]}{{\int }_{-x}^x{t}^{2}d{F}_X\left(t\right)}=0.$$

(11)

3. CLT for a Population of Gaussian Random Variables

We reviewed the fundamental theorems related to the classical CLT, having the Gaussian distribution as limit distribution of the sum of random variables $S_n$. The recurrent and sufficient (but not necessary) condition leading to the classical CLT description is that the variance of the i.i.d. random variables that are summed should be finite. However, there exist distributions with infinite variance that fall in the Gaussian domain of attraction [15,25]. In this paragraph, we provide a preparatory example to introduce the role of the CLT in the HEBP. The sum of a population of Gaussian variables with random variances (which may tend to infinity), is here rewritten as the sum of i.i.d. random variables defined as the mixture of random Gaussian components with random variances. If such a mixture has finite variance, the standard CLT conditions are satisfied. In fact, as it will be explained in more details hereafter, the convergence in distribution of the sum to a Gaussian is not always guaranteed when some extreme distribution of the random variance is considered.

Let us consider partial sums of independent Gaussian random variables

$$S_n=\sum _{k=1}^n{X}_k,$$

(12)

where, denoting with $f_k\left(x_k\right)$ the PDF of $X_k$, we have:

$$f_k\left(x\right)\sim N(0,{\sigma }_k^2).$$

(13)

The distribution of the sum of n random variables can be exploited in term of a convolution integral. Thus, we can derive explicitly the limit distribution of Equation (12) by inverting the characteristic function $\varphi \left(\omega \right)$ of $S_n$, which corresponds to the product of the characteristics ${\varphi }_k\left(\omega \right)$ of $X_k$:

$$\varphi \left(\omega \right)={\mathsf{\Pi }}_{k=1}^n{\varphi }_k\left(\omega \right),$$

(14)

which gives

$$\begin{array}{ccc}\hfill \varphi \left(\omega \right)& =& {\mathsf{\Pi }}_{k=1}^n\left(e^{-\frac{{\omega }^2}{2}{\sigma }_k^2}\right)\hfill \end{array}$$

(15)

$$\begin{array}{ccc}\hfill & =& e^{-\frac{{\omega }^2}{2}{\sum }_{k=1}^n{\sigma }_k^2}.\hfill \end{array}$$

(16)

Let us assume ${\sigma }_k\sim \sqrt{\mathsf{\Lambda }}$, with $\mathsf{\Lambda }$ distributed according to a generic PDF $f\left(\lambda \right)$. If the first moment of $\mathsf{\Lambda }$ exists in the limit of large n, by applying the law of large numbers, we can well approximate the Equation (16) in terms of $\mathbb{E}\mathsf{\Lambda }$:

$$\varphi \left(\omega \right)=e^{-\frac{{\omega }^2}{2}·n·\mathbb{E}\mathsf{\Lambda }},$$

(17)

which is indeed the characteristic function of a Gaussian distribution with variance $n·\mathbb{E}\mathsf{\Lambda }$ for finite expectation of $f\left(\lambda \right)$ even if the supremum of $\mathsf{\Lambda }$ does not exists.

The convergence of $S_n$ can be proven using the CLT for the sequence of independent, identically distributed random variables $X,X_1,X_2,\dots$ with $X\sim \mathcal{N}(0,\mathsf{\Lambda })$. These variables, in general, will not have a Gaussian shape and can equivalently be defined as the product of the independent random variables:

$$X=\sqrt{\mathsf{\Lambda }}·Z,$$

(18)

where $Z\sim f_1\left(z\right)=N(0,1)$, $\mathsf{\Lambda }\sim f_2\left(\lambda \right)$, $\mathsf{\Lambda }\in {\mathbf{R}}_{+}$. The PDF $f\left(x\right)$ of X can be represented by the integral form [26]

$$f\left(x\right)={\int }_0^{\infty }{f}_1(x/\sqrt{\lambda })f_2\left(\lambda \right)\frac{d\lambda }{\sqrt{\lambda }}.$$

(19)

Since Z is a Gaussian distribution, it follows that $\frac{1}{\sqrt{\lambda }}{f}_1(x/\sqrt{\lambda })=N(0,\lambda )$. Using Fubini’s theorem, now it is easy to compute the second moment of X:

$$\begin{array}{ccc}\hfill VarX& =& {\int }_{-\infty }^{\infty }{x}^2{\int }_0^{\infty }{f}_1(x/\sqrt{\lambda })f_2\left(\lambda \right)\frac{d\lambda }{\sqrt{\lambda }}dx\hfill \end{array}$$

(20)

$$\begin{array}{cc}=& {\int }_0^{\infty }\lambda f_2\left(\lambda \right)d\lambda =\mathbb{E}\mathsf{\Lambda }.\hfill \end{array}$$

(21)

If $\mathbb{E}\mathsf{\Lambda }$ is finite the partial sums $S_n=X_1+\cdots +X_n$ of i.i.d. random variables $X_k$ converge in distribution to a Gaussian

$$\frac{S_n}{\sqrt{n}}\stackrel{d}{\to }\mathcal{N}(0,\mathbb{E}\mathsf{\Lambda }).$$

(22)

If $\mathbb{E}\mathsf{\Lambda }$ is not finite, the distribution $f\left(x\right)$ may fall out of the Gaussian domain of attraction. For example, by choosing $\mathsf{\Lambda }$ (the random variance) to be the extremal Lévy density distribution, it follows that $f\left(x\right)$ (the mixture defined by Equation (19)) corresponds to the symmetric Lévy stable distribution [27]. In fact in the case $f\left(x\right)$ is itself a stable distribution, like the Levy stable distribution is, its sum belongs to its own domain of attraction.

However, infinite variance is not a synonym of stability. In fact, despite the presence of infinite $\mathbb{E}\mathsf{\Lambda }$, under certain constraints on the tail of the distribution $f\left(x\right)$, $f\left(x\right)$ still satisfies (11) and falls in Gaussian domain of attraction, for example if its PDF for large x is proportional to $x^{-3},x^{-3}log\left(x\right),x^{-3}/log\left(x\right)$ [15].

4. Application of the CLT in the HEBP

In the HEBP Langevin model [1] the anomalous time scaling of the ensemble-averaged MSD is generated by the superposition of a population of Bm processes in a similar way to equation (12), where each single process is characterized by its own independent timescale, and with frequency of appearance of such timescale described by the same PDF.

CLT applicability guarantees that after averaging over a properly chosen timescale distribution the shape of the PDF will remain Gaussian despite the time scaling will change from being linear in time to be a power low of time in the long time limit, following Equation (1). In order to show this applicability let first introduce the HEPB construction.

Let us start with the classic Langevin equation describing the dynamics of a free particle moving in a viscous medium (or Bm):

$$dV=-\frac{1}{\tau }Vdt+\sqrt{2\nu }dW,$$

(23)

where V is the random process representing the particle velocity, $\tau$ in classical approach corresponds to the characteristic time scale of the process, i.e., the scale of decorrelation of the system. In the classic Langevin description the timescale is defined by the ratio $\frac{m}{\gamma }$, with m being the mass of the diffusing particle and $\gamma$ the Stoke’s drag force coefficient of the velocity. The multiplicative constant of the Wiener noise increment $dW$ in the square root, $\nu$, represents the velocity diffusivity and is related to the drag term by the fluctuation-dissipation theorem (FDT) [28]. This relation does not depend on the mass of the particle but on the average energy of the environment (the fluid) and the cross-sectional interaction between the medium and the particle moving. The Wiener increment $dW$ is the increment per infinitesimal time induced by the presence of a Gaussian white noise with unit variance and is hence fully characterized by its first two moments:

$$\langle dW\left(t\right)\rangle =0,\langle dW{\left(t\right)}^2\rangle =t.$$

(24)

The presence of Gaussian increments in the stochastic equation leads to the stationary state $V\sim N(0,kT/m)$ and, being $V=dX/dt$, to the stationary increments process $X\left(t\right)\sim N(x_0,{\sigma }_x^2\left(t\right))$, with ${\sigma }_x^2\left(t\right)=\nu {\tau }^{2}t$.

Let now the parameters $\nu$ and $\tau$ be time independent random variables: $\nu \sim p_{\nu }\left(\nu \right)$ and $\tau \sim p_{\tau }\left(\tau \right)$. The way it will affect $V\left(t\right)$ and $X\left(t\right)$ is clear in the case of $\nu$, but more complicated to specify in the case of $\tau$.

Let us consider the velocity defined as a product of random variables $V=\sqrt{\nu }{V}^{\prime }$. It is easy to show that $\sqrt{\nu }$ factorizes out from the stochastic differential equation, resulting in the following description of the evolution of $V^{\prime }$:

$$d{V}^{\prime }=-\frac{1}{\tau }{V}^{\prime }dt+\sqrt{2}dW.$$

(25)

Therefore, the PDF associated to the processes $V\left(t\right)$ and $X\left(t\right)$ can be derived by applying the same integral formula of Equation (19), eventually producing non-Gaussian PDF and weak ergodicity breaking stochastic processes as result [29,30,31].

Dealing with random timescales is much more tricky because the variable $\tau$ is embedded in the correlation functions and is not possible to factorize it out without simultaneously transforming the time variable. Furthermore, because of the time variable transformation different realizations of the process would not be comparable directly anymore without reverse transformation.

To avoid these complications, we define $V^{\prime }$ as the superposition of $N_{\tau }$ independent Bm processes each with its own timescale:

$$V^{\prime }\left(t\right)=\frac{1}{N_{\tau }}\sum _{\tau }{V}^{″}\left(t\right|\tau ),$$

(26)

where $V^{″}\left(t\right)$ can still be described by the Equation (25). If the resulting process $V^{\prime }\left(t\right)$ is still a Gaussian process, the previously described approach to derive $V=\sqrt{\nu }{V}^{\prime }$ can be applied without further changes. However, all the correlation functions of $V^{\prime }$ and moments will become the sum of the correlation functions of the single processes $V^{″}\left(t\right|\tau )$, which is equivalent to averaging with respect to $p_{\tau }\left(\tau \right)$. A careful choice of this distribution permits to obtain non-linear time scaling of the MSD of $V^{\prime }$.

Let us demonstrate the applicability of the CLT explicitly making use of the Equation (17). Assuming that a global stationary state (in the sense of stationary increments) has been reached, the relation between the MSD and the VACF determined by the free particle Langevin dynamics can be expressed by:

$${\sigma }_x^2(t,\tau )=2{\int }_0^t(t-s)R(s,\tau )ds,$$

(27)

where $R(t,\tau )=\nu \tau e^{-t/\tau }$, with $\nu$ being an arbitrary constant, is the stationary VACF of the process associated to the realization $\tau$ of the timescale, $V^{″}\left(t\right|\tau )$.

By omitting time dependence for sake of conciseness, we can define $\lambda ={\sigma }_x^2=f\left(\tau \right)$, which can be considered as a random variable itself distributed according to the PDF $P\left(\lambda \right)=p_{\tau }\left(f^{-1}\left(\lambda \right)\right)·{\partial }_{\lambda }\left(f^{-1}\left(\lambda \right)\right)$. The average over $\lambda$ is thus equivalent to computation of the expectation $\langle f\left(\tau \right)\rangle$ with respect to $\tau$.

In principle we may compute the expectation after the integration of Equation (27), however, it is much easier to compute it before performing the integration to avoid self-canceling terms:

$$\langle \lambda \rangle =2{\int }_0^t(t-s){\langle R(s,\tau )\rangle }_{\tau }ds,$$

(28)

For a generic PDF $p_{\tau }\left(\tau \right)$ we obtain:

$${\langle R(t,\tau )\rangle }_{\tau }=\nu {\int }_0^{\infty }\tau e^{-t/\tau }{p}_{\tau }\left(\tau \right)d\tau .$$

(29)

This expression is finite for any value of time only if $\langle \tau \rangle$ is finite. Moreover, this is a very important physical condition. In fact, when times goes to zero, ${\langle R(t=0,\tau )\rangle }_{\tau }$ is proportional to the average kinetic energy of the system.

The distribution $p_{\tau }\left(\tau \right)$ should have a power-law tail to introduce the desired anomalous time scaling of $\lambda$ but a finite value of the first moment of $\tau$ to maintain CLT applicability. The importance of this assumption can be seen explicitly by solving the integral in Equation (29) for the distribution employed in [1]:

$$p_{\tau }\left(\tau \right)=\frac{\alpha }{\mathsf{\Gamma }(1/\alpha )}\frac{1}{\tau }{L}_{\alpha }^{\alpha }\left(\frac{\tau }{C}\right),$$

(30)

where the constant $C=\langle \tau \rangle \frac{\mathsf{\Gamma }(1/\alpha )}{\alpha }$ serves to control the value of the average and $L_{\alpha }^{\alpha }(·)$ is the extreme Levy density distribution [10].

By considering the integral representation of the extremal Lévy density distribution and the Euler’s gamma function with some more simplification (see Section 3.5, equation 3.109 in [32]), the result in (29) can be represented by the integral form:

$$R\left(t\right)=\nu \langle \tau \rangle \frac{1}{2\pi i}{\int }_{\gamma -i\infty }^{\gamma +i\infty }\frac{\mathsf{\Gamma }(z/\alpha +1)\mathsf{\Gamma }(-z)}{\mathsf{\Gamma }(z+1)}{\left(\frac{t}{C}\right)}^{z}dz.$$

(31)

This expression can be solved through the residues theorem considering the poles $z/\alpha +1=-n$ or $z=n$, with $n=0,1,2,\dots ,\infty$, to obtain the short or the long time scaling of the variable. An interested reader can verify the explicit derivation in [1,32]. By plugging this result in Equation (28), without any assumption about time values, we observe that the condition of finite $\langle \tau \rangle$ is necessary to guarantee $\langle \lambda \rangle$ to be finite too, ensuring the Gaussian form of the PDF.

5. Discussions

The CLT has a fundamental role in the HEBP approach and, generally, in the theory of stochastic processes. The domain of attraction of the distribution of the increments determines the shape of the PDF of the stochastic process in the long time limit. In this work, we reviewed the main conditions of the classical CLT, by including also the less known case of distributions with infinite variance which fall in the Gaussian domain (with slower convergence). We proposed and analyzed a preparatory exercise to give the mathematical foundations to understand the approach used in HEBP to generate PDFs with Gaussian shape and non-linear scaling of the variance in time for the long time limit. It is shown that the sum of such a population of Gaussian random variables is mathematically defined by the sum of a more complex, and in general non-Gaussian, i.i.d. random variables. The population of Gaussian distributions can be interpreted, within a Bayesian approach, as the likelihood modulated by the prior distribution of a parameter of the model. The formal randomization of the parameter of the distribution (Equation (19)) is equivalent to the computation of the marginal likelihood, which corresponds indeed to the PDF of the i.i.d. random variables. This approach could be easily generalized to other distributions and parameters for statistical application purposes. The role of CLT in HEPB is then clarified. After recalling the derivation of the model, the conditions obtained in the preparatory example have been explicitly proven.