Fractional Calculus, Langevin Equation and Quaternions 
in Mathematica

Sadataka Furui, Teikyo University, Utsunomiya Japan

#WolframTechConf

ABSTRACT

In analysis of nonlinear differential equations, fractional order calculus is attracting interests of researchers.  Langevin equation is a stochastic differential equation, and It can be reduced to Ito’s stochastic differential equation.  We take a specific Wiener process and obtain the Fokker-Planck equation, which can be regarded as 1/2 order fractional differential equation. 
There are several Fractional order differential equations. We consider the Caputo fractional differential, with initial conditions and with boundary conditions. The left derivative and the right derivative can be calculated by Mathematica and they are different for f(x)=(x-a)^(n q)(b-x)^(n q),  a<x<b, where q is not integer.
In fractional differential equation, Complex variables can be extended to quaternion and biquaternion variables. 
In analysis of (2+1)D ultrasonic wave spectra,  a solution of quaternions differential equation reproduces the conformality.
In analysis of (3+1)D ultrasonic wave propagation, Tomonaga’s spacelike time based split-quaternions can be used for the fixed point action lattice simulation.

1. Introduction of Langevin equation

In 1997, Nakajima and I studied the Langevin equation 

where <...> is averaging over the product of Gaussian distribution and
a solution of the Ito’s stochastic differential equation. We introduced a parameter s=Sqrt[t] and for u^i=k x^i, obtained the corresponding Fokker-Planck equation P[t,k,η], where t is time and k is momentum.
The Fokker - Planck equation P[t, k, η] =(p1,p2,p3)  satisfies D[pi,t]=D[D[pi,x],x]+D[k x pi,x] (i=1,2,3).

1.1 Fokker - Planck equation

The Fokker - Planck equation can be solved exactly in the case of  u^i(x)=-k x. We introduce a parameter α[t, k].

α[t_,k_]:=(1-E^(-2k t))/(2k)

We choose k = 1/2, α[t, 1/2] = 1 - E^(-1)

Plot[α[t,1/2],{t,0.01,10},PlotRange->All,AxesLabel->{"t","α"}]

The initial condition x (0) = x01 yields :

X1t=x01 E^(-k t);
p1=E^((-1/(4α[t,k]))(x1-X1t)^2)/Sqrt[4Pi α[t,k] ]

We check the Fokker - Planck equation D[pi, t] = D[D[pi, x], x] + D[k x pi, x] (i = 1, 2, 3)

D[p1,t]

D[p1, t] is plotted as a function of t and x .

Plot3D[-(((E^(-((k (((x01-(E^(k t)) x1))^2))/(2 (-1+(E^(2 k t))))))) k (-1+(E^(k t)) k x01 x1+(E^(3 k t)) k x01 x1+(E^(2 k t)) (1-k ((x01^2)+(x1^2)))))/((((-1+(E^(2 k t))))^2) (Sqrt[((2-2 (E^(-2 k t)))/k)]) (Sqrt[π])))/.{k->0.5,x01->0.5},{t,0.01,0.1},{x1,0,2},AxesLabel->{"t","x"},PlotRange->All]

D[D[p1, x], x] and  D[k x p1, x] are plotted as a function of t and x .

D[D[p1,x1],x1]

Plot3D[-(((E^(-((k (((-(E^(-k t)) x01+x1))^2))/(2 (1-(E^(-2 k t))))))) k)/((1-(E^(-2 k t))) (Sqrt[((1-(E^(-2 k t)))/k)]) (Sqrt[2 π])))+(((E^(-((k (((-(E^(-k t)) x01+x1))^2))/(2 (1-(E^(-2 k t))))))) (k^2) (((-(E^(-k t)) x01+x1))^2))/((((1-(E^(-2 k t))))^2) (Sqrt[((1-(E^(-2 k t)))/k)]) (Sqrt[2 π])))/.{k->0.5,x01->0.5},{t,0.01,0.1},{x1,0,2},AxesLabel->{"t","x"},PlotRange->All,PlotStyle->Gray]

D[k x1 p1,x1]

