Fourier Transform on the Homogeneous Space of 3D Positions and Orientations for Exact Solutions to Linear PDEs

Fokker–Planck PDEs (including diffusions) for stable Lévy processes (including Wiener processes) on the joint space of positions and orientations play a major role in mechanics, robotics, image analysis, directional statistics and probability theory. Exact analytic designs and solutions are known in the 2D case, where they have been obtained using Fourier transform on SE(2). Here, we extend these approaches to 3D using Fourier transform on the Lie group SE(3) of rigid body motions. More precisely, we define the homogeneous space of 3D positions and orientations R3⋊S2:=SE(3)/({0}×SO(2)) as the quotient in SE(3). In our construction, two group elements are equivalent if they are equal up to a rotation around the reference axis. On this quotient, we design a specific Fourier transform. We apply this Fourier transform to derive new exact solutions to Fokker–Planck PDEs of α-stable Lévy processes on R3⋊S2. This reduces classical analysis computations and provides an explicit algebraic spectral decomposition of the solutions. We compare the exact probability kernel for α=1 (the diffusion kernel) to the kernel for α=12 (the Poisson kernel). We set up stochastic differential equations (SDEs) for the Lévy processes on the quotient and derive corresponding Monte-Carlo methods. We verified that the exact probability kernels arise as the limit of the Monte-Carlo approximations.


Introduction
The Fourier transform has had a tremendous impact on various fields of mathematics including analysis, algebra and probability theory. It has a broad range of applied fields such as signal and image processing, quantum mechanics, classical mechanics, robotics and system theory. Thanks to Jean-Baptiste Joseph Fourier (1768-1830), who published his pioneering work "Théory analytique de la chaleur" in 1822, the effective technique of using a Fourier transform to solve linear PDE-systems (with appropriate boundary conditions) for heat transfer evolutions on compact subsets Ω of R d was born. The Fourier series representations of the solutions helped to understand the physics of heat transfer. Due to the linearity of the evolution operator that maps the possibly discontinuous square integrable initial condition to the square integrable solution at a fixed time t > 0, one can apply a spectral decomposition which shows how each eigenfunction is dampened over time. Thanks to contributions of Johann Peter Gustav Lejeune Dirichlet (1805-1859), completeness of the Fourier basis could then be formalized for several boundary conditions. Indeed, separation of variables (also known as "the Fourier method") directly provides a Sturm-Liouville problem [1] and an orthonormal basis of eigenfunctions for L 2 (Ω), which is complete due to compactness of the associated self-adjoint kernel intriguing, from the perspective of geometric theory of information and heat [76], to study optimal entropy on R 3 S 2 and (Fourier) Cramér Transforms building on results [77,78] on R n . However, such investigations first require a good grip on the spectral decompositions of the PDE-evolution operators for α-stable Lévy processes via a Fourier transform on R 3 S 2 , which is our primary focus here.

Structure of the Article
The structure of the article is as follows. In the first part of the Introduction, we briefly discuss the history of the Fourier transform, and its generalization to other groups that are the semi-direct product of the translation group and another matrix group, where we provide an overview of related works. Then, we specify our domain of interest-the Fourier transform on the homogeneous space G/H of positions and orientations, which is a Lie group quotient of the Lie group G = SE(3) with a subgroup H isomorphic to {0} × SO (2). Then, we address its application of solving PDE systems on G/H, motivated from applications in image analysis, robotics and probability theory.
There are four remaining subsections of the Introduction. In Section 1.2, we provide basic facts on the homogeneous space G/H of positions and orientations and we provide preliminaries for introducing a Fourier transform on G/H. In Section 1.3, we formulate the PDEs of interest on G/H that we solve. In Section 1.4, we formulate the corresponding PDEs on the group G. In Section 1.5, we relate the PDE for α = 1 2 to a Poisson system and quantify monotonic increase of entropy for α ∈ { 1 2 , 1}. In Section 1.6, we provide a roadmap on the spectral decomposition of the PDE evolutions.
In Section 2, based on previous works, we collect the necessary prior information about the PDEs of interest and the corresponding kernels. We also describe how to extend the case α = 1 (the diffusion case) to the general case α ∈ (0, 1]. In Section 3, we describe the Fourier transform on the Lie group SE (3), where we rely on UIRs of SE (3). In particular, by relating the UIRs to the dual orbits of SO(3) and by using a decomposition with respect to an orthonormal basis of modified spherical harmonics, we recall an explicit formula for the inverse Fourier transform.
In Section 4, we present a Fourier transform F G/H on the quotient G/H = R 3 S 2 . Our construction requires an additional constraint-an input function must be bi-invariant with respect to subgroup H, as explained in Remark 3. This extra symmetry constraint is satisfied by the PDE kernels of interest. We prove a theorem, where we present: (1) a matrix representation for the Fourier transform on the quotient; (2) an explicit inversion formula; and (3) a Plancherel formula.
In Section 5, we apply our Fourier transform on the quotient to solve the PDEs of interest. The solution is given by convolution of the initial condition with the specific kernels (which are the probability kernels of α-stable Lévy process). We find the exact formulas for the kernels in the frequency domain relying on a spectral decomposition of the evolution operator (involved in the PDEs). We show that this result can be obtained either via conjugation of the evolution operator with our Fourier transform on R 3 S 2 or (less efficiently) via conjugation of the evolution operator with the Fourier transform acting only on the spatial part R 3 . Then, we present a numerical scheme to approximate the kernels via Monte-Carlo simulation and we provide a comparison of the exact solutions and their approximations. Finally, in Section 6, we summarize our results and discuss their applications. In the appendices, we address the probability theory and stochastic differential equations (SDEs) regarding Lévy processes on R 3 S 2 .
The main contributions of this article are: • We construct F R 3 S 2 -the Fourier transform on the quotient R 3 S 2 , in Equation (43).

