Point Transformations and Relationships among Linear Anomalous Diffusion, Normal Diffusion and the Central Limit Theorem

Donald Kouri^1,2*, Nikhil Pandya^1,2, Cameron L. Williams², Bernhard G. Bodmann², Jie Yao³

¹Department of Physics, University of Houston, Houston, Texas, USA.
²Department of Mathematics, University of Houston, Houston, Texas, USA.
³Department of Mechanical Engineering, Texas Tech University, Lubbock, Texas, USA.

DOI: 10.4236/am.2018.92013                

Abstract

We present new connections among linear anomalous diffusion (AD), normal diffusion (ND) and the Central Limit Theorem (CLT). This is done by defining a point transformation to a new position variable, which we postulate to be Cartesian, motivated by considerations from super-symmetric quantum mechanics. Canonically quantizing in the new position and momentum variables according to Dirac gives rise to generalized negative semi-definite and self-adjoint Laplacian operators. These lead to new generalized Fourier transformations and associated probability distributions, which are form invariant under the corresponding transform. The new Laplacians also lead us to generalized diffusion equations, which imply a connection to the CLT. We show that the derived diffusion equations capture all of the Fractal and Non-Fractal Anomalous Diffusion equations of O’Shaughnessy and Procaccia. However, we also obtain new equations that cannot (so far as we can tell) be expressed as examples of the O’Shaughnessy and Procaccia equations. The results show, in part, that experimentally measuring the diffusion scaling law can determine the point transformation (for monomial point transformations). We also show that AD in the original, physical position is actually ND when viewed in terms of displacements in an appropriately transformed position variable. We illustrate the ideas both analytically and with a detailed computational example for a non-trivial choice of point transformation. Finally, we summarize our results.

Keywords

Generalized Fourier Analysis, Normal Diffusion, Anomalous Diffusion, Point Transformations, Canonical Quantization, Super Symmetric Quantum Mechanics

1. Introduction

The Central Limit Theorem (CLT) is closely related to the “normal diffusion” (ND) process and is relevant to many random processes [1] [2] [3] [4] [5] . The characteristic of ND is the fact that the mean square displacement (MSD) of a diffusing particle scales with time as $t^1$ . 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 [6] . 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, $\frac{{\partial }^2}{\partial x^2}$ [7] . 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) [5] [8] - [35] . 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., [4] [5] [15] ). 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 [15] . In the case of AD, the characteristic feature is the fact that the MSD scales with time as $t^{1/\beta }$ , where β is a real number. If $0.25\le \beta <1$ , the process is termed “super diffusion” and if $1<\beta$ , the process is termed “sub-diffusion”. If $0<\beta <0.25$ , 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”: $\rho \propto \exp \left[-x^{2+\vartheta }/4Dt\right]$ [15] . 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, $W\left(x\right)=x+x^3$ , 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 [15] . 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 [24] . 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

$\frac{{\text{d}}^2}{\text{d}{x}^2}{\text{e}}^{\pm ikx}=-k^2{\text{e}}^{\pm ikx}$ (1)

where ${\text{e}}^{ikx}$ is the x-representation of the momentum eigenket, $k^2\ge 0$ and ${\text{e}}^{-ikx}$ 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 [7] . The Gaussian minimizes the HUP for position and momentum and this provides a rigorous basis for the FT [36] [37] . The diffusion equation,

$\frac{\partial }{\partial t}\rho \left(x,t\right)=D\frac{{\partial }^2}{\partial x^2}\rho \left(x,t\right),$ (2)

where D is constant, becomes

$\frac{\partial }{\partial t}\hat{\rho }\left(k,t\right)=-k^{2}D\hat{\rho }\left(k,t\right),$ (3)

under the FT $\left(\hat{f}\left(k\right)=\frac{1}{\sqrt{\text{2π}}}{\displaystyle {\int }_{-\infty }^{\infty }\text{d}x\exp \left(-ikx\right)f\left(x\right)}\right)$ . This is easily integrated (in

time), yielding

$\hat{\rho }\left(k,t\right)=\exp \left(-D{k}^{2}t\right)\hat{\rho }\left(k,0\right).$ (4)