Plot3D[(((E^(-((k (((-(E^(-k t)) x01+x1))^2))/(2 (1-(E^(-2 k t))))))) k)/((Sqrt[((1-(E^(-2 k t)))/k)]) (Sqrt[2 π])))-(((E^(-((k (((-(E^(-k t)) x01+x1))^2))/(2 (1-(E^(-2 k t))))))) (k^2) x1 (-(E^(-k t)) x01+x1))/((1-(E^(-2 k t))) (Sqrt[((1-(E^(-2 k t)))/k)]) (Sqrt[2 π])))/.{k->0.5,x01->0.5},{t,0.01,0.1},{x1,0,2},AxesLabel->{"t","x"},PlotRange->All,PlotStyle->Red]

Check the Fokker - Planck equation .

D[D[p1,x1],x1]+D[k x1 p1,x1] -D[p1,t]

Simplify[%]

p2, p3 also satisfy the equation. Exact solution of the Fokker-Planck equation is verified.

1.2 Langevin equation as fractional order differential equation

As the parameter of Brownian step s = Sqrt[t], where t is time, the Langevin equation is a Fractional order differential equation.    The Riemann - Liouville integral can be calculated by Mathematica.

In the flat 3 D space, the Langevin  equation is in the form
       D[Xi, t] = ui[X] + ei_α[X] η_α[t], (i=1,2,3).
Langevin step is given by (X1, X2, X3), where Xi is a function of s = Sqrt[t], (η1, η2, η3) are Gaussian Random variable,  and ui[X] is the input.

X1[s_,η1_,η2_,η3_,u1_,u2_,u3_]:=s η1+s^2 u1+(s^3/2)(η1 D[u1,x1]+η2 D[u1,x2]+η3 D[u1,x3])+s^4((1/6)(η1  η1 D[D[u1 ,x1],x1]+η1  η2 D[D[u1 ,x1],x2]+η1  η3 D[D[u1 ,x1],x3]+η2  η1 D[D[u1 ,x2],x1]+η2  η2 D[D[u1 ,x2],x2]+η2  η3 D[D[u1 ,x2],x3]+η3  η1 D[D[u1 ,x3],x1]+η3  η2 D[D[u1 ,x3],x2]+η3  η3 D[D[u1 ,x3],x3])+(1/2)(u1 D[u1,x1]+u2 D[u1,x2]+u3 D[u1,x3])+(1/6)Laplacian[u1,{x1,x2,x3}])+s^5((1/3)(η1 u1 D[D[u1,x1],x1]+η2 u1 D[D[u1,x2],x1]+η3 u1 D[D[u1,x3],x1] +η1 u2 D[D[u1,x1],x2]+η2 u2 D[D[u1,x2],x2]+η3 u2 D[D[u1,x3],x2]+η1 u3 D[D[u1,x1],x3]+η2 u3 D[D[u1,x2],x3]+η3 u3 D[D[u1,x3],x3])+(1/6)(η1 D[u1,x1]D[u1,x1]+η2 D[u1,x2]D[u1,x2]+η3 D[u1,x3]D[u1,x3]+η1 D[u2,x1]D[u1,x2]+η2 D[u2,x2]D[u1,x2]+η3 D[u2,x3]D[u1,x2]+η1 D[u3,x1]D[u1,x3]+η2 D[u3,x2]D[u1,x3]+η3 D[u3,x3]D[u1,x3])+(1/24)(η1 D[u1,x1]D[u1,x1]+η2 D[u2,x1]D[u1,x2]+η3 D[u3,x1]D[u1,x3]+η1 D[u1,x2]D[u1,x1]+η2 D[u2,x2]D[u1,x2]+η3 D[u3,x2]D[u1,x3]+η1 D[u1,x3]D[u1,x1]+η2 D[u2,x3]D[u1,x2]+η3 D[u3,x3]D[u1,x3])+(1/3)(η1 D[Laplacian[u1,{x1,x2,x3}],x1]+η2 D[Laplacian[u1,{x1,x2,x3}],x2]+η3 D[Laplacian[u1,{x1,x2,x3}],x3]))+s^6((1/3)(D[u1,x1]D[D[u1,x1],x1] +D[u1,x2]D[D[u1,x1],x2]+D[u1,x3]D[D[u1,x1],x3]+D[u2,x1]D[D[u1,x2],x1] +D[u2,x2]D[D[u1,x2],x2]+D[u2,x3]D[D[u1,x2],x3]+D[u3,x1]D[D[u1,x3],x1] +D[u3,x2]D[D[u1,x3],x2]+D[u3,x3]D[D[u1,x3],x3])+(1/6)(u1*u1 D[D[u1,x1],x1]+u1*u2 D[D[u1,x2],x1]+u1*u3 D[D[u1,x3],x1]+u2*u1 D[D[u1,x1],x2]+u2*u2 D[D[u1,x2],x2]+u2*u3 D[D[u1,x3],x2]+u3*u1 D[D[u1,x1],x3]+u3*u2 D[D[u1,x2],x3]+u3*u3 D[D[u1,x3],x3])+(1/6)(Laplacian[u1,{x1,x2,x3}]D[u1,x1]+Laplacian[u2,{x1,x2,x3}] D[u1,x2]+Laplacian[u3,{x1,x2,x3}]D[u1,x3])+(1/6)(u1 D[u1,x1]D[u1,x1]+u2 D[u1,x2]D[u1,x1]+u3 D[u1,x3]D[u1,x1]+u1 D[u2,x1]D[u1,x2]+u2 D[u2,x2]D[u1,x2]+u3 D[u2,x3]D[u1,x2]+u1 D[u3,x1]D[u1,x3]+u2 D[u3,x2]D[u1,x3]+u3 D[u3,x3]D[u1,x3])+(1/3)(u1 D[Laplacian[u1,{x1,x2,x3}],x1]+u2 D[Laplacian[u1,{x1,x2,x3}]+u3 D[Laplacian[u1,{x1,x2,x3}]]])+(1/6)Laplacian[Laplacian[u1,{x1,x2,x3}],{x1,x2,x3}]
);