•
The matrix representations for F R 3 S 2 , explicit inversion and Plancherel formulas are shown in Theorem 1.

•
The explicit spectral decompositions of PDE evolutions for α-stable Lévy process on R 3 S 2 , in the Fourier domains of both R 3 S 2 and R 3 , are shown in Theorem 2; here, the new spectral decomposition in the Fourier domain of R 3 S 2 is simpler and involves ordinary spherical harmonics.
• The quantification of monotonic increase of entropy of PDE solutions for α-stable Lévy processes on R 3 S 2 for α ∈ { 1 2 , 1} in terms of Fisher information matrices is shown in Proposition 1. • the exact formulas for the probability kernels of α-stable Lévy processes on R 3 S 2 , in Theorem 3. This also includes new formulas for the heat kernels (the case α = 1), that are more efficient than the heat kernels presented in previous work [40].

•
Simple formulation and verifications (Monte-Carlo simulations) of discrete random walks for α-stable Lévy processes on R 3 S 2 in Proposition 3. The corresponding SDEs are in Appendix A.

Introduction to the Fourier Transform on the Homogeneous Space of Positions and Orientations
Let G = SE(3) denote the Lie group of rigid body motions, equipped with group product: Here, σ −1 g = σ * g for all g ∈ G; and 3.
there does not exist a closed subspace V of H σ other than {0, H σ } such that σ g V ⊂ V.
We denote byĜ the dual group of G. Its elements are equivalence classes of UIRs, where one identifies elements via σ 1 ∼ σ 2 ⇔ there exists a unitary linear operator υ, s.t.
is a unimodular Lie group of type I, which means that the left and right-invariant Haar measure coincide, and that its dual group and its quasi dual group coincide. Thereby it admits a Plancherel theorem [22,24]. Definition 3. The Fourier transform F G ( f ) = ((F G f )(σ)) σ∈Ĝ of a square-integrable, measurable and bounded function f on G is a measurable field of bounded operators indexed by unitary irreducible representations (UIR's) σ. Now,Ĝ can be equipped with a canonical Plancherel measure ν and the Fourier transform F G admits an extension unitary operator from L 2 (G) to the direct-integral space For details, see [22,24], and, for detailed explicit computations, see [4].
In this article, we constrain and modify the Fourier transform F G on G = SE(3) such that we obtain a suitable Fourier transform F G/H defined on a homogeneous space of left cosets, where Stab SO(3) (a) = {R ∈ SO(3) | Ra = a} denotes the subgroup of SO(3) that stabilizes an a priori reference axis a ∈ S 2 , say a = e z = (0, 0, 1) T . In the remainder of this article, we set this choice a = e z .

Remark 1.
Although the semi-direct product notation R 3 S 2 is formally not correct as S 2 is not a Lie group, it is convenient: it reminds that G/H denotes the homogeneous space of positions and orientations. (3) denote equivalence classes of rigid body motions g = (x, R n ) ∈ SE(3) that map (0, a) to (x, n):

Remark 2. (notation and terminology) Elements in Equation
under the (transitive) action Therefore, we simply denote the equivalence classes [g] by (x, n). This is similar to the conventional writing n ∈ S 2 = SO(3)/SO (2). Throughout this manuscript, we refer to G/H as "the homogeneous space of positions and orientations" and henceforth R n denotes any rotation that maps the reference axis a into n.
The precise definition of the Fourier transform F G/H on the homogeneous space G/H is presented in Section 4. It relies on the decomposition into unitary irreducible representations in Equation (2), but we must take both a domain and a range restriction into account. This is explained in Section 4. Next, we address an a priori domain constraint that is rather convenient than necessary.

Remark 3.
We constrain the Fourier transform F G/H to This constraint is convenient in view of the PDEs of interest (and the symmetries of their kernels) that we formulate in the next subsection, and that solve via Fourier's method in Section 5.

Introduction to the PDEs of Interest on the Quotient
Our main objective is to use the Fourier transform F G/H to solve the following PDEs on R 3 S 2 : where (x, n) ∈ R 3 S 2 , t ≥ 0, α ∈ (0, 1] and the generator is expressed via with D 33 > D 11 ≥ 0, D 44 > 0, and with ∆ S 2 n the Laplace-Beltrami operator on S 2 = n ∈ R 3 n = 1 .
Note that the generator Q is a self-adjoint unbounded operator with domain where H 2 denotes the Sobolev space W 2 2 . The semigroup for α = 1 is a strongly continuous semigroup on L 2 (R 3 S 2 ) with a closed generator, and by taking the fractional power of the generator one obtains another strongly continuous semigroup, as defined and explained in a more general setting in the work by Yosida ([66], ch:11). The fractional power is formally defined by In Section 1.6, we show that the common technical representation Equation (8) is not really needed for our setting. In fact, it is very easy to account for α ∈ (0, 1] in the solutions; by a spectral decomposition, we only need to take fractional powers of certain eigenvalues in the Fourier domain. For the moment, the reader may focus on the case α = 1, where the system in Equation (6) becomes an ordinary elliptic diffusion system which is hypo-elliptic (in the sense of Hörmander [79]) even in the degenerate case where D 11 = 0.
The PDEs in Equation (6) have our interest as they are Forward-Kolmogorov equations for α-stable Lévy processes on G/H. See Appendix A for a precise formulation of discrete and continuous stochastic processes. This generalizes previous works on such basic processes [61,64] with applications in financial mathematics [80] and computer vision [65,78,81,82], from Lie group R 3 to the Lie group quotient R 3 S 2 . See Figure 1 for a visualization of sample paths from the discrete stochastic processes explained in Appendix A. They represent "drunk man's flights" rather than "drunk man's walks". (Top) random walks (or rather "drunk man's drives") and an iso-contour of the limiting diffusion kernel, for the case d = 2 studied in previous works (see, e.g., [15,37,83]); and (Bottom) random walks (or rather "drunk man's flights") and a visualization of the limiting distribution for the case d = 3. This limiting distribution is a degenerate diffusion kernel (x, n) → K α=1 t (x, n) that we study in this article. We visualize kernel K α=1 t by a spatial grid of surfaces, where all surfaces are scaled by the same µ > 0.

Reformulation of the PDE on the Lie Group SE(3)
Now, we reformulate and extend our PDEs in Equation (6) to the Lie group G = SE(3) of rigid body motions, equipped with group product in Equation (1). This helps us to better recognize symmetries, as we show in Section 2.1. To this end, the PDEs are best expressed in a basis of left-invariant vector fields {g → A i | g } 6 i=1 on G. Such left-invariant vector fields are obtained by push forward from the left-multiplication L g 1 g 2 := g 1 g 2 as where A i := A i | e form an orthonormal basis for the Lie algebra T e (G). We choose such a basis typically such that the first three are spatial generators A 1 = ∂ x , A 2 = ∂ y , A 3 = ∂ z = a · ∇ R 3 and the remaining three are rotation generators, in such a way that A 6 is the generator of a counter-clockwise rotation around the reference axis a. For allŨ ∈ C 1 (G) and g ∈ G, one has where A → e A denotes the exponent that maps Lie algebra element A ∈ T e (G) to the corresponding Lie group element. The explicit formulas for the left-invariant vector fields in Euler-angles (requiring two charts) can be found in Appendix B, or in [4,84]. Now we can re-express the PDEs in Equation (6) on the group G = SE(3) as follows: where the generatorQ is again a fractional power (α ∈ (0, 1]) of the diffusion generatorQ given bỹ . . , 5}. The initial condition in Equation (10) is given bỹ Similar to the previous works [40,85], one has that holds for all t ≥ 0, (x, R) ∈ SE(3).
so that the generator of the PDE in Equation (10) on G and the generator of the PDE in Equation (6) on G/H indeed stay related viaQ 1.5. Increase of Entropy for the Diffusion System (α = 1) and the Poisson System (α = 1 2 ) on G/H The PDE-system in Equation (6) on G/H relates to the PDE-system in Equation (10) on G via Equation (14). Next, we show that for α = 1 2 the PDE-system boils down to a Poisson system. For α = 1 the PDE-system in Equation (10) is a diffusion system on Lie group G, for which one has monotonic increase of entropy [75]. The next theorem quantifies the monotonic increase of entropy for α ∈ { 1 2 , 1} in terms of Fisher matrices.