The inverse FT of a product of functions of k yields the exact convolution solution

$\rho \left(x,t\right)={\left[4\text{π}Dt\right]}^{-1/2}{\displaystyle {\int }_{-\infty }^{\infty }\text{d}{x}^{\prime }}\exp \left[-{\left(x-x^{\prime }\right)}^2/4Dt\right]\rho \left(x^{\prime },0\right).$ (5)

Clearly, if $\hat{\rho }\left(k,0\right)$ is equal to 1, $\rho \left(x^{\prime },0\right)=\delta \left(x^{\prime }-x_0\right)$ (the diffusing particle is localized initially at $x=x_0$ , with time dependence given by $\rho \left(x,t\right)=C\exp \left[-{\left(x-x_0\right)}^2/4Dt\right]$ ). For any initial probability distribution whose FT is differentiable, the long-time behavior will be dominated by the overall Gaussian envelope $\exp \left(-x^2/4Dt\right)$ (simply expand the square in the exponent in Equation (5) and factor out $\exp \left(-x^2/4Dt\right)$ ). The $t^{-1/2}$ in the normalization arises from requiring that the probability distribution be normalized to 1 under the measure dx. Then the MSD varies as $2Dt$ (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) [36] [37] [38] . We briefly summarize the relevant details here. A mathematically rigorous derivation of the FT (and its generalization) from the HUP has recently been given [36] [37] . The fundamental starting point is that the minimum uncertainty state, $|{\psi }_{\min }\rangle$ , satisfies

$\hat{x}|{\psi }_{\min }\rangle =-i\hat{k}|{\psi }_{\min }\rangle .$ (6)

The state $|{\psi }_{\min }\rangle$ cannot be an eigenstate simultaneously of both the position and momentum operators since they do not commute ( $\left[\hat{x},\hat{k}\right]=i\hat{1}$ ). Explicitly, we note that Equation (6) reminds one of the property of simultaneous eigenvectors of commuting operator observables, $\hat{A},\hat{B}$ :

$\hat{A}|a,b\rangle =a|a,b\rangle ,$ (7)

$\hat{B}|a,b\rangle =b|a,b\rangle ,$ (8)

$\Rightarrow \hat{A}|a,b\rangle =\frac{a}{b}\hat{B}|a,b\rangle .$ (9)

The uncertainty product for such operators as $\hat{A},\hat{B}$ 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 $\hat{x}|{\psi }_{\min }\rangle$ nor $-i\hat{k}|{\psi }_{\min }\rangle$ can be simply proportional to the state $|{\psi }_{\min }\rangle$ . 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 [39] [40] .

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

$\frac{1}{\sqrt{2}}\left[\hat{x}+i\hat{k}\right]|{\psi }_{\min }\rangle ={\hat{a}}_{HO}|{\psi }_{\min }\rangle =0\Rightarrow |{\psi }_{\min }\rangle =|{\psi }_0\rangle$ (10)

In the coordinate representation, ${\hat{a}}_{HO}$ is $\left(x+\text{d}/\text{d}x\right)/\sqrt{2}$ .

In 1-D SUSY, the ground state is restricted to the form

${\psi }_0\left(x\right)={\psi }_0\left(0\right)\exp \left[-{\displaystyle {\int }_0^x\text{d}{x}^{\prime }W\left(x^{\prime }\right)}\right]$ (11)

for W’s such that the ground state is $L^2$ . Choosing the zero of energy to be that of the ground state, one finds that the potential of the system is given by

$V=\left[W^2-\frac{\text{d}W}{\text{d}x}\right]/2,$ (12)

and in terms of the position representation, the ground state satisfies

$\frac{1}{\sqrt{2}}\left[W\left(x\right)+\text{d}/\text{d}x\right]{\psi }_0={\hat{a}}_W{\psi }_0=0.$ (13)