X2[s_,η1_,η2_,η3_,u1_,u2_,u3_]:=s η2+s^2 u2+(s^3/2)(η1 D[u2,x1]+η2 D[u2,x2]+η3 D[u2,x3])+s^4((1/6)(η1  η1 D[D[u2 ,x1],x1]+η1  η2 D[D[u2 ,x1],x2]+η1  η3 D[D[u2 ,x1],x3]+η2  η1 D[D[u2 ,x2],x1]+η2  η2 D[D[u2 ,x2],x2]+η2  η3 D[D[u2 ,x2],x3]+η3  η1 D[D[u2 ,x3],x1]+η3  η2 D[D[u2 ,x3],x2]+η3  η3 D[D[u2 ,x3],x3])+(1/2)(u1 D[u2,x1]+u2 D[u2,x2]+u3 D[u2,x3])+(1/6)Laplacian[u2,{x1,x2,x3}])+s^5((1/3)(η1 u1 D[D[u2,x1],x1]+η2 u1 D[D[u2,x2],x1]+η3 u1 D[D[u2,x3],x1] +η1 u2 D[D[u2,x1],x2]+η2 u2 D[D[u2,x2],x2]+η3 u2 D[D[u2,x3],x2]+η1 u3 D[D[u2,x1],x3]+η2 u3 D[D[u2,x2],x3]+η3 u3 D[D[u2,x3],x3])+(1/6)(η1 D[u1,x1]D[u2,x1]+η2 D[u1,x2]D[u2,x2]+η3 D[u1,x3]D[u2,x3]+η1 D[u2,x1]D[u2,x2]+η2 D[u2,x2]D[u2,x2]+η3 D[u2,x3]D[u2,x2]+η1 D[u3,x1]D[u2,x3]+η2 D[u3,x2]D[u2,x3]+η3 D[u3,x3]D[u2,x3])+(1/24)(η1 D[u1,x1]D[u2,x1]+η2 D[u2,x1]D[u2,x2]+η3 D[u3,x1]D[u2,x3]+η1 D[u1,x2]D[u2,x1]+η2 D[u2,x2]D[u2,x2]+η3 D[u3,x2]D[u2,x3]+η1 D[u1,x3]D[u2,x1]+η2 D[u2,x3]D[u2,x2]+η3 D[u3,x3]D[u2,x3])+(1/3)(η1 D[Laplacian[u2,{x1,x2,x3}],x1]+η2 D[Laplacian[u2,{x1,x2,x3}],x2]+η3 D[Laplacian[u2,{x1,x2,x3}],x3]))+s^6((1/3)(D[u1,x1]D[D[u2,x1],x1] +D[u1,x2]D[D[u2,x1],x2]+D[u1,x3]D[D[u2,x1],x3]+D[u2,x1]D[D[u2,x2],x1] +D[u2,x2]D[D[u2,x2],x2]+D[u2,x3]D[D[u2,x2],x3]+D[u3,x1]D[D[u2,x3],x1] +D[u3,x2]D[D[u2,x3],x2]+D[u3,x3]D[D[u2,x3],x3])+(1/6)(u1*u1 D[D[u2,x1],x1]+u1*u2 D[D[u2,x2],x1]+u1*u3 D[D[u2,x3],x1]+u2*u1 D[D[u2,x1],x2]+u2*u2 D[D[u2,x2],x2]+u2*u3 D[D[u2,x3],x2]+u3*u1 D[D[u2,x1],x3]+u3*u2 D[D[u2,x2],x3]+u3*u3 D[D[u2,x3],x3])+(1/6)(Laplacian[u1,{x1,x2,x3}]D[u2,x1]+Laplacian[u2,{x1,x2,x3}] D[u2,x2]+Laplacian[u3,{x1,x2,x3}]D[u2,x3])+(1/6)(u1 D[u1,x1]D[u2,x1]+u2 D[u1,x2]D[u2,x1]+u3 D[u1,x3]D[u2,x1]+u1 D[u2,x1]D[u2,x2]+u2 D[u2,x2]D[u2,x2]+u3 D[u2,x3]D[u2,x2]+u1 D[u3,x1]D[u2,x3]+u2 D[u3,x2]D[u2,x3]+u3 D[u3,x3]D[u2,x3])+(1/3)(u1 D[Laplacian[u2,{x1,x2,x3}],x1]+u2 D[Laplacian[u2,{x1,x2,x3}]+u3 D[Laplacian[u2,{x1,x2,x3}]]])+(1/6)Laplacian[Laplacian[u2,{x1,x2,x3}],{x1,x2,x3}]
);