Definition 4.
Let α ∈ (0, 1]. LetW α be the solution to Equation (10) with positive initial conditionŨ > 0 withŨ ∈ L 2 (G) and GŨ (g)dg = 1. Then, we define the entropy E α (t) at evolution time t ≥ 0 as Proposition 1. For α = 1 2 , the PDE system in (10) yields the same solutions as the following Poisson system: The entropy in Equation (15) For all t > 0, one has for the diffusion matrix D = diag{D ii } 6 i=1 > 0, where D 11 = D 22 , D 33 and D 44 = D 55 are the coefficients iñ Q, and with Fisher matrix Proof. For α = 1 2 , one has by the square integrability constraint in Equation (16) and application of the unitary Fourier transform on G that ∂ 2 = 0 and thereby the PDE system in Equation (10) can be replaced by the Poisson system in Equation (16) on G × R + . The formula for the entropy follows from a product decomposition of the (bi-invariant) haar measure on G into measure on the quotient G/H and a measure on the subgroup H ≡ {0} × SO(2) and the fact thatW α (gh, t) =W α (g, t) for all h ∈ H, α ∈ (0, 1], due to Equation (14). For α = 1 2 , we have thatW α satisfies Equation (16) and where we use integration by parts and short notation with ∂ t := ∂ ∂t . Now, E 1 2 < 0 and E 1 2 is continuous (due to the Lebesgue dominated convergence principle and continuity of each mapping t → ∂ tW (g, t) indexed by g ∈ G) and E 1 For α = 1, we follow ( [75], Thm.2) and compute (again using the PDE and integration by parts) Regarding the strict positivity in Equation (17), we note thatŨ > 0 ⇒W α > 0 and if E α (t) = 0 then this would imply thatW α (·, t) is constant, which violatesW α (·, t) ∈ L 2 (G) as G is not compact.