As noted previously, this immediately implies that the ground state, ${\psi }_0$ , minimizes the uncertainty product $\Delta W\Delta k$ ( $\hslash =1$ ) and suggests that one ought to be able to obtain a new transform under which ${\psi }_0$ will be form invariant [36] [37] [41] . Unfortunately, since the commutator of W and $\text{d}/\text{d}x$ is proportional to $\text{d}W/\text{d}x$ , 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) ${\psi }_0$ minimizes $\Delta W\Delta k$ and 2) the quintessential choice of W is that for the harmonic oscillator, where W is precisely the particle’s “generalized position” [36] [40] [41] [42] [43] .

Recently, in [37] , we explored choosing W to represent a point transformation of the usual position, such that 1) the domain of W is $-\infty <W<\infty$ , (extension to the half line can be done) 2) the domain of the canonically conjugate momentum is also $-\infty <P_W<\infty$ , 3) the transformation is invertible and 4) both W and $P_W$ can be interpreted as Cartesian-like coordinates. One expression that satisfies the above is the polynomial [37]

$W\left(x\right)={\displaystyle \underset{j=1}{\overset{2J+1}{\sum }}{a}_j}{x}^j,$ (14)

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 [37] , the classical canonically conjugate momentum, $P_W$ , is required to satisfy

$\left\{W,P_W\right\}=1,$ (15)

where $\left\{,\right\}$ is the Poisson bracket with respect to the original Cartesian-like position and momentum [37] [38] . It is easily shown that the classical canonical momentum can be expressed in the form

$P_W=\frac{1}{{\left(\frac{\text{d}W}{\text{d}x}\right)}^{1-\alpha }}{p}_x\frac{1}{{\left(\frac{\text{d}W}{\text{d}x}\right)}^{\alpha }}+g\left(x\right).$ (16)

Invoking Occam’s razor, we set the integration constant along a constant-x integration path, g(x), equal to zero. Since $P_W$ 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 $i\hslash$ times the identity operator [38] :

$\left[\hat{W},{\hat{P}}_W\right]=i\hat{1},\hslash =1,$ (17)

$\hat{W}=W,$ (18)

${\hat{P}}_W=-i\frac{\text{d}}{\text{d}W}.$ (19)

We stress that these operators are manifestly self adjoint in the W-representation, with the measure dW. In addition, the minimizer of $\Delta W\Delta P_W$ is obviously the Gaussian, $\exp \left(-W^2/2\right)$ , since

$W{\psi }_{\min }=-\frac{\text{d}{\psi }_{\min }}{\text{d}W}.$ (20)

As in the case of the original position and momentum operators, ${\psi }_{\min }$ is invariant under the “W-Fourier transform (W-FT) kernel”:

$\langle W|K\rangle =\exp \left(iKW\right)/\sqrt{\text{2π}},\langle K|W\rangle =\exp \left(-iKW\right)/\sqrt{\text{2π}},$ (21)

$-\infty <W<\infty ,-\infty <K<\infty .$ (22)

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 ${\hat{P}}_W$ is self adjoint, we have the four equivalent, exact expressions:

$-{\hat{P}}_W^{+}{\hat{P}}_W=-{\hat{P}}_W{\hat{P}}_W^{+}=-{\hat{P}}_W{\hat{P}}_W=-{\hat{P}}_W^{+}{\hat{P}}_W^{+}=\frac{{\text{d}}^2}{\text{d}{W}^2},$ (23)

${\hat{P}}_W={\hat{P}}_W^{+},$ (24)

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 $\left(\exp \left(-W^2/2\right)\right)$ , Equation (17), satisfies

$\frac{{\text{d}}^2}{\text{d}{W}^2}\exp \left(\pm iKW\right)=-K^2\exp \left(\pm iKW\right).$ (25)

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 [36] [37]

${\hat{a}}_W^{+}{\hat{a}}_W|{\psi }_0\rangle =\frac{1}{2}\left(-\frac{\text{d}}{\text{d}W}+W\right)\left(\frac{\text{d}}{\text{d}W}+W\right)\exp \left(-W^2/2\right)=0.$ (26)

The diffusion equation describing independent random motion in W is obviously

$\frac{\partial \rho }{\partial t}=D\frac{{\partial }^2\rho }{\partial W^2},$ (27)

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

