The Central Limit Theorem (CLT) is closely related to the “normal diffusion” (ND) process and is relevant to many random processes      . The characteristic of ND is the fact that the mean square displacement (MSD) of a diffusing particle scales with time as . The CLT together with the assumed zero mean and independence of increments implies that regardless of the microscopic model, the transition probability at large scales is well approximated by a Gaussian  . When this Gaussian kernel is convolved with an initial probability distribution, it evolves under time and after scaling tends to a Gaussian. We associate the evolution of the initial probability distribution in time with the semi-group generated by the Laplacian,  . Because of the spectral properties of the standard Laplacian, the time dependent (Gaussian) solution is an “attractor” solution to the ND equation. However, in recent years, great interest has been focused on the phenomena of “anomalous diffusion” (AD)   -  . As can be seen from the above references, such processes are observed for a host of physical, chemical and biological phenomena. The AD processes of interest to us here are those for systems with an “effective position dependent diffusion coefficient”, resulting in time evolution governed by a generalized Laplacian (e.g.,    ). Herein, we provide the derivation of a general, infinite family of such Laplacians, including new ones not included in the landmark study of O’Shaughnessy and Procaccia  . In the case of AD, the characteristic feature is the fact that the MSD scales with time as , where β is a real number. If , the process is termed “super diffusion” and if , the process is termed “sub-diffusion”. If , the diffusion is termed ballistic. O’Shaughnessy and Procaccia were the first to show that many AD systems have exact, attractor solutions, called “stretched Gaussians”:  . Here, D is a “constant effective diffusion coefficient”. We postulate the fact that these are computationally observed to be attractor solutions of the relevant differential equations implies that the standard CLT is “hidden” and applies also to AD associated with the diffusion equations we derive. In this paper, we show that this is indeed true.
The paper is organized as follows. In Section II, we summarize the solution of the ND equation, emphasizing its dependence on the Fourier transform (FT) of the probability function (the characteristic function) and noting that the characteristic function plays a key role in the proof of the CLT. In Section III, we discuss important connections among super symmetric quantum mechanics (SUSY), Heisenberg’s uncertainty principle (HUP), point transformations (PT) (in order to deduce generalizations of the FT and the related Laplacians), and the CLT for diffusion in the new canonically conjugate “position”. In Section IV, we explore the relations among AD, ND, scaling laws and PT’s. In Section V, we provide a computational example for the point transformation, , that illustrates the relation of the time dependence of the MSD to the functional form of the PT, W(x). Section VI discusses the detailed relationship between our diffusion equations and those in  . Finally, in Section VII, we summarize our results.
2. The Normal Diffusion Equation, Laplacian Semigroup and the Central Limit Theorem
In the case of ND as a random process, the proof of the CLT is most readily based on the characteristic function, i.e., the FT of the probability distribution  . At its heart is the fact that the Gaussian is invariant under the FT. This is interesting because the ND equation is also exactly solved using the FT of the diffusion equation. This makes use of the fact that the Laplacian satisfies
where is the x-representation of the momentum eigenket, and is the k-representation of the position eigenket. The FT kernel is therefore also an improper eigenstate of the standard Laplacian. This spectral relation ensures that the Laplacian generates a semigroup. It results in the “directional” restriction of solutions to the ND equation: the solution is stable only for increasing time. The semigroup does not possess a bounded inverse operator. It is also related to the attractor character of the solution to the ND equation  . The Gaussian minimizes the HUP for position and momentum and this provides a rigorous basis for the FT   . The diffusion equation,
where D is constant, becomes
under the FT . This is easily integrated (in
The inverse FT of a product of functions of k yields the exact convolution solution
Clearly, if is equal to 1, (the diffusing particle is localized initially at , with time dependence given by ). For any initial probability distribution whose FT is differentiable, the long-time behavior will be dominated by the overall Gaussian envelope (simply expand the square in the exponent in Equation (5) and factor out ). The in the normalization arises from requiring that the probability distribution be normalized to 1 under the measure dx. Then the MSD varies as (note, the slope of the MSD is not necessarily equal to one). We thus observe the role of the CLT as implying that the time dependence of any initial distribution, evolving according to Equation (2), will tend to a Gaussian distribution envelope after a sufficiently long time. We stress that the characteristic function proof of the CLT also rests on the facts that the Gaussian is invariant under the FT and the FT satisfies the convolution theorem. These properties are related to the semigroup structure of the standard Laplacian.
3. Super Symmetric Quantum Mechanics, Heisenberg’s Uncertainty Principle, Point Transformations, Canonical Quantization, Generalized Fourier Transforms and Laplacians
Recently, we have explored connections among super symmetric quantum mechanics (SUSY), Heisenberg’s uncertainty principle (HUP), point transformations (PT) and generalized Fourier transforms (GFT)    . We briefly summarize the relevant details here. A mathematically rigorous derivation of the FT (and its generalization) from the HUP has recently been given   . The fundamental starting point is that the minimum uncertainty state, , satisfies
The state cannot be an eigenstate simultaneously of both the position and momentum operators since they do not commute ( ). Explicitly, we note that Equation (6) reminds one of the property of simultaneous eigenvectors of commuting operator observables, :
The uncertainty product for such operators as has course, the absolute minimum value of zero and the simultaneous eigenvector is the only vector appearing in the above equations. The fact that the uncertainty product for position and momentum is positive definite results from the fact that neither nor can be simply proportional to the state . The best one can achieve for the two non-commuting operators is that their actions on the minimizing vector produce the same abstract vector, different from the minimizing state, and the new vector produced must have the same functional form in either the x- or k-representations. This immediately implies the form invariance of the Gaussian under the FT. In Appendix A, we show that this is a general property of the minimizing state vector for any pair of non-commuting, canonically conjugate Cartesian-like operators   .
The 1-D SUSY formalism is built on a program to make all systems (on the domain −∞ < x < ∞) look as much as possible like the harmonic oscillator, whose ground state also satisfies Equation (6). In the ladder operator notation, this equation is
In the coordinate representation, is .
In 1-D SUSY, the ground state is restricted to the form
for W’s such that the ground state is . Choosing the zero of energy to be that of the ground state, one finds that the potential of the system is given by
and in terms of the position representation, the ground state satisfies
As noted previously, this immediately implies that the ground state, , minimizes the uncertainty product ( ) and suggests that one ought to be able to obtain a new transform under which will be form invariant    . Unfortunately, since the commutator of W and is proportional to , which is constant only for W = x + constant, the desired transform kernel cannot be the simple FT kernel. As a consequence, we do not actually use the SUSY expressions directly to derive the desired transform. We remark that in the SUSY community, W is interpreted as a “super potential”, based on Equation (9). In our approach, we shall interpret W to be a generalized “position” variable, based the facts that: 1) minimizes and 2) the quintessential choice of W is that for the harmonic oscillator, where W is precisely the particle’s “generalized position”      .
Recently, in  , we explored choosing W to represent a point transformation of the usual position, such that 1) the domain of W is , (extension to the half line can be done) 2) the domain of the canonically conjugate momentum is also , 3) the transformation is invertible and 4) both W and can be interpreted as Cartesian-like coordinates. One expression that satisfies the above is the polynomial 
in which all coefficients are positive semi-definite and each even-power coefficient is always less than the next higher odd power coefficient. This results in a monotonic increasing behavior that is invertible almost everywhere (if one requires that the coefficient of the first power of x is positive definite, then the expression is invertible everywhere).
Assuming that W(x) satisfies the above conditions  , the classical canonically conjugate momentum, , is required to satisfy
where is the Poisson bracket with respect to the original Cartesian-like position and momentum   . It is easily shown that the classical canonical momentum can be expressed in the form
Invoking Occam’s razor, we set the integration constant along a constant-x integration path, g(x), equal to zero. Since above satisfies Equation (15) for all α, we then apply Dirac’s canonical quantization to obtain the generalized position and momentum operators (in the W-representation). This quantization consists of replacing the observables by operators, the Poisson bracket by the commutator of the operators and the scalar 1 by times the identity operator  :
We stress that these operators are manifestly self adjoint in the W-representation, with the measure dW. In addition, the minimizer of is obviously the Gaussian, , since
As in the case of the original position and momentum operators, is invariant under the “W-Fourier transform (W-FT) kernel”:
This is analogous to the FT in terms of x and k. This transform satisfies the convolution theorem under the measure dW or dK. We note that in the W-representation, since is self adjoint, we have the four equivalent, exact expressions:
which are all obviously self adjoint under the measure dW. The above relations are at the heart of our approach. The transformation that preserves the W-Gaussian , Equation (17), satisfies
Again, the W-Laplacian possesses a negative semi-definite spectrum so it is evident that this transform kernel possesses all the properties of the FT, including the fact that it supports a semigroup property.
We also note that there exists a generalized W-harmonic oscillator, with ground state satisfying  
The diffusion equation describing independent random motion in W is obviously
where here, D is a constant diffusion coefficient. It follows that the exact solution of Equation (27) is obtained using the W-FT kernel Equation (21), yielding
where is any proper, initial distribution.
The “characteristic function” underlying the above expression is
If , then equals one and the distribution is a Gaussian centered at . Note that in general, involves a positive function of time times a Gaussian envelope, . The analysis of the CLT for probability distributions in W is identical to that for the normal probability distribution  .
We next consider the situation where we quantize Equations (23)-(24) in the x-representation (we do not invoke the chain rule). Because of the ambiguity of operator ordering, we consider the following explicit expressions for a generalized momentum operator and its adjoint, involving the parameter α:
But these operators are obviously not self adjoint under the measure dx (except in the case where α = 1/2)! Normally, they are discarded and replaced by their arithmetic average.
However, we know that any self adjoint operator remains such, regardless of the representation. Of course, the point is that the above are self adjoint in the x-representation but the measure now includes the effect of the point transformation. Thus, Equation (30) is self adjoint under the measure and Equation (31) is self adjoint under the measure . When alpha equals 1/2, the measure for self adjointness reduces to simply dx.
We remark that, as shown above in Equation (23), in the W-domain, all four Laplacians possess the W-Gaussian as the ground state solution of the W-HO. Also, all four Laplacians satisfy Equation (25) in the W-domain. Thus, all four Laplacians are not only self adjoint in the x-representation but they also are all negative semi-definite operators which have an underlying semigroup structure regardless of their explicit form in the x-domain. Each of the four Laplacians will be self adjoint in the x-domain but only with appropriate measures. We now explicitly explore the forms of these Laplacians in the x-domain. It is convenient to group them in two pairs of two. This is because the forms
are already manifestly self adjoint under the measure dx. The other pair is
The above two Laplacians, and , require different measures from one another as well as from the first two Laplacians (except for ). We stress that both are self adjoint operators in the x-domain, with the correct measure for being and the correct measure for being . We next note that all four of these Laplacians support HO-like Hamiltonians. We define
We note that the first two Hamiltonians are self adjoint under the measure dx. The next two are self adjoint under the measures discussed above for Laplacians 3 and 4. It is easily seen that Hamiltonians 1 and 3 will have identical ground states, satisfying
Hamiltonians 2 and 4 will also have identical ground states, satisfying
These follow from the fact that the ground states are annihilated by the right hand operator in the respective factored Hamiltonians. We stress that, in the general case, none are simple Gaussians. Rather, the influence of the point transformation is clearly evident. Of course, if one sets or 1, then we see that the “stretched Gaussian” is one of the ground states. However, there is a second (bi-orthogonal partner) ground state which, in general, will be a multi-modal distribution. Because of the monotonic increasing character of the PT, is never negative and the ground state is always positive semi-definite. This is new and it reflects the fact that our formalism is closely related to the “Coupled SUSY” formalism introduced elsewhere   .
The next question of importance is: What are the eigenstates of the four x-domain Laplacians? These will provide the transformations under which the various generalized HO ground states are invariant (in analogy to the FT and the Gaussian distribution). Of course, in the “W-world”, there is only the W-FT, resulting from the W-Laplacian, having only the simple W-Gaussian as the solution of the diffusion and HO equations. We begin by considering the first two Laplacians. In this case, the explicit, analytical transforms have been derived in  for monomial W’s and α = 0. In that case, the first equation is
where (the Laplacian is negative semi-definite and self adjoint under the measure dx). The unitary transformation was derived restricting W to the monomial form . The result is
which is unitary under the measure dx. The stretched Gaussian, is invariant under this transform. The second transform equation is
and the unitary transform (under the measure dx), , is given by
In this case, it is the, in general, multi-modal distribution which is invariant.
The transformations for the second two Laplacians are easily obtained for general values of . For the Laplacian , the result is
The result for the Laplacian is given by
It is easily seen that these are bi-orthogonally related, since
and that (in Dirac notation),
Thus, the bi-orthogonal transforms serve to carry out Fourier-like transformations of the ground states of the GHOs. The ground state solution, , when transformed using under the measure dx, will produce the stretched Gaussian. The same is true when is transformed by under the measure dx. This explicitly shows how the underlying CLT for the W-domain Laplacian and W-Gaussian will be manifested in the x-representation.
4. Anomalous Diffusion and Normal Diffusion
We now discuss the implications of the above for diffusion. We consider the four linear diffusion equations
In the cases , the Laplacians are self adjoint under dx. When we choose monomial W’s with and, transformation under leads to
When we choose the same alpha and j = 2, transformation under again yields
In the cases of j = 3 and 4, we take care to recall the biorthogonal character of the transforms but it again results in equations of the same exact form as above, but now for all α. Once the time dependent equation is in the K-domain, the solution can be obtained using the original W-domain transform, Equation (17). (We note, however, that the transform of is not the same as that resulting from the inverses of). This transformation possesses a convolution theorem and in the W-domain, all behave in time like normal diffusion! We conclude that all of the linear diffusion equations derived herein are subject to the CLT associated with Gaussian distributions. The above analysis holds in general for any proper choice of PT but the above analysis for the j = 1, 2 cases has only been explicitly realized for the monomial case with alpha equal to zero. In the case of j = 3, 4, everything applies for general point PT’s, including polynomials  .
For simplicity, we shall illustrate these ideas using the specific example of the j = 3 case for alpha equal to zero. The transform equation is then
The diffusion equation is
The corresponding x-domain GFT kernel from Equation (21) is given by
When applied in the form Equation (56), the measure is. We note
that requiring K to have the same functional dependence on k as W(x) on x ensures the argument of the exponential is dimensionless and that the measure for
integration over K in the k-domain is. After applying the GFT
and using Equation (25), the exact solution (in the k-domain) is
The inverse GFT (which satisfies the convolution theorem in the variable W) then leads to the probability distribution at time t:
We stress that as derived, this is a probability distribution in the variable W(x), not the variable x. The measure for normalization will be, the normalization constant will be and the MSD will be proportional to. The diffusion is normal in the variable W.
However, Equation (58) can also be interpreted as a probability distribution in the variable x. It is a positive definite function of x and can be normalized to 1 under the measure dx:
If, e.g., , the integration will be transformed to the variable and the normalization will be proportional to. The average value of x will be zero (for symmetric initial distributions) and the MSD will scale as
This is, of course, the scaling for AD: the probability distribution is interpreted as applying to the physical position variable, x, rather than to the canonical variable,. The range corresponds to super diffusion, the range corresponds to ballistic diffusion and corresponds to sub-diffusion. However, because the Laplacian satisfies Equation (25), the time evolution will still be such that the CLT is satisfied. The normalization has nothing to do with the fact that the time evolution is governed by the attractor, whether we express it in the W- or the x-domain. Similar scaling results apply in the general cases.
5. A Computational Example for the Polynomial
We illustrate our results by the numerical example of a polynomial choice of PT:
We assume the initial distribution to be
This is very localized at W = 0, implying that it is centered at x = 0. The diffusion equation expressed in the x-domain is given by
We note that in the region, W is, to a very good approximation, equal to x but as the particle diffuses to larger magnitude values of x, it quickly transitions to being dominated by the cubic x-dependence. The W-domain results are shown in Figure 1 and Figure 2, where it is obvious that the diffusion is normal, since Equation (63) becomes Equation (27).
It is clear that at the earliest time, the linear x-dependence dominates the Laplacian and the diffusion is essentially normal, with the normal MSD t-scaling becoming more and more accurate as t approaches zero. At longer times, the cubic x-dependence dominates the Laplacian and we see that the MSD behavior tends to that of, just as one expects. While for more complicated polynomial PTs, such analysis will be more difficult, these results suggest that experimental
Figure 1. Time evolution in W coordinate is normal. The different color graphs correspond to times of diffusion, as indicated in the inset.
Figure 2. MSD in W coordinate. Diffusion is normal. The solid line is the result of solving the W-domain diffusion equation. The dashed line is the curve resulting from the MSD scaling as the time, t.
measurements of the scaling as a function of time can provide insight into the appropriate PT for which the diffusion will be normal. More important, this will
Figure 3. Evolution in x coordinate is anomalous. The different color graphs correspond to times of diffusion, as indicated in the inset.
Figure 4. Short and long time behavior of MSD in x. The early time shows ND and the later time AD. The red dashed curve beginning at t = 0 is the linear in time scaling law and the dashed blue curve at large time is the 1/3 power of t scaling law.
give information as to the best “effective” displacement variable with which to characterize and interpret the diffusion dynamics.
6. Relationship between Our Diffusion Equations and the O’Shaughnessy-Procaccia Equations
As formulated, our analysis and diffusion equations are very general and applicable both to fractal and non-fractal AD processes. Furthermore, for specific choices of our parameters, our equations exactly incorporate the radial diffusion equations of  . This might appear somewhat strange since our equations are defined on the entire real line while those of O’Shaughnessy and Procaccia are for radial variables on the half line. In fact, it is only a subset of our equations that incorporate the O’Shaughnessy and Procaccia equations  . The Laplacians are, explicitly, Equation (32) for α = zero. They were first derived on the half line and then extended to the full real axis and these are the particular equations that capture those of reference  . We explicitly prove this below. The equations of interest are given by
To transform this to our form of equations, we define the auxilliary variable
in the equation above. It easily follows that
Here, c is the dimension (non-integer c’s correspond to fractal cases). If one next defines the parameter
then Equation (66) becomes
The parameters, remain independent of one another. Of course, this is exactly our diffusion equation for the Laplacian in Equation (32) with. (In addition, we note that Equation (33) also captures the O’Shaughnessy and Procaccia equations for the specific value.) It is important to note, however, that even for the simplest monomial cases, with, all four of our Laplacians in Equations (32)-(35) differ from those of reference  . These new diffusion equations are also currently under computational study in our group.
First, we note our analysis shows that for diffusion according to any of the four diffusion equations, AD will be observed in the variable x but ND will be observed in the variable W. For ND to be manifested, one must analyze the results using the relevant generalized coordinate. AD as a function of x is simply “disguised ND” in W.
Second, we have shown that when linear anomalous diffusion is analyzed in terms of the appropriate PT displacement variable, the usual CLT applies. Of course, this in no way alters the fact that the situation is much more complex in the case of nonlinear diffusion    .
Third, our results show that, for such AD systems, experimental determination of the anomalous scaling leads directly to the identification of the relevant generalized coordinate, W(x) (in the case of the monomial choice of W(x)). Certainly, for the example case of, it is also quite easy to extract the relevant PT from the numerical data. For more general polynomial PTs, it will, naturally, be more complicated but the experimental scaling should still give information as to the relevant PT and effective displacement variable.
Fourth, we recall that the CLT is related to approximating the semigroup generated by the Laplacian operator. Because of this, we expect that any diffusion process for which a “proper” Laplacian operator exists will have an attractor solution that is invariant under the generalized transform and is an eigenfunction of the corresponding HO. It will be of interest to explore treatment of Levy processes using fractional powers of the generalized Laplacians.
Fifth, we have obtained linear diffusion equations that not only capture all those of reference  , but include infinite families of new equations for the values. It will be of interest to explore whether any experimental data correspond to these new linear diffusion equations.
Sixth, in our analysis, we considered diffusion in a 1D “Cartesian system” with a constant diffusion coefficient, D. It readily generalizes to any number of Cartesian random variables by simply summing the generalized Laplacians for each degree of freedom. We also point out that the scaling need not be the same in each degree of freedom. This will be important for anisotropic diffusion processes. It is also possible to introduce non-Cartesian coordinates (e.g., in a 3-D Cartesian system, one can define, along with angular variables to obtain a “spherically symmetric diffusion operator”). As in quantum mechanics, these should be derived by a coordinate transformation of the Cartesian-like Laplacians.
Finally, we are currently exploring the more general and difficult case of non-linear diffusion equations using our generalized Fourier transform methods.
This research was supported by the R. A. Welch Foundation under Grant E-0608. The author, DJK, gratefully acknowledges J. R. Klauder and G. H. Gunaratne for helpful conversations. The authors declare that no conflict of interests exists with the results and conclusions presented in this paper. Publication ethics have been observed.
In this Appendix we show the simple result that for any pair of non-commuting, canonically conjugate Cartesian-like variables, such as, their quantum operators generate a Fourier-like transform for which the minimizing quantum state is a Gaussian that is invariant under the transform. Since the variables are assumed to be Cartesian-like and canonically conjugate, they result in the commutator given in Equation (17). The minimizing condition is
Inserting the resolutions of the identity in the W- and -representations, one obtains
If we project this onto an eigenket of, we obtain
Projecting Equation (A.1) onto results in the analogous equation
By Equations (21)-(22), Equations (A.3) and (A.4) result in the W-wavefunction and the K-wavefunction being Gaussians of the exact same form.