A Preview on the Spectral Decomposition of the PDE Evolution Operator and the Inclusion of α
Let U be in the domain of the generator Q α given by Equation (7), of our evolution Equation (6). For a formal definition of this domain, we refer to ( [86], Equation 9). Let its spatial Fourier transform be given by To the operator Q α , we associate the corresponding operator −(−B) α in the spatial Fourier domain by Then, direct computations show us that where, for each fixed ω ∈ R 3 , the operator − In this article, we employ Fourier transform techniques to derive a complete orthonormal basis Then, clearly, this basis is also an ONB of eigenfunctions for −(−B ω ) α , as we only need to take the fractional power of the eigenvalues. Indeed, once the eigenfunctions in Equation (22) and the eigenvalues are known, the exact solution of Equation (6) is given by (shift-twist) convolution with a probability kernel on R 3 S 2 . More precisely, the solutions of Equation (6) can be expressed as follows: with the probability kernel given by Here, the inner product in L 2 (S 2 ) is given by where µ S 2 is the usual Lebesgue measure on the sphere S 2 .
Remark 5. The eigenvalues λ l,m r only depend on r = ω due to the symmetry Φ l,m Rω (Rn) = Φ l,m ω (n) that one directly recognizes from Equations (21) and (23). Remark 6. The kernels K α t are the probability density kernels of stable Lévy processes on R 3 S 2 , see Appendix A.1. Therefore, akin to the R n -case [61,65], we refer to them as the α-stable Lévy kernels on R 3 S 2 .

Symmetries of the PDEs of Interest
Next, we employ the PDE formulation in Equation (10) on the group G = SE(3) to summarize the symmetries for the probability kernels K α t : R 3 S 2 → R + . For details, see [40,87].

PDE Symmetries
Consider the PDE system in Equation (10) on the group G = SE (3). Due to left-invariance (or rather left-covariance) of the PDE, linearity of the mapŨ(·) →W α (·, t), and the Dunford-Pettis theorem [88], the solutions are obtained by group convolution with a kernelK α t ∈ L 1 (G): where we take the convention that the probability kernel acts from the left. In the special case, U = δ e with unity element e = (0, I) we getW α (g, t) =K α t (g). Thanks to the fundamental relation in Equation (13) that holds in general, we have in particular that ∀ t≥0 ∀ (x,R)∈G :K α To avoid confusion between the Euler angle α and the α indexing the α-stable Lévy distribution, we put an overline for this specific angle. Henceforth, R v,ψ denotes a counter-clockwise rotation over axis v with angle ψ. This applies in particular to the case where the axis is the reference axis v = a = (0, 0, 1) T and ψ = α. Recall that R n (without an angle in the subscript) denotes any 3D rotation that maps reference axis a onto n.
We write the symbol· above a function to indicate its Fourier transform on G and G/H; we use the symbol · for strictly spatial Fourier transform; the symbol· above a function/operator to indicate that it is defined on the group G and the function/operator without symbols when it is defined on the quotient G/H.

2.2.
Obtaining the Kernels with D 11 > 0 from the Kernels with D 11 = 0 In ( [40], cor.2.5), it was deduced that for α = 1 the elliptic diffusion kernel (D 11 > 0) directly follows from the degenerate diffusion kernel (D 11 = 0) in the spatial Fourier domain via For the general case α ∈ (0, 1], the transformation from the case D 11 = 0 to the case D 11 > 0 is in Equation (24) for the kernel. Henceforth, we set D 11 = 0.

The Fourier Transform on SE(3)
The group G = SE(3) is a unimodular Lie group (of type I) with (left-and right-invariant) Haar measure dg = dxdµ SO(3) (R) being the product of the Lebesgue measure on R 3 and the Haar measure Then, for all f ∈ L 1 (G) ∩ L 2 (G), the Fourier transform F G f is given by Equation (2). For more detailsm see [22,24,26]. One has the inversion formula: In our Lie group case of SE (3), we identify all unitary irreducible representations σ p,s having non-zero dual measure with the pair (p, s) ∈ R + × Z. This identification is commonly applied (see, e.g., [4]). Using the method ( [26], Thm. 2.1, [25]) of induced representations, all unitary irreducible representations (UIRs) of G, up to equivalence, with non-zero Plancherel measure are given by: where pS 2 denotes a 2D sphere of radius p = u ; ∆ s is a unitary irreducible representation of SO(2) (or rather of the stabilizing subgroup Stab SO(3) (a) ⊂ SO(3) isomorphic to SO (2)) producing a scalar. In Equation (32), R u p denotes a rotation that maps a onto u p . Thus, direct computation shows us that it is a rotation around the z-axis (recall a = e z ), e.g. about angle α. This yields character Mackey's theory [25] relates the UIR σ p,s to the dual orbits pS 2 of SO(3). Thereby, the dual measure ν can be identified with a measure on the family of dual orbits of SO(3) given by {pS 2 | p > 0}, and for all p > 0, s ∈ Z. For details, see ( [24], ch. 3.6.). The matrix elements off = F G f with respect to an orthonormal basis of modified spherical harmonics {Y l,m s (p −1 ·)}, with |m|, |s| ≤ l (see ([4], ch.9.8)) for L 2 (pS 2 ) are given bŷ where the L 2 inner product is given by (25)).
For an explicit formula for the modified spherical harmonics Y l,m s see [4], where they are denoted by h l m,s . The precise technical analytic expansion of the modified spherical harmonics is not important for this article. The only properties of Y l,m s that we need are gathered in the next proposition. (1) For s = 0 or m = 0, they coincide with standard spherical harmonics Y l,m , cf. ( [89], eq.4.32): with P m l the normalized associated Legendre polynomial and m = (−1) 1 2 (m+|m|) .
with P l m m a generalized associated Legendre polynomial given in ( [4], eq.9.21). Moreover, we have inversion formula ( [4], Equation 10.46): with matrix coefficients (independent of f ) given by Note that σ p,s is a UIR so we have