$\rho \left(W,t\right)={\displaystyle {\int }_{-\infty }^{\infty }\text{d}{W}^{\prime }}\exp \left[-{\left(W-W^{\prime }\right)}^2/4Dt\right]\rho \left(W^{\prime },0\right),$ (28)

where $\rho \left(W^{\prime },0\right)$ is any proper, initial distribution.

The “characteristic function” underlying the above expression is

$\hat{\rho }\left(K,t\right)=\exp \left(-D{K}^{2}t\right)\hat{\rho }\left(K,0\right).$ (29)

If $\rho \left(W^{\prime },0\right)=\delta \left(W^{\prime }-W_0\right)$ , then $\hat{\rho }\left(K,0\right)$ equals one and the distribution is a Gaussian centered at $W_0$ . Note that in general, $\rho \left(W,t\right)$ involves a positive function of time times a Gaussian envelope, $\exp \left[-W^2/4Dt\right]$ . The analysis of the CLT for probability distributions in W is identical to that for the normal probability distribution [24] .

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 α:

${\hat{P}}_W=\frac{-i}{{\left(\frac{\text{d}W}{\text{d}x}\right)}^{1-\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\frac{\text{d}W}{\text{d}x}\right)}^{\alpha }},$ (30)

${\hat{P}}_W^{+}=\frac{-i}{{\left(\frac{\text{d}W}{\text{d}x}\right)}^{\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\frac{\text{d}W}{\text{d}x}\right)}^{1-\alpha }}.$ (31)

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 $\text{d}x{\left(\text{d}W/\text{d}x\right)}^{1-2\alpha }$ and Equation (31) is self adjoint under the measure $\text{d}x{\left(\text{d}W/\text{d}x\right)}^{2\alpha -1}$ . 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

${\Delta }_1=-{\hat{P}}_W^{+}{\hat{P}}_W=\frac{1}{{\left(\text{d}W/\partial x\right)}^{\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\text{d}W/\partial x\right)}^{2-2\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\text{d}W/\partial x\right)}^{\alpha }}$ (32)

and

${\Delta }_2=-{\hat{P}}_W{\hat{P}}_W^{+}=\frac{1}{{\left(\text{d}W/\partial x\right)}^{1-\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\text{d}W/\partial x\right)}^{2\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\text{d}W/\partial x\right)}^{1-\alpha }}$ (33)

are already manifestly self adjoint under the measure dx. The other pair is

${\Delta }_3=-{\hat{P}}_W{\hat{P}}_W=\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{1-\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{\left(\text{d}W/\text{d}x\right)}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{\alpha }}$ (34)

and

${\Delta }_4=-P_W^{+}{P}_W^{+}=\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{\left(\text{d}W/\text{d}x\right)}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{1-\alpha }}.$ (35)

The above two Laplacians, $-{\hat{P}}_W^{+}{\hat{P}}_W^{+}$ and $-{\hat{P}}_W{\hat{P}}_W$ , require different measures from one another as well as from the first two Laplacians (except for $\alpha =1/2$ ). We stress that both are self adjoint operators in the x-domain, with the correct measure for ${\Delta }_3$ being $\text{d}x{\left(\text{d}W/\text{d}x\right)}^{1-2\alpha }$ and the correct measure for ${\Delta }_4$ being $\text{d}x{\left(\text{d}W/\text{d}x\right)}^{2\alpha -1}$ . We next note that all four of these Laplacians support HO-like Hamiltonians. We define

${\hat{H}}_1=\frac{1}{2}\left(-\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{1-\alpha }}-W\right)\left(\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{1-\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{\alpha }}+W\right)$ (36)

${\hat{H}}_2=\frac{1}{2}\left(-\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{1-\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{\alpha }}+W\right)\left(\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{1-\alpha }}+W\right)$ (37)

${\hat{H}}_3=\frac{1}{2}\left(-\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{1-\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{\alpha }}+W\right)\left(\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{1-\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{\alpha }}+W\right)$ (38)

${\hat{H}}_4=\frac{1}{2}\left(-\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{1-\alpha }}+W\right)\left(\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{1-\alpha }}+W\right)$ (39)

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