X3[s_,η1_,η2_,η3_,u1_,u2_,u3_]:=s η3+s^2 u3+(s^3/2)(η1 D[u3,x1]+η2 D[u3,x2]+η3 D[u3,x3])+s^4((1/6)(η1  η1 D[D[u3 ,x1],x1]+η1  η2 D[D[u3 ,x1],x2]+η1  η3 D[D[u3 ,x1],x3]+η2  η1 D[D[u3 ,x2],x1]+η2  η2 D[D[u3 ,x2],x2]+η2  η3 D[D[u3 ,x2],x3]+η3  η1 D[D[u3 ,x3],x1]+η3  η2 D[D[u3 ,x3],x2]+η3  η3 D[D[u3 ,x3],x3])+(1/2)(u1 D[u3,x1]+u2 D[u3,x2]+u3 D[u3,x3])+(1/6)Laplacian[u3,{x1,x2,x3}])+s^5((1/3)(η1 u1 D[D[u3,x1],x1]+η2 u1 D[D[u3,x2],x1]+η3 u1 D[D[u3,x3],x1] +η1 u2 D[D[u3,x1],x2]+η2 u2 D[D[u3,x2],x2]+η3 u2 D[D[u3,x3],x2]+η1 u3 D[D[u3,x1],x3]+η2 u3 D[D[u3,x2],x3]+η3 u3 D[D[u3,x3],x3])+(1/6)(η1 D[u1,x1]D[u3,x1]+η2 D[u1,x2]D[u3,x2]+η3 D[u1,x3]D[u3,x3]+η1 D[u2,x1]D[u3,x2]+η2 D[u2,x2]D[u3,x2]+η3 D[u2,x3]D[u3,x2]+η1 D[u3,x1]D[u3,x3]+η2 D[u3,x2]D[u3,x3]+η3 D[u3,x3]D[u3,x3])+(1/24)(η1 D[u1,x1]D[u3,x1]+η2 D[u2,x1]D[u3,x2]+η3 D[u3,x1]D[u3,x3]+η1 D[u1,x2]D[u3,x1]+η2 D[u2,x2]D[u3,x2]+η3 D[u3,x2]D[u3,x3]+η1 D[u1,x3]D[u3,x1]+η2 D[u2,x3]D[u3,x2]+η3 D[u3,x3]D[u3,x3])+(1/3)(η1 D[Laplacian[u3,{x1,x2,x3}],x1]+η2 D[Laplacian[u3,{x1,x2,x3}],x2]+η3 D[Laplacian[u3,{x1,x2,x3}],x3]))+s^6((1/3)(D[u1,x1]D[D[u3,x1],x1] +D[u1,x2]D[D[u3,x1],x2]+D[u1,x3]D[D[u3,x1],x3]+D[u2,x1]D[D[u3,x2],x1] +D[u2,x2]D[D[u3,x2],x2]+D[u2,x3]D[D[u3,x2],x3]+D[u3,x1]D[D[u3,x3],x1] +D[u3,x2]D[D[u3,x3],x2]+D[u3,x3]D[D[u3,x3],x3])+(1/6)(u1*u1 D[D[u3,x1],x1]+u1*u2 D[D[u3,x2],x1]+u1*u3 D[D[u3,x3],x1]+u2*u1 D[D[u3,x1],x2]+u2*u2 D[D[u3,x2],x2]+u2*u3 D[D[u3,x3],x2]+u3*u1 D[D[u3,x1],x3]+u3*u2 D[D[u3,x2],x3]+u3*u3 D[D[u3,x3],x3])+(1/6)(Laplacian[u1,{x1,x2,x3}]D[u3,x1]+Laplacian[u2,{x1,x2,x3}] D[u3,x2]+Laplacian[u3,{x1,x2,x3}]D[u3,x3])+(1/6)(u1 D[u1,x1]D[u3,x1]+u2 D[u1,x2]D[u3,x1]+u3 D[u1,x3]D[u3,x1]+u1 D[u2,x1]D[u3,x2]+u2 D[u2,x2]D[u3,x2]+u3 D[u2,x3]D[u3,x2]+u1 D[u3,x1]D[u3,x3]+u2 D[u3,x2]D[u3,x3]+u3 D[u3,x3]D[u3,x3])+(1/3)(u1 D[Laplacian[u3,{x1,x2,x3}],x1]+u2 D[Laplacian[u3,{x1,x2,x3}]+u3 D[Laplacian[u3,{x1,x2,x3}]]])+(1/6)Laplacian[Laplacian[u3,{x1,x2,x3}],{x1,x2,x3}]
);