A Specific Fourier Transform on the Homogeneous Space R 3 S 2 of Positions and Orientations
Now that we have introduced the notation of Fourier transform on the Lie group G = SE (3), we define the Fourier transform F G/H on the homogeneous space G/H = R 3 S 2 . Afterwards, in the subsequent section, we solve the Forward-Kolmogorov/Fokker-Planck PDEs in Equation (6)  Throughout this manuscript, we rely on a Fourier transform on the homogeneous space of positions and orientations that is defined by the partition of left-cosets: R 3 S 2 := G/H, given by Equation (3).
Note that subgroup H can be parameterized as follows: where we recall that R a,α denotes a (counter-clockwise) rotation around the reference axis a = e z . The reason behind this construction is that the group SE(3) acts transitively on R 3 S 2 by (x , n ) → g (x , n ) given by Equation (4). Recall that by the definition of the left-cosets one has The latter equivalence simply means that for g 1 = (x 1 , R 1 ) and g 2 = (x 2 , R 2 ) one has The equivalence classes [g] = {g ∈ SE(3) | g ∼ g} are often just denoted by (x, n) as they consist of all rigid body motions g = (x, R n ) that map reference point (0, a) onto (x, n) ∈ R 3 S 2 : where we recall R n is any rotation that maps a ∈ S 2 onto n ∈ S 2 .

Fourier Transform on R 3 S 2
Now we can define the Fourier transform F G/H on the homogeneous space G/H. Prior to this, we specify a class of functions where this transform acts.

Definition 5.
Let p > 0 be fixed and s ∈ Z. We denote the subspace of spherical functions that have the prescribed axial symmetry, with respect to the subgroup H (recall Equation (38)).

Definition 6.
We denote the orthogonal projection from L 2 (pS 2 ) onto the closed subspace L sym 2 (pS 2 ) by P sym p .

Definition 7.
To the group representation σ p,s : SE(3) → B(L 2 (pS 2 )) given by Equation (32), we relate a "representation" σ p,s : R 3 S 2 → B(L 2 (pS 2 )) on R 3 S 2 , defined by σ p,s To each function U : G/H → C, we relate an axially symmetric functionŨ : Definition 9. We define the Fourier transform of function U on G/H = R 3 S 2 bŷ Standard properties of the Fourier transform F G on SE(3) such as the Plancherel theorem and the inversion formula [4,26] naturally carry over to F G/H with "simpler formulas". This is done by a domain and range restriction via the projection operators P sym p in Equation (43). The reason for the specific construction Equation (43)  Conversely, ifŨ = F −1 G (Û) and thenŨ satisfies the axial symmetry in Equation (41).
Proof. Item 1: Uniqueness of U follows by the fact that the choice of R n of some rotation that maps a onto n does not matter. Indeed, U(x, n) =Ũ(x, R n R a,α ) =Ũ(x, R n ). (41) can be rewritten asŨ(g) =Ũ(gh α ) for all h α ∈ H, g ∈ G. This gives:Û p,s l,m,l ,m

Item 2: Assumption Equation
where we recall that σ is a UIR and that the Haar measure on G is bi-invariant. In the first step, we used the third property, whereas in the final step we used the second property of Proposition 2 together with We conclude that (1 − e −imα )Û Thus, it remains to show why Equation (48) and thereby Equation (48) follows by Equation (36).

Lemma 2.
IfK ∈ L 2 (G) is real-valued and satisfies the axial symmetry in Equation (41), and moreover the following holdsK (g −1 ) =K(g) Proof. The proof follows by Equation (37) and inversion invariance of the Haar measure on G (see [86]).
The next lemma shows that Equation (50) is a sufficient but not a necessary condition for the Fourier coefficients to vanish for both the cases m = 0 and m = 0. Lemma 3. LetK ∈ L 2 (G) and K ∈ L 2 (G/H) be related by Equation (42). Then, we have the following equivalences: ,α x, R a,α n), for all α ∈ [0, 2π), (x, n) ∈ G/H K (gh) =K(g) =K(hg), for all g ∈ G, h ∈ H The Fourier coefficientsK p,s l,m,l ,m vanish for m = 0 and for m = 0.
= K(R a,α x, R a,α Ra) =K(R a,α x, R a,α R) =K(h α g). b ⇒ c: By Lemma 1, we know that the Fourier coefficients vanish for m = 0. Next, we show they also vanish for m = 0. Similar to Equation (49) we have which gives the following relation for the matrix-coefficients: The implication can be directly verified by Proposition 2, Equations (34) and (52), and From Equation (53), we deduce that: c ⇒ a: By inversion of Equation (35), where the only contributing terms have m = 0 and m = 0, we see thatK(gh) =K(hg) =K(g) for all h = (0, R a,α ). Thereby,K is axially symmetric and by Lemma 1 it relates to a unique kernel on G/H via K(x, n) =K(x, R n ) and the result follows by Equation (30). Now that we have characterized all functions K ∈ L 2 (G/H) for which the Fourier coefficientŝ K p,s l,m,l ,m vanish for m = 0 and m = 0 in the above lemma, we considerably simplify the inversion and Plancherel formula for Fourier transform F G on the group G = SE(3) to the Fourier transform F G/H on the homogeneous space G/H = R 3 S 2 in the next theorem. This is important to our objective of deriving the kernels for the linear PDEs in Equation (6) that we address in the next section.  (34). Furthermore, we have the following Plancherel and inversion formula: l,0,l ,0 | 2 p 2 dp, for g = (x, R n ).
Proof. The above formulas are a direct consequence of Lemma 3 and the Plancherel and inversion formulas (see [4], ch:10.8, [26]) for Fourier transform on SE(3). Recall that a coordinate-free definition of σ p,s is given in Equation (40). Its matrix coefficients are given by Equation (54), where we recall the first item of Proposition 2 and where we note that they are independent on the choice of R n ∈ SO(3) mapping a onto n. Corollary 1. Let K 1 , K 2 ∈ L sym 2 (G/H). Then, for shift-twist convolution on G/H = R 3 S 2 given by Proof. SetK 1 (g) = K 1 (g (0, a)). Standard Fourier theory [5] gives where the first equality is given by Equation (43) and the third equality follows by Lemma 3 and Equation (47).