$\left[\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{1-\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{\alpha }}+W\right]{\left(\text{d}W/\text{d}x\right)}^{\alpha }\exp \left(-W^2/2\right)=0.$ (40)

Hamiltonians 2 and 4 will also have identical ground states, satisfying

$\left[\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{\alpha }}\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\text{d}W/\text{d}x\right)}^{1-\alpha }}+W\right]{\left(\text{d}W/\text{d}x\right)}^{1-\alpha }\exp \left(-W^2/2\right)=0.$ (41)

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 $\alpha =0$ 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, $\text{d}W/\text{d}x$ 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 [42] [43] .

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 [36] for monomial W’s and α = 0. In that case, the first equation is

$-{\hat{P}}_W^{+}{\hat{P}}_W{\Phi }_K=\frac{\text{d}}{\text{d}x}\frac{1}{{\left(\text{d}W/\text{d}x\right)}^2}\frac{\text{d}{\Phi }_K}{\text{d}x}=-K^2{\Phi }_K$ (42)

where $K^2\ge 0$ (the Laplacian is negative semi-definite and self adjoint under the measure dx). The unitary transformation was derived restricting W to the monomial form $W\left(x\right)={\left|x\right|}^{\beta }$ . The result is

${\Phi }_K=\frac{1}{2}{\left|kx\right|}^{\beta -1/2}\left[J_{-1+1/2\beta }\left({\left|kx\right|}^{\beta }\right)-i\mathrm{sgn}\left({\left[kx\right]}^{\beta }\right)J_{1-1/2\beta }\left({\left|kx\right|}^{\beta }\right)\right]$ (43)

which is unitary under the measure dx. The stretched Gaussian, $\exp \left(-{\left|x\right|}^{2\beta }/2\right)$ is invariant under this transform. The second transform equation is

$\frac{1}{\left(\text{d}W/\text{d}x\right)}\frac{{\text{d}}^2}{\text{d}{x}^2}\frac{1}{\left(\text{d}W/\text{d}x\right)}{\tilde{\Phi }}_K=-K^2{\tilde{\Phi }}_K$ (44)

and the unitary transform (under the measure dx), ${\tilde{\Phi }}_K$ , is given by

${\tilde{\Phi }}_K=\frac{1}{2}\left[\mathrm{sgn}{\left(\eta \right)}^{\beta -1}{\left|\eta \right|}^{\left(\beta -1\right)/2}{J}_{-1/2\beta }\left({\left|\eta \right|}^{\beta }\right)+i\mathrm{sgn}{\left(\eta \right)}^{\beta }{\left|\eta \right|}^{\left(\beta -1\right)/2}{J}_{1/2\beta }\left({\left|\eta \right|}^{\beta }\right)\right]$ (45)

where $\eta =kx$ .

In this case, it is the, in general, multi-modal distribution $\left(\text{d}W/\text{d}x\right)\exp \left(-W^2/2\right)$ which is invariant.

The transformations for the second two Laplacians are easily obtained for general values of $\alpha$ . For the Laplacian ${\Delta }_3$ , the result is

${\Delta }_3{\phi }_K=-K^2{\phi }_K,{\phi }_K={\left(\text{d}W/\text{d}x\right)}^{\alpha }\exp \left(\pm iKW\right)/\sqrt{\text{2π}}.$ (46)

The result for the Laplacian ${\Delta }_4$ is given by

${\Delta }_4{\tilde{\phi }}_K=-K^2{\tilde{\phi }}_K,{\tilde{\phi }}_K={\left(\text{d}W/\text{d}x\right)}^{1-\alpha }\exp \left(\pm iKW\right)/\sqrt{\text{2π}}.$ (47)

It is easily seen that these are bi-orthogonally related, since

${\Delta }_3^{+}={\Delta }_4$ (48)

and that (in Dirac notation),

$\langle {\tilde{\phi }}_{K^{\prime }}|{\phi }_K\rangle =\langle {\phi }_K|{\tilde{\phi }}_{K^{\prime }}\rangle =\delta \left(K-K^{\prime }\right),$ (49)