The triality symmetry can be checked by the command Simplify .

Simplify[X3[Sqrt[t],η1,η2,η3,-k x1,-k x2,-k x3]-X2[Sqrt[t],η1,η2,η3,-k x1,-k x2,-k x3]]

Simplify[X3[Sqrt[t],η1,η2,η3,-k x1,-k x2,-k x3]-X1[Sqrt[t],η1,η2,η3,-k x1,-k x2,-k x3]]

Simplify[X2[Sqrt[t],η1,η2,η3,-k x1,-k x2,-k x3]-X1[Sqrt[t],η1,η2,η3,-k x1,-k x2,-k x3]]

Symmetry under interchanges of indices 1, 2, 3 can be checked by the command Simplify .

2. Caputo Fractional Integral

Kiskinov et al . (2024) discussed whether the Caputo Fractional derivative 

should be equal for continuous function f(x). 
He considered the integral of f[s_]:=s^0.25 (1-s)^0.25.  
The former is 1/Γ(1-0.25) times the following integral.

f[s_]:=s^0.25 (1-s)^0.25;
Integrate[(t-s)^(-0.25)D[f[s],s],{s,0,t}]

We choose t = 0.2 and obtain

1.1107207345395917` Hypergeometric2F1[-0.25`,0.25`,1.`,t]-(0.27768018363489794` +2.7755575615628914`*^-17 I) (t^(1.`)) Hypergeometric2F1[0.75`,1.25`,2.`,t]/.{t->0.2}

1.03457/Gamma[0.75]

The latter is 1/Γ (1 - 0.25) times the following integral .

Integrate[(-t+s)^(-0.25)D[f[s],s],{s,t,1}]

The condition is complicated, and I replace 1 to b .

Integrate[(-t+s)^(-0.25)D[f[s],s],{s,t,b}]

Choosing b = 1, I obtain

-(1/(((1.` -1.` t))^(1/4)))0.19999999999999998` (((1.` -1.` b))^(5/4)) AppellF1[1.25`,0.75`,0.25`,2.25`,1.` -1.` b,((1.` -1.` b)/(1.` -1.` t))]-(0.14142135623730953` +0.14142135623730953` I) b ((((b/t)))^(1/4)) AppellF1[1.25`,0.75`,0.25`,2.25`,b,(b/t)]+(0.2776801836348979` -0.2776801836348979` t) Hypergeometric2F1[0.75`,1.25`,2.`,1.` -1.` t]+(0.19634954084936207` +0.19634954084936207` I) t Hypergeometric2F1[0.75`,1.25`,2.`,t]/.{b->1,t->0.2}

-0.607571*(-1/Gamma[0.75])

Kiskinov et al . obtained the same value with opposite sign. I agree that the left side integral and the right side integral are not equal.

2.1 Gegenbauer Polynomial

Lane - Emden equation is a well - known fractional differential equation . 
   Ergindy and Refat (2018) proposed to solve the equation using the Gegenbauer Polynomial . Mathematica has GegenbauerC[n, m, z] that satisfy.
   (1 - x^2) y "-(2m+1) x y'+n(n+2m) y=0
The polynomial is real valued orthogonal function with respect to the weight function :

w[α_,x_]:=(1-x^2)^(α-1/2)

Orthogonality of Gegenbauer polynomial can be checked by using

dot[f_,g_]:=(2/Pi)Integrate[ Conjugate[f]g Sqrt[1-x^2],{x,-1,1}]

and Orthogonalize command .

Orthogonalize[{1,x,x^2,x^3,x^4,x^5,x^6},dot]//Expand

The polynomials agree with GegenbauerC[n, 1, x] .

Table[GegenbauerC[n,1,x],{n,0,6}]

Using  the expression xx1[Sqrt[t], η1, k, x1] obtained from X1[Sqrt[t], η1, η2, η3, -k x1, -k x2, -k x3],

xx1[st_,η1_,k_,x1_]:=-k st^2 x1+(k^2/2)st^4 x1-(k^3/6)st^6 x1+st η1-(k/2)st^3 η1+(5k^2/24)st^5 η1;

Choosing η1 = 0.2, k = 0.5, I obtain 6 th order polynomial of st :

xx1[st,0.2,0.5,x]

Each term can be expanded by Gegenbauer polynomials, e.g.
-0.0208333 st^6=-0.0016276041666666665`-0.0029296874999999996`*(-1+4 (st^2))-0.0016276041666666672`*(1-12 (st^2)+16 (st^4))-0.0003255208333443612`*(-1+24 (st^2)-80 (st^4)+64 (st^6))

2.2 ProductLog Function

Ergindy transformed the parameter t to t-τ =t y^(1/α) and considered branching cut by a function z=w e^w.

Reduce[E^(u+1)(u+1)==E/2,u,Complexes]

N[ProductLog[0,E/2]]-1

Ergindy considered Gegenbauer polynomial of real variables x, but it can be complex z .

GegenbauerC[5,z]

ComplexPlot3D[GegenbauerC[5,z],{z,-2-2I,2+2I},AxesLabel->{"Re z","Im z","GegenbauerC[5,z]"}]

3. Quaternion Differential Equation

Quaternion H is an extension of complex number  expressed as 
       q=q0 +q1 i+q2 j +q3 k ,  i j=-j i=k , j k=-k j=i, ki=-ik=j, ii=jj=kk=ijk=-1.

Dynamics of (2 + 1) dimensional space can be expressed by quaternions, but for (3 + 1) D dynamics, we need biquaternions or split - quaternions .

3.1 Biquaternion Function

Dirac (1945) studied Lorentz transformation of bi-quaternion and expressed quaternion as the ratio of u to v, and assigned Q1=u vbar, where vbar is the conjugate of v, Q2=u ubar, and Q3=v vbar.
Mathematical basis of quaternions was studied by Connes’ Noncommutative Geometry (1994) and Woronowicz’s Non-Commutative Differential Calculus (1986).

The ProductLog function has a branch cut on the complex plane from -1/e to -infinity.
When u is replaced by v, Q1 changes to Q3’, Q3 changes to Q1’, Q2 changes to Q1bar=v ubar, which can be assigned as a branch of Q2’. 
When vbar is replaced by ubar, Q1 changes to Q2’’, Q2 changes to Q1’’, Q3 changes to v ubar, which can be assigned as a branch of Q3’’.
The replacement does not change elements of quaternions. Only the branch lines change to -ω/e to infinity and -ω^2/e to infinity.

ComplexPlot3D[ProductLog[z],{z,-2-2 I,2+2 I},PlotLegends->Automatic]

(Legended[, (BarLegend)])

It is a new kind of  triality symmetry.

3.2 Split-Quaternion Function

Godberashvili (2014) defined split-Quaternion
q=q0+q1 i+(q2+q3 i)j, q*=q0-q1 i-(q2+q3 i)j.
We consider also
q=q0+q2 j+(q3+q1 j)k, q*=q0-q2 j-(q3+q1 j)k,
q=q0+q3 k+(q1+q2 k)i, q*=q0-q3 k-(q1+q2 k)i.

The fact that ij = k can be taken as the time axis induces the conformal  O (2, 1) symmetry.
Time Reversal based Nonlinear Elastic Wave Spectroscopy (TR-NEWS) experiments confirm the conformality of (2+1)D spectra.

3.3 Spacelike time based quaternion

Tomonaga (1946) proposed Lorentz covariant dynamics using space - like time. 
In the propagation of ultrasonic waves in (3+1)D, we consider paths of the Fixed Point action of lattice QCD adopted by DeGrand et al. Restricting the length of path be less than or equal to 8 lattice units in one direction, there are 7 paths which contain time delay defined by the space-like time. We place Weyl fermions on the lattice cites, and find the weight function of the 7 paths which yield minimum action using Machine Learning techniques.

Since 1+ω+ω^2=0, chaos plays an important role in spectra of (3+1)D. The action depends on chirality of paths. By comparison of experimental spectra and simulation results, one may be able to verify chiral symmetry violation in ultrasonic wave propagation.

4. Summary and Outlook

1. Langevin equation of Ornstein - Uhlenbeck process is a fractional differential equation.    The Fokker - Planck equation for this process gives a probability distribution  following the Ito' s stochastic differential equation . Hysteresis effects can be formulated from this equation .
  
  2. The Caputo Fractional differential has left-differential and right-differential, which are not same in general. Mathematica calculates the terms
 properly.
     
  3. The Gegenbauer polynomial has a branch cut in a complex plane. The split - quaternion and Dirac' s quaternion representation by u vbar suggest that there are three quaternions Q1’,Q2’ and Q3’ on extended complex planes which have  new type of triality . 
  In (2 + 1) D ultrasonic wave spectra show conformality, which can be reproduced by quaternions.
 In (3 + 1) D ultrasonic wave propagation, lattice method using split - quaternions is promising.  For medical applications, fast optimization techniques are necessary .

Thank you for your attention .

References
[1] H.Nakajima and S.Furui, ‘A new algorithm for numerical simulation of Langevin equations’, Nucl. Phys. B(Proc. Suppl.)53, 983-986 (1997).
[2] K.Ito, ‘The Probability Theory’(In Japanese), Iwanami Shoten (1999).
[3] A.S. Vatsala et al. ‘Analysis of Sequential Caputo Fractional Differential Equations versus Non-Sequential Caputo Fractional Differential Equations with Applications,  Foundations, 2, 1129-1142 (2022).
[4] H.Kiskinov et al., ‘Comment on Valsera et al. Analysis of Sequential Caputo Fractional Differential Equations versus Non-Sequential Caputo Fractional Differential Equations with Applications. Foundations 2022,2,1129-1142’, Foundations 2024,4, 488-490.
[5] K.T.Elgindy and H.M. Refat, ‘High-order shifted Gegenbauer integral pseudo-spectral method for solving differential equations of Lane-Emden type’, Applied Numerical Mathematics, 128, 98-124 (2018).
[6] K.T.Elgindy, ‘Super-Exponential Approximation of the Riemann-Liouville Fractional Integral via Gegenbauer-Based FRactional Approximation Methods’, Algorithms , 18 (7) 395 (2025).
[7] P.A.M.Dirac, ‘Application of Quaternions to Lorentz Transformations’, Proc.R.I.A. vol. L. sect. A, 261-270 (1945).
[8] M.Godberashvili, ‘Split quaternions and particles in (2+1)D-space’, Eur.Phys.J. C 74, 3200 (2014).
[9] S.Tomonaga, ‘On a Relativistically Invariant Formulation of the Quantum Theory of Wave Fields, Prog. Theor. Phys. 1 (2) 27-40 (1946).
[11] S.Furui and S. Dos Santos, ‘Quaternion Differential Equation for the (2+1)D and (3+1)D NDT, ICU2025 congress proceedings (2025).
[12] A.Connes, ‘Noncommutative Geometry’ , Academic Press, San Diego (1990).
[13] S.I. Woronowicz, ‘Twisted SU(2) Group, An Example of a Non-Commutative Differential Calculus’, Publ. RIMS, Kyoto Univ. 23, 117-181 (1987).