Application of the Fourier Transform on R 3 S 2 for Explicit Solutions of the Fokker-Planck PDEs of α-stable Lévy Processes on R 3 S 2
Our objective is to solve the PDE system in Equation (6) on the homogeneous space of positions and orientations G/H. Recall that we extended this PDE system to G in Equation (10). As the cases D 11 > 0 follow from the case D 11 = 0 (recall Section 2.2), we consider the case D 11 = 0 in this section. From the symmetry consideration in Section 2, it follows that the solution of Equation (10) is given bỹ W α (g, t) = (K α t * Ũ)(g) with a probability kernelK α t : G → R + , whereas the solution of Equation (6) is given by where the kernels K α t are invariant with respect to left-actions of the subgroup H (recall Equation (30)). This invariance means that the condition for application of the Fourier transform F G/H on R 3 S 2 is satisfied (recall Lemma 3) and we can indeed employ Theorem 1 to keep all our computations, spectral decompositions and Fourier transforms in the 5D homogeneous space R 3 S 2 = G/H rather than a technical and less direct approach [40] in the 6D group G = SE(3).

Remark 8.
For the underlying probability theory, and sample paths of discrete random walks of the α-Stable Lévy stochastic processes, we refer to Appendix A. To get a general impression of how Monte Carlo simulations of such stochastic processes can be used to approximate the exact probability kernels K α t , see Figure 1. In essence, such a stochastic approximation is computed by binning the endpoints of the random walks. A brief mathematical explanation follows in Section 5.2.
For now, let us ignore the probability theory details and let us first focus on deriving exact analytic solutions to Equation (6) and its kernel K α t via Fourier transform F G/H on G/H = R 3 S 2 .

Exact Kernel Representations by Spectral Decomposition in the Fourier Domain
Let us consider the evolution in Equation (6) for α-stable Lévy process on the quotient G/H = R 3 S 2 . Then, the mapping from the initial condition W(·, 0) = U(·) ∈ L 2 (G/H) to the solution W(·, t) at a fixed time t ≥ 0 is a bounded linear mapping. It gives rise to a strongly continuous (holomorphic) semigroup [66]. We conveniently denote the bounded linear operator on L 2 (G/H) as follows: W α (·, t) = (e tQ α U)(·), for all t ≥ 0.
In the next main theorem, we provide a spectral decomposition of the operator using both a direct sum and a direct integral decomposition. Note that definitions of direct integral decompositions (and the underlying measure theory) can be found in ( [24], ch:3.3 and 3.4).

Eigenfunctions and Preliminaries
To formulate the main theorem, we need some preliminaries and formalities. First, let us define Recalling Equation (19), we re-express the generator in the spatial Fourier domain: where β ω denotes the angle between n and r −1 ω (see Figure 2). This re-expression is the main reason for the following definitions. Instead of the modified spherical Harmonics Y l,m s in Proposition 2, which are commonly used as a standard basis to represent each operator in the Fourier transform on SE(3), we use our generalized spherical harmonics, depending on a spatial frequency vector, as this is in accordance with Equation (57).

Definition 10.
Let l ∈ N 0 . Let m ∈ Z such that |m| ≤ l. Let ω ∈ R 3 be a frequency vector. We define where we take the rotation which maps a onto r −1 ω whose matrix representation in the standard basis is:
Recall the standard spherical angle formula n(β, γ) = (sin β cos γ, sin β sin γ, cos β) T from Proposition 2. These are Euler-angles relative to the reference axis a = e z . For the Euler-angles relative to the (normalized) frequency r −1 ω one has (see also Figure 2): Definition 11. Let l ∈ N 0 . Let m ∈ Z such that |m| ≤ l. We define the functions Φ l,m where r = ω and d l,m (r) := d l,m j (r) ∞ j=0 are coefficients such that where S l,m ρ (·) denotes the L 2 -normalized spheroidal wave function.

Remark 9.
The spheroidal wave function arises from application of the method of separation on operator B ω in Equation (57) where basic computations (for details, see [40]) lead to the following singular Sturm-Liouville problem: 44 . In this formulation, p(x) vanishes at the boundary of the interval, which makes our problem a singular Sturm-Liouville problem. It is sufficient to require boundedness of the solution and its derivative at the boundary points to have nonnegative, distinct, simple eigenvalues λ l,m r and existence of a countable, complete orthonormal basis of eigenfunctions {y j } ∞ j=0 [91] for the spheroidal wave equation.