$\hat{1}={\displaystyle {\int }_{-\infty }^{\infty }\text{d}K|{\phi }_K\rangle }\langle {\tilde{\phi }}_{K^{\prime }}|={\displaystyle {\int }_{-\infty }^{\infty }\text{d}K|{\tilde{\phi }}_K\rangle }\langle {\phi }_{K^{\prime }}|.$ (50)

Thus, the bi-orthogonal transforms serve to carry out Fourier-like transformations of the ground states of the GHOs. The ground state solution, ${\left(\text{d}W/\text{d}x\right)}^{\alpha }\exp \left(-W^2/2\right)$ , when transformed using ${\tilde{\phi }}_K$ under the measure dx, will produce the stretched Gaussian. The same is true when ${\left(\text{d}W/\text{d}x\right)}^{1-\alpha }\exp \left(-W^2/2\right)$ is transformed by ${\phi }_K$ 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

$\frac{\partial {\rho }_j}{\partial t}=D{\Delta }_j{\rho }_j,j=1\text{-}4.$ (51)

In the cases $j=1,2$ , the Laplacians are self adjoint under dx. When we choose monomial W’s with $\alpha =0$ and, transformation under leads to

(52)

When we choose the same alpha and j = 2, transformation under again yields

(53)

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 [37] .

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

(54)

The diffusion equation is

(55)

The corresponding x-domain GFT kernel from Equation (21) is given by

(56)

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

(57)

The inverse GFT (which satisfies the convolution theorem in the variable W) then leads to the probability distribution at time t:

(58)

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:

(59)

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

(60)

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:

(61)

We assume the initial distribution to be

(62)

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

(63)

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).

Of greater interest are the results where we interpret our density in the x-domain and use the Laplacian resulting in the diffusion Equation (63). The results are displayed in Figure 3 and Figure 4.

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 [15] . 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 [36] . 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 [15] . We explicitly prove this below. The equations of interest are given by

(64)

To transform this to our form of equations, we define the auxilliary variable

(65)

in the equation above. It easily follows that

(66)

Here, c is the dimension (non-integer c’s correspond to fractal cases). If one next defines the parameter

(67)

then Equation (66) becomes

(68)

where

(69)

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 [15] . These new diffusion equations are also currently under computational study in our group.

7. Conclusions

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 [4] [5] [6] .

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 [15] , 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.

Acknowledgements

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.

Appendix A

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

(A.1)

Inserting the resolutions of the identity in the W- and -representations, one obtains

(A.2)

If we project this onto an eigenket of, we obtain

(A.3)

Projecting Equation (A.1) onto results in the analogous equation

(A.4)

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.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Gillespie, D.T. and Seitaridou, E. (2013) Simple Brownian Motion. Oxford University Press, Cambridge, U.K.
[2]	Feller, W. (1971) An Introduction to Probability Theory and Its Applications. Vol. 2, Wiley, New York.
[3]	Harrison, J.M. (1971) Brownian Motion and Stochastic Flow Systems. Wiley, New York.
[4]	Tsallis, C. (2005) Nonextensive Statistical Mechanics, Anomalous Diffusion and Central Limit Theorems. Milan Journal of Mathematics, 73, 145-176. https://doi.org/10.1007/s00032-005-0041-1
[5]	Plastino, A. and Rocca, M.C. (2011) Inversion of Tsallis’ q-Fourier Transform and the Complex Plane Generalization. arXiv:1112.1985v1 [math-phys]
[6]	Einstein, A. (1905) Uber die von der molecular kinetischen Theorie der Warmegeforderte Bewegung von in ruhenden Flüssigkeitensuspendierten Teilchen. Ann. der Phys., 17, 549. https://doi.org/10.1002/andp.19053220806
[7]	Goldstein, J. (1985) Semigroups of Linear Operators and Applications. Oxford University Press, New York.
[8]	Tsallis, C. (2009) Introduction to Nonextensive Statistical Mechanics. Springer, New York.
[9]	Metzler, R. and Klafter, J. (2000) The Random Walk’s Guide to Anomalous Diffusion: A Fractional Dynamics Approach. Physics Reports, 339, 1-77. https://doi.org/10.1016/S0370-1573(00)00070-3
[10]	Metzler, R. and Klafter, J. (2004) The Restaurant at the End of the Random Walk: Recent Developments in the Description of Anomalous Transport by Fractional Dynamics. Journal of Physics A: Mathematical and General, 37, R161-R208. https://doi.org/10.1088/0305-4470/37/31/R01
[11]	Metzler, R., Jeon, J.H., Cherstvy, A.G. and Barkai, E. (2014) Anomalous Diffusion Models and Their Properties: Non-Stationarity, Non-Ergodicity, and Ageing at the Centenary of Single Particle Tracking. Physical Chemistry Chemical Physics, 16, 24128-24164. https://doi.org/10.1039/C4CP03465A
[12]	Ben-Avraham, D. and Havlin, S. (2000) Diffusion and Reactions in Fractals and Disordered Systems. Cambridge University Press, London. https://doi.org/10.1017/CBO9780511605826
[13]	Sornette, D. (2001) Critical Phenomena in Natural Sciences. Series in Synergetics, Springer, New York.
[14]	Risken, H. (1984) The Fokker-Planck Equation. Springer, New York, 63-95. https://doi.org/10.1007/978-3-642-96807-5_4
[15]	O’Shaughnessy, B. and Procaccia, I. (1985) Analytical Solutions for Diffusion on Fractal Objects. Physical Review Letters, 54, 455-458. https://doi.org/10.1103/PhysRevLett.54.455
[16]	Plerou, V., Gopikrishnan, P., Nunes Amaral, L.A., Gabaix, X. and Stanley, H.E. (2000) Economic Fluctuations and Anomalous Diffusion. Physical Review E, 62, R3023-R3026.
[17]	Barkai, E., Aghion, E. and Kessler, D.A. (2014) From the Area under the Bessel Excursion to Anomalous Diffusion of Cold Atoms. Physical Review X, 4, Article ID: 021036. https://doi.org/10.1103/PhysRevX.4.021036
[18]	Wang, K.G. (1992) Long-Time-Correlation Effects and Biased Anomalous Diffusion. Physical Review A, 45, 833-837. https://doi.org/10.1103/PhysRevA.45.833
[19]	Karger, J., Pfeifer, H. and Vojta, G. (1988) Time Correlation during Anomalous Diffusion in Fractal Systems and Signal Attenuation in NMR Field-Gradient Spectroscopy. Physical Review A, 37, 4514-4517. https://doi.org/10.1103/PhysRevA.37.4514
[20]	Grebenkov, D.S. (2007) NMR Survey of Reflected Brownian Motion. Reviews of Modern Physics, 79, 1077-1137. https://doi.org/10.1103/RevModPhys.79.1077
[21]	Gefen, Y., Aharony, A. and Alexander, S. (1983) Anomalous Diffusion on Percolating Clusters. Physical Review Letters, 50, 77-80. https://doi.org/10.1103/PhysRevLett.50.77
[22]	Zanette, D.H. and Alemany, P.A. (1995) Thermodynamics of Anomalous Diffusion. Physical Review Letters, 75, 366. https://doi.org/10.1103/PhysRevLett.75.366
[23]	Bohr, T. and Pikovsky, A. (1993) Anomalous Diffusion in the Kuramoto-Sivashinsky Equation. Physical Review Letters, 70, 2892-2895. https://doi.org/10.1103/PhysRevLett.70.2892
[24]	Kleuke, A. (2014) Probability Theory. Springer, New York.
[25]	Klafter, J. and Sokolov, I.M. (2005) Anomalous Diffusion Spreads Its Wings. Physics World, 18, 29-32. https://doi.org/10.1088/2058-7058/18/8/33
[26]	Sancho, J.M., Lacastra, A.M., Lindenberg, K., Sokolov, I.M. and Romero, A.H. (2004) Diffusion on a Solid Surface: Anomalous Is Normal. Physical Review Letters, 92, Article ID: 250601. https://doi.org/10.1103/PhysRevLett.92.250601
[27]	Lapas, L.C., Morgado, R., Vainstein, M.H., Rubi, J.M. and Oliveira, F.A. (2008) Khinchin Theorem and Anomalous Diffusion. Physical Review Letters, 101, Article ID: 230602. https://doi.org/10.1103/PhysRevLett.101.230602
[28]	Liu, B. and Goree, J. (2008) Superdiffusion and Non-Gaussian Statistics in a Driven-Dissipative 2D Dusty Plasma. Physical Review Letters, 100, Article ID: 055003. https://doi.org/10.1103/PhysRevLett.100.055003
[29]	He, Y., Burov, S., Metzler, R. and Barkai, E. (2008) Random Time-Scale Invariant Diffusion and Transport Coefficients. Physical Review Letters, 101, Article ID: 058101. https://doi.org/10.1103/PhysRevLett.101.058101
[30]	Amir, A., Oreg, Y. and Imry, Y. (2010) Localization, Anomalous Diffusion, and Slow Relaxations: A Random Distance Matrix Approach. Physical Review Letters, 105, Article ID: 070601. https://doi.org/10.1103/PhysRevLett.105.070601
[31]	Sagi, Y., Brook, M., Almog, I. and Davidson, N. (2012) Observation of Anomalous Diffusion and Fractional Self-Similarity in One Dimension. Physical Review Letters, 108, Article ID: 093002. https://doi.org/10.1103/PhysRevLett.108.093002
[32]	Hansen, Y.V., Gekle, S. and Netz, R.R. (2013) Anomalous Anisotropic Diffusion Dynamics of Hydration Water at Lipid Membranes. Physical Review Letters, 111, Article ID: 118103. https://doi.org/10.1103/PhysRevLett.111.118103
[33]	Agrawal, K., Gopalakrishnan, S., Knap, M., Müller, M. and Demler, E. (2015) Anomalous Diffusion and Griffiths Effects Near the Many-Body Localization Transition. Physical Review Letters, 114, Article ID: 160401. https://doi.org/10.1103/PhysRevLett.114.160401
[34]	Cisternas, J., Descalzi, O., Albers, T. and Radons, G. (2016) Anomalous Diffusion of Dissipative Solitons in the Cubic-Quintic Complex Ginzburg-Landau Equation in Two Spatial Dimensions. Physical Review Letters, 116, Article ID: 203901. https://doi.org/10.1103/PhysRevLett.116.203901
[35]	Golestanian, R. (2009) Anomalous Diffusion of Symmetric and Asymmetric Active Colloids. Physical Review Letters, 102, Article ID: 188305. https://doi.org/10.1103/PhysRevLett.102.188305
[36]	Williams, C.L., Bodmann, B.G. and Kouri, D.J. (2017) Fourier and Beyond: Invariance Properties of a Family of Integral Transforms. Journal of Fourier Analysis and Applications, 23, 660-678. https://doi.org/10.1007/s00041-016-9482-x
[37]	Kouri, D.J., Williams, C.L. and Pandya, N. (2017) Canonical Transformations, Quantization, Mutually Unbiased and Other Complete Bases. Applied Mathematics, 8, 901-919. https://doi.org/10.4236/am.2017.87071
[38]	Dirac, P.A.M. (1958) Quantum Mechanics. 4th Edition, Oxford U. P., London.
[39]	Junker, G. (1996) Supersymmetric Methods in Quantum and Statistical Physics. Springer, New York. https://doi.org/10.1007/978-3-642-61194-0
[40]	Klauder, J.R. (2015) Enhanced Quantization. World Scientific, Singapore. https://doi.org/10.1142/9452
[41]	Chou, C.-C., Biamonte, M.T., Bodmann, B.G. and Kouri, D.J.J. (2012) New System-Specific Coherent States for Bound State Calculations. Journal of Physics A: Mathematical and Theoretical, 45, Article ID: 505303. https://doi.org/10.1088/1751-8113/45/50/505302
[42]	Williams, C.L. (2017) Article Title. PhD Thesis, University of Houston, Department of Mathematics, Houston.
[43]	Williams, C.L., Pandya, N.N., Bodmann, B.G. and Kouri, D.J. (2017) Coupled Supersymmetry and Ladder Structures beyond the Harmonic Oscillator.