The Explicit Spectral Decomposition of the Evolution Operators
In Theorem 2, we present the explicit spectral decompositions both in the spatial Fourier domain and in the Fourier domain of the homogeneous space of positions and orientations.
Prior to this theorem, we explain the challenges that appear when we apply F G/H to the PDE of interest in Equation (6) on the quotient G/H. To get a grip on the evolution operator and the corresponding kernel, we set the initial condition equal to a delta distribution at the origin, i.e., we consider In this case, the necessary condition in Equation (51) in Lemma 3 for application of F G/H is indeed satisfied, due to the symmetry property of the kernel, given by Equation (30). Now, due to linearity we just need to study the generator in the Fourier domain.
For the moment, we set α = 1 (the degenerate diffusion case) and return to the general case later on (recall Sections 1.6 and 2.2). Then, it follows that (for details, see ( [40], App.D)) Here, ∆ pS 2 u denotes the Laplace-Beltrami operator on a sphere pS 2 = {u ∈ R 3 u = p} of radius p > 0.
We recall that u ∈ pS 2 is the variable of the functions on which σ p,s acts. Recalling Equation (32), the first part in the righthand side of Equation (65) denotes a multiplier operator M given by (Mφ)(u) := − (a · u) 2 φ(u), for all φ ∈ L 2 (pS 2 ), and almost every u ∈ pS 2 .
As a result, we obtain the following PDE system forK α t (now for general α ∈ (0, 1]):

Remark 12.
There is a striking analogy between the operators F G/H • Q α • F −1 G/H and F R 3 • Q α • F −1 R 3 given by Equation (57), where the role of rR T ω/r n corresponds to u. This correspondence ensures that the multipliers of the multiplier operators in the generator coincide and that the roles of p and r coincide: Proof. Recall Equation (63) that defines matrix M m (for analytic formulas of this tri-diagonal matrix, see [40]). This may be re-written as follows: Now, fix s ∈ Z and set m = s and n = p −1 u and we have: where again l = |s| + j, l = |s| + j and j, j ∈ N 0 .
Finally, we note that operator D 33 M + D 44 ∆ pS 2 u is negative definite and maps each subspace span {Y l,s (p −1 ·)} ∞ l=|s| for fixed s ∈ Z onto itself, which explains direct sum decomposition in Equation (66).
Next, we formulate the main result, where we apply a standard identification of tensors a ⊗ b with linear maps: Theorem 2. We have the following spectral decompositions for the Forward-Kolomogorov evolution operator of α-stable Lévy-processes onthe homogeneous space G/H = R 3 S 2 : • In the Fourier domain of the homogeneous space of positions and orientations, we have: • In the spatial Fourier domain, we have where W(ω, ·, t) = F R 3 W(ω, ·, t) and U(ω, ·) = F R 3 U(ω, ·) (recall Equation (56)).
In both cases, the normalized eigenfunctions Φ l,m ω are given by Equation (60) in Definition 11. The eigenvalues λ l,m r are the eigenvalues of the spheroidal wave equation, as explained in Remark 9.
Proof. The first identity in Equation (68) follows by: In the last equality, we use the fact that Φ l,m a = Y l,m . By applying the identification in Equation (67), one observes that Equation (69) is a reformulation of Equation (24), was already been derived for α = 1 in previous work by the first author with J.M. Portegies ([40], Thm.2.3 and Equation31). The key idea behind the derivation, the expansion and the completeness of the eigenfunctions {Φ l,m ω } is summarized in Remark 9. The general case α ∈ (0, 1] then directly follows by Section 1.6. Recently, exact formulas for the (degenerate) heat-kernels on G = SE(3) and on G/H = R 3 S 2 (i.e., the case α = 1) have been published in [40]. In the next theorem: • We extend these results to the kernels of PDE in Equation (6), which are Forward-Kolmogorov equations of α-stable Lévy process with α ∈ (0, 1].

•
We provide a structured alternative formula via the transform F G/H characterized in Theorem 1.
Equation (71) follows similarly by δ a and the result follows from the second part of Theorem 1 (again by taking U a sequence that is a bounded approximation of the unity centered around (0, a)).

Remark 13.
There also exist Gaussian estimates for the heat kernel K α=1 t that use a weighted modulus on the logarithm on G, [92]. Such Gaussian estimates can account for the quotient structure G/H [87], and can be reasonably close (cf. [44], Figure 4.4, [93]) to the exact solutions for practical parameter settings in applications [48,94,95].

Monte-Carlo Approximations of the Kernels
A stochastic approximation for the kernel K α t is computed by binning the endpoints of discrete random walks simulating α-stable processes on the quotient R 3 S 2 that we explain next. Let us first consider the case α = 1. For M ∈ N fixed, we have the discretization stochastically independent Gaussian distributed on R with t = 1; with uniformly distributed γ k ∼ Unif (R/(2πZ) ≡ [−π, π)); and β k ∼ g, where g : R → R + equals g(r) = |r| 2 e − r 2 4 in view of the theory of isotropic stochastic processes on Riemannian manifolds by Pinsky [96]. By the central limit theorem for independently distributed variables with finite variance it is only the variances of the distributions for the random variables g and G R t=1 that matter. One may also take k ∼ √ 3 Unif − 1 2 , 1 2 and β k ∼ These processes are implemented recursively; for technical details and background, see Appendix A. (72) can be re-expressed, up to order 1 M for M 0, as follows:

Proposition 3. The discretization of Equation
with i k ∼ G R t=1 stochastically independent normally distributed variables with t = 1 2 σ 2 = 1, and D 44 = D 55 .
Proof. In our construction, β k and γ k can be seen as the polar radius and the polar angle (on a periodic square [−π, π] × [−π, π]) of a Gaussian process with t = 1 on a plane spanned by rotational generators A 4 and A 5 . The key ingredient to obtain Equation (73) from Equation (72) is given by the following relation: e u cos vA 5 −u sin vA 4 = e vA 6 e uA 5 e −vA 6 , for all u, v ∈ R, For the binning, we divide R 3 into cubes c ijk , i, j, k ∈ Z, of size ∆s × ∆s × ∆s: We divide S 2 into bins B l , l = {1, . . . , b} for b ∈ N, with surface area σ B l and maximal surface area σ B . The number of random walks in a simulation with traveling time t that have their end point x M ∈ c ijk with their orientation n M ∈ B l is denoted with # ijkl t . Furthermore, we define the indicator function When the number of paths N → ∞, the number of steps in each path M → ∞ and the bin sizes tend to zero, the obtained distribution converges to the exact kernel: The convergence is illustrated in Figure 3.

Comparison of Monte-Carlo Approximations of the Kernels to the Exact Solutions
In this section, we compute the probability density kernels K α t via the analytic approach of Section 5.1.2 (Equation (71), Theorem 3) and via the Monte-Carlo approximation of Section 5.2. The kernels are computed on a regular grid with each (x i , y j , z k ) at the center of the cubes c ijk of Equation (75) with i, j = −3, . . . , 3, k = −5, . . . , 5, and ∆s = 0.5. The Monte-Carlo simulations also require spherical sampling which we did by a geodesic polyhedron that sub-divides each mesh triangle of an icosahedron into n 2 new triangles and projects the vertex points to the sphere. We set n = 4 to obtain 252 (almost) uniformly sampled points on S 2 .
The exact solution is computed using (truncated) spherical harmonics with l ≤ 12. To obtain the kernel, we first solve the solution in the spatial Fourier domain and then do an inverse spatial Fast Fourier Transform. The resulting kernel K α t (where we literally follow Equation (71)) is only spatially sampled and provides for each (x i , y j , z k ) an analytic spherical distribution expressed in spherical harmonics.
For the Monte-Carlo approximation, we follow the procedure described in Section 5.2. The kernel K α t is obtained by binning the end points of random paths on the quotient R 3 S 2 (cf. Equation (72)) and thereby approximate the limit in Equation (76). Each path is discretized with M = 40 steps and in total N = 10 10 random paths were generated. The sphere S 2 is divided into 252 bins with an average surface area of σ B l ≈ 4π 252 . In Figures 1 and 3-5, we set D 33 = 1, D 44 = 0.2. In the comparison between the kernels K α=1 t with K α=0.5 t , we set t = 2 and t = 3.5, respectively, to match the full width at half maximum value of the distributions. In Figures 1, 3 and 5, we set α = 1 and t = 2. In Figures 1, 3 and 4, we sample the grid in Equation (75) with |i|, |j| ≤ 4, |k| ≤ 8. Figure 5 shows that the Monte-Carlo kernel closely approximates the exact solution and since the exact solutions can be computed at arbitrary spherical resolution, it provides a reliable way to validate numerical methods for α-stable Lévy processes on R 3 S 2 .

Conclusions
We set up a Fourier transform F G/H on the homogeneous space of positions and orientations. The considered Fourier transform acts on functions that are bi-invariant with respect to the action of subgroup H. We provide explicit formulas (relative to a basis of modified spherical harmonics) for the transform, its inverse, and its Plancherel formula, in Theorem 1.
Then, we use this Fourier transform to derive new exact solutions to the probability kernels of α-stable Lévy processes on G/H, including the diffusion PDE for Wiener processes, which is the special case α = 1. They are obtained by spectral decomposition of the evolution operator in Theorem 2. New formulas for the probability kernels are presented in Theorem 3. There, the general case 0 < α < 1 follows from the case α = 1 by taking the fractional power of the eigenvalues. In comparison to previous formulas in [40] for the special case α = 1 obtained via a spatial Fourier transform, we have more concise formulas with a more structured evolution operator in the Fourier domain of G/H, where we rely on ordinary spherical harmonics, and where we reduce the dimension of the manifold over which it is integrated from 3 to 1 (as can be seen in Theorem 3).
We introduce stochastic differential equations (or rather stochastic integral equations) for the α-stable Lévy processes in Appendix A.1, and we provide simple discrete approximations where we rely on matrix exponentials in the Lie group SE(3) in Proposition 3. We verified the exact solutions and the stochastic process formulations, by Monte-Carlo simulations that confirmed to give the same kernels, as shown in Figure 5. We also observed the expected behavior that the probability kernels for 0 < α < 1 have heavier tails, as shown in Figure 4.
The PDEs and the probability kernels have a wide variety of applications in image analysis (crossing-preserving, contextual enhancement of diffusion-weighted MRI, cf. [45,46,49,94,97,98] or in crossing-preserving diffusions in 3D scalar images [18]), robotics [4,5,57] and probability theory [56,61]. The generalizations to α ∈ (0, 1] allow for longer range interactions between local orientations (due to the heavy tails). This is also of interest in machine learning, where convolutional neural networks on the homogeneous space of positions and orientations [9,12] can be extended to 3D [67,68], which may benefit from the PDE descriptors and the Fourier transform presented here.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations and symbols are used in this manuscript: