Spectral decomposition of discrepancy kernels on the Euclidean ball, the special orthogonal group, and the Grassmannian manifold
Josef Dick Martin Ehler Manuel Gräf Christian Krattenthaler




To numerically approximate Borel probability measures by finite atomic measures, we study the spectral decomposition of discrepancy kernels when restricted to compact subsets of ^d. For restrictions to the Euclidean ball in odd dimensions, to the rotation group (3), and to the Grassmannian manifold _2,4, we compute the kernels' Fourier coefficients and determine their asymptotics. The L_2-discrepancy is then expressed in the Fourier domain that enables efficient numerical minimization based on the nonequispaced fast Fourier transform. For (3), the nonequispaced fast Fourier transform is publicly available, and, for _2,4, the transform is derived here. We also provide numerical experiments for (3) and _2,4. 




[2010]Primary 42C10; Secondary
42B10, 43A65, 43A85, 65T40













Viscosity solutions for
controlled McKean–Vlasov jump-diffusions

    Matteo BurzoniMathematical Institute, University of Oxford. Partly supported by the Hooke Research Fellowship from the University of Oxford.1[4]
    Vincenzo IgnazioDepartment of Mathematics, ETH Zürich.1[4]
A. Max ReppenDepartment of Operations Research and Financial Engineering, Princeton University. Supported by the Swiss National Science Foundation grant SNF 181815.
H. Mete Soner[2]1Partly supported by the Swiss National Science Foundation grant SNF 200020_172815.1Partly supported by the ETH Foundation and the Swiss Finance Institute.

    December 30, 2023
=====================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================






§ INTRODUCTION

Consider a Borel probability measure μ:ℬ(^d)→[0,1] on ^d, where ℬ(^d) denotes the Borel sigma algebra on ^d. For fixed n∈, we aim to allocate a suitable n-point set {x_1,…,x_n}⊂^d such that the normalized atomic measure 



    ν_n:=1/n∑_j=1^nδ_x_j

approximates μ. Here, δ_x_j:ℬ(^d)→{0,1} denotes the point measure localized at x_j. To quantify the L_2-discrepancy between μ and ν_n, select a measure β on ℬ(^d) with μ,δ_x∈ L_2(ℬ(^d),β), for all x∈^d, and consider

    𝒟_β(μ,ν_n):=μ-ν_n^2_L_2(ℬ(^d),β)=∫_ℬ(^d)| μ(B)- ν_n(B)|^2 dβ(B),

cf. <cit.>, see Section <ref> for explicit examples. For fixed n∈, we aim to minimize 𝒟_β(μ,ν_n) among all n-point sets {x_1,…,x_n}⊂^d. The present manuscript is concerned with discretizations of (<ref>) that facilitate numerical minimization.

The associated discrepancy kernel K_β:^d×^d→ is defined by 

    K_β(x,y):= ⟨δ_x,δ_y⟩_L_2(ℬ(^d),β)=∫_ℬ(^d)δ_x(B)δ_y(B)  dβ(B),

and we assume it is continuous. Fubini's Theorem and μ(B)=∫_^dδ_x(B) dμ(x)
applied to 

    μ-ν_n^2_L_2(ℬ(^d),β)=μ^2_L_2(ℬ(^d),β)-2⟨μ,ν_n⟩_L_2(ℬ(^d),β)+ν_n^2_L_2(ℬ(^d),β)

yields that (<ref>) is identical to

    𝒟_β(μ,ν_n)=∬_^d×^d K_β(x,y) dμ(x)  dμ(y) -2∑_j=1^n∫_^dK_β(x,x_j)/n dμ(x)
    +∑_i,j=1^n K_β(x_i,x_j)/n^2 .


If a compact set ⊂^d is known in advance such that (μ)⊂, then we shall restrict the minimization to {x_1,…,x_n}⊂, so that only the restricted kernel K_β|_× matters. By endowing  with a finite Borel measure σ_ having full support, Mercer's Theorem yields an orthonormal basis {ϕ_l}_l=0^∞ for L_2(,σ_) and coefficients (a_l)_l=0^∞ such that the spectral decomposition 

    K_β|_×(x,y) = ∑_l=0^∞ a_l ϕ_l(x)ϕ_l(y),   x,y∈,

holds with absolute and uniform convergence. We call (a_l)_l=0^∞ the Fourier coefficients of the kernel K_β|_×. If (μ),(ν_n)⊂, then the Fourier expansion of the L_2-discrepancy (<ref>) is

    𝒟_β(μ,ν_n) = ∑_l=0^∞ a_l  | μ̂_l-ν̂_n,l|^2,
        μ̂_l:=∫_ϕ_l(x) dμ(x),  ν̂_n,l:=1/n∑_j=1^nϕ_l(x_j),

where the Fourier coefficients μ̂_l and ν̂_n,l of the measures μ and ν_n, respectively, are well-defined if a_l≠ 0.  
Truncation of the discretization (<ref>) enables the use of the nonequispaced fast Fourier transform, thereby offering more efficient minimization of 𝒟_β(μ,ν_n), cf. <cit.>. Thus, we aim to 

A) compute (a_l)_l=0^∞ and (ϕ_l)_l=0^∞ in the Fourier expansion (<ref>) of K_β|_×.







The L_2-discrepancy 𝒟_β(μ,ν_n) also coincides with the worst case integration error 

    𝒟_β(μ,ν_n) = sup_f_ℋ_β()≤ 1 |∫_ f(x) dμ(x) - 1/n∑_j=1^n f(x_j)  |^2

with respect to the reproducing kernel Hilbert space ℋ_β() generated by K_β|_×, 
cf. <cit.>. 

To specify ℋ_β(), we aim to 

B) identify ℋ_β() with a classical function space.

Fourier decay properties generally quantify Sobolev smoothness. To accomplish B), we aim to determine the asymptotics of K_β|_×'s Fourier coefficients (a_l)_l=0^∞ in (<ref>). 

For =^d-1 and a particular choice of β, the kernel K_β|_^d-1×^d-1 essentially coincides with the Euclidean distance, see <cit.>. The Fourier expansion is determined in <cit.>, and the 
decay of the Fourier coefficients yields that K_β|_^d-1×^d-1 reproduces the Sobolev space ℋ_β(^d-1)=ℍ^d/2(^d-1). For the sphere and the torus, the nonequispaced fast Fourier transform is available, and both A) and B) are discussed in <cit.>.  

This manuscript is dedicated to derive analogous results for other compact sets . We focus on the unit ball, the special orthogonal group, and the Grassmannian manifold, 

    𝔹^d   :={x∈^d:x≤ 1},
    (d)   :={x∈^d× d : (x)=1,  x^-1=x^⊤},
    𝒢_k,d   :={x∈^d× d : x^⊤=x,  x^2=x, (x)=k}.


We achieve goal A) for 𝕏=𝔹^d with odd d. Both goals,  A) and B), are achieved for (3) and _2,4. We also provide numerical experiments. For (3), the computations are based on the nonequispaced fast Fourier transform designed in <cit.>. For _2,4, we derive the nonequispaced fast Fourier transform by parametrization via the double covering 𝕊^2×𝕊^2 and developing the respective transform there. We also accomplish B) for the general cases (d) and _k,d. 











§ TWO INTRODUCTORY EXAMPLES


We first present a well-known elementary example on the interval [0,s], for which both aims A) and B) are achieved. Second, to support our perspective on discrepancy, we prove that the so-called 
Askey function is a discrepancy kernel of the form (<ref>).


 §.§ The Brownian motion kernel on [0,s]


Let dr be the Lebesgue measure on [0,∞). The mapping
h:[0,∞)→ℬ() defined by r↦ [r,∞) 
induces the pushforward measure β:=h_*( dr) that induces the discrepancy 

    𝒟_β(μ,ν_n)= ∫_0^∞ |μ([r,∞))-ν_n([r,∞))|^2  dr.

The associated discrepancy kernel K_β:×→ is[For r∈, we use the notation r_+=
r,    r≥ 0,
 0,   otherwise.] 

    K_β(x,y) = ∫_0^∞δ_x([r,∞))δ_y([r,∞))  dr =min(x,y)_+,

so that 𝒟_β(δ_x,δ_y)=|x-y| for x,y∈[0,∞). The restriction of the kernel K_β to [0,s]×[0,s] has the Fourier expansion 

    K_β(x,y) 
       =∑_m∈
     m  odd4s^2/m^2π^2·sin(π/2s mx)/√(s/2)·sin(π/2s my)/√(s/2),   x,y∈[0,s],

with respect to the Lebesgue measure σ_[0,s] on [0,s]. The reproducing kernel Hilbert space is

    ℋ_K_β([0,s])={f:[0,s]→ :  f  is absolutely continuous,  f(0)=0,  f'∈ L_2([0,s])},

where the inner product between f and g is given by ⟨ f',g'⟩_L_2([0,s]), cf. <cit.> and <cit.>. Note that K_β|_[0,1]×[0,1] is often called the Brownian motion kernel and ℋ_K_β([0,s]) is continuously embedded into the Sobolev space ℍ^1([0,s]).  





 §.§ Askey's function and its restrictions

Many positive definite kernels in the literature are of the form (<ref>) and, hence, are discrepancy kernels. For odd d, Askey's kernel function (x,y)↦ (1-x-y)^d+1/2_+ is positive definite, cf. <cit.>. In the following, we shall check that it is of the form (<ref>).

Denote the Euclidean ball of radius s centered at z∈^d by 

    𝔹^d_s(z):={x∈^d:x-z≤s},

with the conventions 𝔹^d_s:=𝔹^d_s(0) and 𝔹^d:=𝔹^d_1. Fix r>0 and consider the discrepancy

    𝒟_d,r(μ,ν_n):=1/(𝔹^d_r/2)
    ∫_^d| μ(𝔹^d_r/2(z))- ν_n(𝔹^d_r/2(z))|^2 d z,

where (𝔹^d_r/2)=π^d/2/Γ(d/2+1)(r/2)^d.  
The associated discrepancy kernel is 

    K_d,r(x,y)=1/(𝔹^d_r/2)∫_^dδ_x(𝔹^d_r/2(z))δ_y(𝔹^d_r/2(z))  d z.

In order to additionally integrate over r, recall the 
(generalized) hypergeometric functions

    klf_1,…,f_kg_1,…,g_lz:=∑_n=0^∞(f_1)_n ⋯ (f_k)_n/(g_1)_n ⋯ (g_l)_nz^n/n!,

where f_1,…,f_k,g_1,…,g_l,z∈ and (f)_n:=f· (f+1)⋯(f+n-1) denotes the Pochhammer symbol with (f)_0:=1. We consider G_d:[0,∞)→ given by

    G_d(r) = 21-d+1/4,-d-1/4-d/2r^2,    0≤ r≤ 1,
    
    0,    1<r.

Since d is odd, either d+1/4 or d-1/4 is a natural number, so that the series terminates and 
G_d is a polynomial in r^2 on [0,1]. 









By integration with respect to G_d, we obtain the L_2-discrepancy and the associated discrepancy kernel 

    𝒟_d(μ,ν_n):=
    
    ∫_0^∞𝒟_d,r(μ,ν_n) d G_d(r),   and    K_d(x,y)= ∫_0^∞K_d,r(x,y)  d G_d(r),

respectively. It turns out that K_d coincides with Askey's function.


Let d be odd. The discrepancy kernel K_d satisfies 

    K_d(x,y)=(1-x-y)^d+1/2_+,   x,y∈^d.


The proof is presented 
in Appendix <ref>. Provided that d≥ 3, Askey's kernel function  reproduces the Sobolev space ℍ^d+1/2(^d) with an equivalent norm, see <cit.>. 











§ THE DISTANCE KERNEL ON ^D-1

This section is dedicated to recall results on discrepancy kernels on the sphere ^d-1⊂^d, for d≥ 2, from <cit.> that shall guide our subsequent investigations.


Denote the geodesic ball of radius r centered at z∈^d-1 by 

    B^^d-1_r(z):= {x∈^d-1 : _^d-1(x,z)≤ r},

where _^d-1(x,z)=arccos(⟨ x,z⟩) is the geodesic distance on ^d-1. We define 

    h:[0,π]×^d-1→ℬ(^d),     (r,z)↦ B^^d-1_r(z)

and endow [0,π] with the weighted Lebesgue measure sin(r)dr, whereas ^d-1 carries the normalized, orthogonal invariant surface measure σ_^d-1. The push-forward β_d:=h_*(sin(r)dr ⊗σ_^d-1) is a measure on ℬ(^d), so that the associated L_2-discrepancy is

    𝒟_β_d(μ,ν_n)= ∫_0^π∫_^d-1 |μ(B^^d-1_r(z))-ν_n(B^^d-1_r(z))|^2  dσ_^d-1(z)sin(r) dr.

The associated discrepancy kernel is 

    K_β_d(x,y)=∫_0^π∫_^d-1δ_x(B^^d-1_r(z))δ_y(B^^d-1_r(z))    dσ_^d-1(z)sin(r) dr,   x,y∈^d.

According to <cit.>, see also <cit.>, K_β_d satisfies

    K_β_d(x,y)=1-Γ(d/2)/2√(π)Γ(d+1/2)x-y,     x,y∈^d-1.

If either x or y is not contained in ^d-1, then K_β_d(x,y)=0. 

Choose σ_:=σ_^d-1 for the decomposition (<ref>) and let {Y^m_l : l=1,…,Z(d,m)}⊂ L_2(^d-1,σ_^d-1) denote the set of orthonormal spherical harmonics of degree m on ^d-1, where 
Z(d,m)=2m+d-2/d-2m+d-3m. 
For τ>(d-1)/2, the Sobolev space ℍ^τ(^d-1) is the reproducing kernel Hilbert space associated with the reproducing kernel

    (x,y)↦∑_m=0^∞ (1+m(m+d-2))^-τ∑_l=1^Z(d,m)Y^m_l(x)Y^m_l(y),   x,y∈^d-1.

The coefficients in the Fourier expansion 

    1-Γ(d/2)/2√(π)Γ(d+1/2)x-y= ∑_m=0^∞ c_m ∑_l=1^Z(d,m)Y^m_l(x)Y^m_l(y),   x,y∈^d-1,

satisfy |c_m|∼ m^-d, cf. <cit.>. This is the same asymptotics as the coefficients in (<ref>) for s=d/2. Therefore, K_β_d|_^d-1×^d-1 reproduces the Sobolev space ℋ_β_d(^d-1)=ℍ^d/2(^d-1) with an equivalent 
norm[The kernel K_β_d|_^d-1×^d-1 generates an inner product ⟨·,·⟩_K_β_d|_^d-1×^d-1 in ℍ^d/2(^d-1), for which it is reproducing and the induced norm is equivalent to the standard norm in ℍ^d/2(^d-1), which is induced by the standard kernel (<ref>).]. 

In order to determine the Fourier coefficients of kernels on the sphere that are polynomial in x-y, such as K_β_d|_^d-1×^d-1, we require the Fourier coefficients of the monomial terms
x-y^p. For any p∈, the Fourier expansion 

    2^-p/2x-y^p = ∑_m=0^∞ a_m(p,^d-1) ∑_l=1^Z(d,m)Y^m_l(x)Y^m_l(y),   x,y∈^d-1,

holds with coefficients 
determined by 

    a_m(p,^d-1):=∬_^d-1×^d-1 2^-p/2x-y^p Y^m_l(x)Y^m_l(y) dσ_^d-1(x) dσ_^d-1(y).

Note that (<ref>) is well-defined for the entire range p>-(d-1) and p is not required to be an integer. For p>0, the following proposition is essentially due to <cit.>, see also <cit.>. Simple continuation arguments cover the full range of p, and the asymptotics  Γ(-p/2+m)/Γ(p/2+d-1+m)=m^-(p+d-1)(1+o(1)) are standard.

Suppose d≥ 2. For any p >-(d-1), we have

    a_m(p,^d-1)  
          =2^dΓ(d/2)/4√(π)·2^p/2Γ(d/2+p/2-1/2)/Γ(-p/2)·Γ(-p/2+m)/Γ(p/2+d-1+m).

In particular, if p∉2, then 

    |a_m(p,^d-1)|   =|2^dΓ(d/2)/4√(π)·2^p/2Γ(d/2+p/2-1/2)/Γ(-p/2)|m^-(p+d-1)(1+o(1)),

and the series (<ref>) terminates if p∈ 2. 

For p∈2, the term Γ(-p/2) is not well-defined and (<ref>) is to be understood with the convention Γ(-p/2+m)/Γ(-p/2)=(-p/2)_m. Hence, we observe a_m(p,^d-1)  = 0 for all m > p/2 if  p∈ 2. 

It is noteworthy that the kernel K_d,r in (<ref>) for d=3 is a discrepancy kernel that does not generate a Sobolev space on ^d but its restriction does. The proof of the following proposition is presented 
in Appendix <ref>.

Let r≥ 1. The reproducing kernel Hilbert space of K_3,r, given by (<ref>) with d=3, is continuously embedded into ℍ^2(^3), but the reverse embedding does not hold. In contrast, K_3,r|_^2×^2 reproduces ℍ^3/2(^2) with an equivalent norm.




To provide numerical examples for d=3, Proposition <ref> provides the coefficients (a_m)_m=0^∞ in the kernel expansion of K_β_3,

    1-1/4x-y= ∑_m=0^∞ a_m ∑_l=1^2m+1Y^m_l(x)Y^m_l(y),   x,y∈^2.

For (μ),(ν_n)⊂^2, the L_2-discrepancy (<ref>) for K_β_3 with 𝕏=^2  becomes

    𝒟_β_3(μ,ν_n) = ∑_m=0^∞ a_m ∑_l=1^2m+1|μ̂^m_l- 1/n∑_j=1^nY^m_l(x_j)|^2,

where μ̂^m_l denotes the Fourier coefficient of μ with respect to Y^m_l, cf. (<ref>).  By truncating this series, the nonequispaced fast Fourier transform on ^2, cf. <cit.>, enables efficient minimization of

    ∑_m=0^M a_m ∑_l=1^2m+1|μ̂^m_l- 1/n∑_j=1^nY^m_l(x_j)|^2

among all n-point sets {x_1,…,x_n}⊂^2 for fixed n. We are most interested in n≫ M. 
See Figure <ref>
for a numerical experiment with M=8 and n=50. 






§ DISCREPANCY KERNELS ON COMPACT SETS

Here we discuss discrepancy kernels that extend the kernels of the previous section in a natural way. 
For d≥ 1, let us define the half-space 

    Ω^d_r(z):={x∈^d : ⟨ z,x⟩≥ r }∈ℬ(^d),   z∈^d-1,   r∈.

For fixed s>0, we consider the mapping
h:[-s,s]×^d-1→ℬ(^d) defined by (r,z)↦Ω^d_r(z) 
and endow [-s,s] with the Lebesgue measure dr. The push-forward measure β_d,s:=h_*(dr⊗σ_^d-1) leads to the associated L_2-discrepancy

    𝒟_β_d,s(μ,ν_n)=∫_-s^s∫_^d-1| μ(Ω^d_r(z)) -  ν_n(Ω^d_r(z))|^2 dσ_^d-1(z)  d r.

The associated discrepancy kernel is

    K_β_d,s(x,y)=∫_-s^s∫_^d-1δ_x(Ω^d_r(z))δ_y(Ω^d_r(z))  dσ_^d-1(z)  d r,     x,y∈^d.

Since B^^d-1_r(z) = Ω^d_cos(r)(z)∩^d-1, for r∈[0,π] and z∈^d-1, 
we deduce 

    K_β_d|_^d-1×^d-1=K_β_d,1|_^d-1×^d-1,     d≥ 2,
 
with K_β_d as in (<ref>). 
In contrast to K_β_d, the kernel K_β_d,s is not identically zero outside of ^d-1×^d-1 and makes also sense for d=1. 

For d=1, we have ^0={± 1}, so that the half-spaces are Ω^1_r(1)=[r,∞) and Ω^1_r(-1)=(-∞,-r]. 
Direct calculation of (<ref>) yields 

    K_β_1,s(x,y)= 
    s-1/2|x-y|,    |x|,|y|≤ s,
    s/2+xy/2|y|,    |x|≤ s≤ |y|,
    
    s(xy),     s≤ |x|,|y|,


where  is the 
Heaviside step function. 








The Fourier expansion of the kernel K_β_1,s|_[-s,s]×[-s,s] with respect to the Lebesgue measure σ_[-s,s] on [-s,s] is 

    K_β_1,s(x,y)   = ∑_m∈
     m  odd4/m^2π^2·1/s·sin(π/2smx)sin(π/2smy) 
    
     + ∑_{u>0  : tan(u)=1/u}1/u^2·1/s(sin(u)^2+1)·cos(u/sx)cos(u/sy),   x,y∈[-s,s].

Its reproducing kernel Hilbert space is

    ℋ_K_β_1,s([-s,s])={f:[-s,s]→ :  f  is absolutely continuous,  f'∈ L_2([-s,s])},

where the inner product between f and g is given by 

    1/2s(f(-s)+f(s))(g(-s)+g(s))+⟨ f',g'⟩_L_2([-s,s]).


Note that ℋ_K_β_1,s([-s,s]) is continuously embedded into ℍ^1([-s,s]). The proof of Proposition <ref> is presented 
in Appendix <ref>.
It uses that, 

up to a constant, K_β_1,s|_[-s,s]×[-s,s] is the Green's function of the 1-dimensional harmonic equation Δ u = f on [-s,s] with the boundary conditions u'(s) = - u'(-s) and u(s) + u(-s) = -2u'(s). 


It turns out that K_β_d,s has a simple form on 𝔹^d_s×𝔹^d_s.

For d≥ 1, the discrepancy kernel K_β_d,s satisfies

    K_β_d,s(x,y) = s-Γ(d/2)/2√(π)Γ(d+1/2)x-y,     x,y∈𝔹^d_s.


The identity (<ref>) for x,y∈^d-1 with s=1 has been established in <cit.>, see also (<ref>). Essentially, 
the same proof still works for the more general situation. 







Theorem <ref> provides a simple form of K_β_d,s|_× with = 𝔹^d_s, which may facilitate further computations. An immediate consequence is 𝒟_β_d,s(δ_x,δ_y)=Γ(d/2)/√(π)Γ(d+1/2)x-y, for x,y∈𝔹^d_s. 











§ THE EUCLIDEAN BALL 𝔹^D

This section is dedicated to derive the Fourier expansion of the discrepancy kernel K_β_d,1 in (<ref>) on 𝔹^d. Proposition <ref> has covered d=1, and we now derive the spectral decomposition of 

    K^d,p:𝔹^d×𝔹^d →,     (x,y)↦x-y^p,

for all odd d with odd p>1-d with respect to the Lebesgue measure σ_𝔹^d on 𝔹^d. The case d=3 with p=-1 is discussed in <cit.>.

Let {𝒞^α_m: m∈, α>-1/2} denote the family of Gegenbauer polynomials with the standard normalization 

    𝒞^α_m(1) =m+2α-1m= Γ(m+2α)/Γ(2α)Γ(m+1),  α≠ 0.

By α=d/2-1, the addition theorem for spherical harmonics yields

    ∑_l=1^Z(d,m) Y^m_l(x)Y^m_l(y)
      =  2m+d-2d-2𝒞^d/2-1_m(⟨ x,y⟩),   x,y∈^d-1.

For m∈, let us define the kernels 𝒦^d,p_m:[0,1]×[0,1]→,

    𝒦_m^d,p(r,s) 
    
      = (-p/2)_m/(d/2 -1 )_m(min(r,s)max(r,s))^m max(r,s)^p 21m-p/2,1-d+p/2m+d/2(min(r,s)max(r,s))^2.


For d≥ 3 and arbitrary real p>1-d, we deduce from <cit.> that 


    x-y^p = ∑_m=0^∞𝒦^d,p_m(x,y)  𝒞^d/2-1_m(⟨xx,yy⟩),   x,y∈𝔹^d,  (x y  if  p < 0).


For x=0 or y=0, the 
right-hand side of (<ref>) is well-defined by 
analytic continuation. 


Using the addition formula (<ref>) we obtain

    K^d,p(x,y) = ∑_m=0^∞d-22m+d-2𝒦^d,p_m(x,y) ∑_l=1^Z(d,m)Y^m_l(xx) Y^m_l(yy),   x,y∈𝔹^d.

The Fourier expansion of 𝒦_m^d,p with respect to the measure r^d-1dr satisfies 

    𝒦_m^d,p(r,s) = ∑_j=1^∞λ_m,j^d,pφ_m,j^d,p(r) φ_m,j^d,p(s),
          ∫_0^1𝒦_m^d,p(r,s) φ_m,j^d,p(r)r^d-1d r = λ_m,j^d,pφ^d,p_m,j(s),

where ∫_0^1 |φ_m,j^d,p(r)|^2 r^d-1dr = 1. Then, by setting

    φ_m',j',l'^d,p(x):= φ_m',j'^d,p(x) Y^m'_l'(xx),     x ∈𝔹^d,

direct computations yield 

    T^d,pφ_m',j',l'^d,p(x)  = λ_m',j'^d,p(d-2) vol(𝕊^d-1)2m'+d-2φ_m',j',l'^d,p(x)

with the scaling ∫_𝔹^d| φ_m',j',l'^d,p(x)|^2d x

  = vol(𝕊^d-1). 


















This leads to the Fourier expansion

    K^d,p (x,y)
    
        = ∑_m=0^∞∑_j=1^∞λ_m,j^d,p(d-2)vol(𝕊^d-1)2m+d-2∑_l=1^Z(d,m)φ^d,p_m,j,l(x) /√(vol(𝕊^d-1))φ^d,p_m,j,l(y)/√(vol(𝕊^d-1))
    ,   x,y∈𝔹^d.






Thus, the original problem is reduced to the spectral decomposition of the sequence of kernels 𝒦^d,p_m, for m∈. The kernel 𝒦^d,p_m induces the integral operator

    T_m^d,p : L_2([0,1],r^d-1 dr) →𝒞([0,1]),   f↦∫_0^1 𝒦^d,p_m(·,r)f(r)r^d-1 dr,

with eigenvalues λ^d,p_m,j and eigenfunctions φ^d,p_m,j. 
We now specify these eigenvalues and eigenfunctions, where J_ν denotes the Bessel function of the first kind of order ν and 
ζ_k:= e^2π i/k is the k-th root of unity.

Suppose that both d≥ 3 and p>1-d are odd and let m∈. Then the following holds:


  a) Any eigenvalue λ≠ 0 of T_m^d,p is in a one-to-one correspondence with 

    ω =   |λ^-12^d+p-2(d+2m-2)(- p2)_d+p/2-1(d+p2-1)!|^1/d+p

with ω satisfying (A(ω))=0, where 

    A(ω) =
      (ζ_d+p^-iℓ J_m+d/2-i-1(ζ_d+p^ℓω)
          )_i=1, ℓ=0^d+p/2, d+p/2-1,     (- p2 )_d+p/2-1λ > 0,
    (ζ_2(d+p)^-i(2ℓ+1)J_m+d/2-i-1(ζ_2(d+p)^2ℓ+1ω)
          )_i=1, ℓ=0^d+p/2, d+p/2-1,     (- p2 )_d+p/2-1λ < 0.


  b) The eigenfunctions are exactly 

    r↦∑_ℓ=1^d+p/2 c_ℓ r^1-d/2J_m+d/2-1(ζ^ℓ-1_d+pω r ),    (- p2 )_d+p/2-1λ > 0,
    ∑_ℓ=1^d+p/2 c_ℓ r^1-d/2J_m+d/2-1(ζ^2ℓ-1_2(d+p)ω r ),    (- p2 )_d+p/2-1λ < 0,

where c∈^d+p/2 is in the nullspace of A(ω). 





Computer experiments seem to indicate that the nullspace of A(ω) is one-dimensional if (A(ω))=0. In that case, the function 

    r↦∑_ℓ=1^d+p/2 (-1)^ℓ A_[1,ℓ](ω)  r^1-d/2J_m+d/2-1(ζ^ℓ-1_d+pω r ),    (- p2 )_d+p/2-1λ > 0,
    ∑_ℓ=1^d+p/2 (-1)^ℓ A_[1,ℓ](ω)  r^1-d/2J_m+d/2-1(ζ^2ℓ-1_2(d+p)ω r ),    (- p2 )_d+p/2-1λ < 0,

where A_[1,ℓ](ω) denotes the (1,ℓ) minor of A(ω), spans the eigenspace  associated 
with λ.
 

Appendix <ref> is dedicated to the proof of Theorem <ref>. The proof reveals strong ties with polyharmonic operators on the unit ball and higher order differential operators on the interval [0,1]. We refer to <cit.> for structurally related spectral decompositions of polyharmonic operators on [0,1] with homogeneous Neumann boundary conditions.   
















The nonzero eigenvalues of T^3,1_m for m∈∖{0} are exactly the positive solutions of the equation
      
    J_m-1/2(ω) J_m-3/2( iω) - i J_m-1/2( iω)  J_m-3/2(ω)
          = 0,

with λ=-ω^-4(4m+2). The corresponding eigenspaces are 1-dimensional with the representative
 
    f_λ(r) 
           = r^-1/2(J_m+1/2(ω r)i J_m-1/2(iω) +  J_m+1/2(iω r)J_m-1/2(ω)).







The formulas in Corollary <ref> are derived from Theorem <ref>. Since J_m-1/2( iω)≠ 0, for all ω∈∖{0}, the eigenspaces are 1-dimensional and (<ref>) is not the zero-function.  


In view of (<ref>) in Theorem <ref> we are particularly interested in kernels of the form c-x-y. In this case, the expansion holds with -𝒦^d,1_m for m≥ 1 and c-𝒦^d,1_0 for m=0. 








































§ THE ROTATION GROUP (3)

In this section we derive the Fourier expansion of the discrepancy kernel on the special orthogonal group (3). The eigenfunctions turn out to be classical functions but the coefficients and their decay rates need to be determined. We 
also provide numerical experiments by using the nonequispaced fast Fourier transform on (3).


 §.§ Fourier expansion on (3)

By identifying ^d× d with ^d^2, Theorem <ref> applies to subsets of ^d× d endowed with the trace inner product 

    ⟨ x,y⟩_:=(x^⊤ y),   x,y∈^d× d,

and the induced Frobenius norm ·_ on ^d× d. 
In this way, (3) is contained in 𝔹^9_√(3), and it is natural to consider s=√(3). We endow (3) with the normalized Haar measure σ_(3). Let {𝒟^m_k,l:k,l=-m,…,m} denote the Wigner 𝒟-functions on (3), which are closely related to the irreducible representations of (3) and provide an orthonormal basis for L_2((3)), cf. <cit.>. For p>0, the Fourier expansion

    2^-p/2x-y^p_ =∑_m=0^∞ a_m(p,(3)) ∑_k,l=-m^m 𝒟^m_k,l(x)𝒟^m_k,l(y),   x,y∈(3),

holds and, analogous to (<ref>), the coefficients are determined by
 
    a_m(p,(3)):=∬_(3)×(3)2^-p/2x-y_^p 𝒟^m_k,l(x)𝒟^m_k,l(y) dσ_(3)(x) dσ_(3)(y).

We now compute these coefficients for the entire range p>-3.

For p > -3, the coefficients (<ref>) are 

    a_m(p,(3)) =   2^pΓ(p/2+3/2)/√(π)Γ(-p/2)·Γ(-p/2+m)/Γ(p/2+2+m)·1/(m+1/2).

In particular, if p∉2, then 

    |a_m(p,(3))| = | 2^pΓ(p/2+3/2)/√(π)Γ(-p/2)|m^-(p+3)(1+o(1)),     m∈,

and the series (<ref>) terminates if p∈ 2. 

The proof is given in Appendix <ref>. For p∈ 2, we again apply the convention Γ(-p/2+m)/Γ(-p/2)=(-p/2)_m in (<ref>), so that a_m(p,(3))=0 for all m>p/2 if p∈ 2. Provided that τ>3/2, the Sobolev space ℍ^τ((3)) is the reproducing kernel Hilbert space associated with the reproducing kernel

    (x,y)↦∑_m=0^∞ (1+m(m+1))^-τ∑_k,l=-m^m 𝒟^m_k,l(x)𝒟^m_k,l(y),   x,y∈(3).

The choice p=1 in Proposition <ref> implies that the kernel K_β_9,s|_(3)×(3) reproduces the Sobolev space ℋ_K_β_9,s((3))=ℍ^2((3)) with an equivalent norm provided that s≥√(3).



For d≥ 2, K_β_d^2,√(d)|_(d)×(d) reproduces the Sobolev space ℍ^d(d-1)+2/4((d)) with an equivalent norm. Indeed, Theorem <ref> and Section <ref> yield that K_β_d^2,d|_^d^2-1×^d^2-1 reproduces ℍ^d^2/2(^d^2-1) with equivalent norms. Rescaling implies that K_β_d^2,d|_d ^d^2-1× d ^d^2-1 reproduces ℍ^d^2/2(d ^d^2-1). Since (d)⊂ d ^d^2-1, the assertion for (d) follows from results on restricting kernels in <cit.>.  




 §.§ Numerical examples on (3)

Proposition <ref> yields the coefficients of the kernel expansion

    K_β_9,√(3)(x,y) = ∑_m=0^∞ a_m ∑_k,l=-m^m 𝒟^m_k,l(x)𝒟^m_k,l(y),   x,y∈(3).

For (μ),(ν_n)⊂(3), the L_2-discrepancy (<ref>) for K_β_9,√(3) becomes

    𝒟_β_9,√(3)(μ,ν_n) = ∑_m=0^∞ a_m ∑_k,l=-m^m  |μ̂^m_k,l- 1/n∑_j=1^n 𝒟^m_k,l(x_j)|^2,

where μ̂^m_k,l denotes the Fourier coefficient of μ with respect to 𝒟^m_k,l, cf. (<ref>). We 
truncate the series (<ref>) at M=8 and minimize 

    ∑_m=0^M a_m ∑_k,l=-m^m  |μ̂^m_k,l - 1/n∑_j=1^n 𝒟^m_k,l(x_j)|^2

among all n-point sets {x_1,…,x_n}⊂(3) for fixed n=30. We efficiently solve the least squares minimization by using the nonequispaced fast Fourier transform on (3), cf. <cit.>. 

Figure <ref> shows the minimizing points mapped onto 𝔹^3. 






§ THE GRASSMANNIAN _2,4

First, the Fourier expansion of the discrepancy kernel on _2,4 is computed. To prepare for developing the nonequispaced fast Fourier transform on _2,4, we then explicitly parametrize the Grassmannian _2,4 by its double covering ^2×^2. Next, we derive the nonequispaced fast Fourier transform on ^2×^2 and provide numerical minimization experiments on _2,4. 


 §.§ Fourier expansion on _2,4

Proposition <ref> also applies to the Grassmannian

    _2,4:={x∈^4× 4 : x^⊤ = x,  x^2=x, (x)=2}

with s=√(2) when ^4× 4 is identified with ^16. To derive the Fourier expansion on _2,4, we require some preparations. We shall use integer partitions λ=(λ_1,λ_2)∈^2 with λ_1≥λ_2≥ 0. We also denote |λ|:=λ_1+λ_2. The orthogonal group (4) acts transitively on _2,4 by conjugation and induces the irreducible decomposition 

    L_2(_2,4,σ__2,4) = ⊕_λ_1≥λ_2≥ 0 H_λ(_2,4),    
    H_λ(_2,4) ⊥ H_λ'(_2,4),   λλ',

where σ__2,4 is the normalized orthogonally invariant measure on _2,4 and H_λ(_2,4) is equivalent to the
irreducible representation ℋ_2λ^4 of (4) with type 
2λ, cf. <cit.>. The normalized eigenfunctions of the Laplace–Beltrami operator on _2,4 form an orthonormal basis for L_2(_2,4), and each H_λ(_k,d) is contained in the eigenspace E_α_λ associated 
with the eigenvalue α_λ = 4(λ_1^2+λ_2^2+λ_1), cf. <cit.>. 


































Let Q_λ be the reproducing kernel of H_λ(_2,4). Any orthonormal basis {φ_λ,l}_l=1^(ℋ^4_2λ) for H_λ(_2,4) yields the spectral decomposition 

    Q_λ(x,y) = ∑_l=1^(ℋ^4_2λ)φ_λ,l(x)φ_λ,l(y),     x,y∈_2,4,  λ_1≥λ_2≥ 0.

The orthogonal decomposition (<ref>) leads to the Fourier expansion 

    2^-p/2 x -y_^p=∑_λ_1≥λ_2≥ 0a_λ(p,_2,4) Q_λ(x,y),     x,y ∈_2,4,   p>0.

The coefficients a_λ(p,_2,4) in (<ref>) are defined by

    a_λ(p,_2,4) := 1/(ℋ^4_2λ)∬__2,4×_2,42^-p/2 x -y_^p Q_λ(x,y)  dσ__2,4(x) dσ__2,4(y).

In order to determine a_λ(p,_2,4), we shall make use of the hypergeometric coefficients (f)_(λ_1,λ_2) := (f)_λ_1 (f-12)_λ_2.

For p>-4, 
we have

    a_λ(p,_2,4) =   4^-|λ||λ|!/(32)_|λ| (32)_λ(-p2)_|λ| 43|λ|+12,|λ|+22,|λ|2-p4,|λ|+12-p4|λ|+32,λ_1+32,λ_2+11.

In particular, if p∉2, then 

    |a_λ(p,_2,4)|= |Γ(p/2+2)/2Γ(-p/2)|λ^-(p+4) (1+o(1)),
        λ_1≥λ_2≥0,

and the series (<ref>) terminates if p∈ 2. 

The proof of this theorem is contained in Section <ref> of
Appendix <ref>. 

If p∈ 2, then a_λ(p,_2,4)=0 for all |λ|>p/2. For τ>2, the Sobolev space ℍ^τ(_2,4) is the reproducing kernel Hilbert space with associated reproducing kernel

    (x,y)↦∑_λ_1≥λ_2≥ 0(1+4λ_1^2+4λ_2^2+4λ_1)^-τQ_λ(x,y),   x,y∈_2,4,

cf. <cit.>. 
 Since the coefficients in (<ref>) behave asymptotically as λ^-2τ, the choice p=1 in Theorem <ref> implies that the kernel K_β_16,s|__2,4×_2,4 reproduces the Sobolev space ℋ_K_β_16,s(_2,4)=ℍ^5/2(_2,4) with an equivalent norm provided that s≥√(2). 


By invoking <cit.>, we deduce that, for d≥ 2, K_β_d^2,√(k)|__k,d×_k,d reproduces the Sobolev 
space ℍ^k(d-k)+1/2(_k,d) with an equivalent norm. 





 §.§ Parametrization of _2,4 by (𝕊^2×𝕊^2)± 1




To derive the nonequispaced fastFourier transform on _2,4, we shall first explicitly construct the parametrization of _2,4 by its double covering ^2×^2. We denote the d× d-identity matrix by I_d, and the cross-product between two vectors x,y∈^2 is denoted by x× y∈^3. The mapping 𝒫:^2×^2→_2,4 given by

    (x,y) ↦
    
    
    
    
    
    
    
    
           1/2[                     1 + x^⊤ y                   - (x × y)^⊤;                       - x × y                               x y^⊤ + y x^⊤ + (1-x^⊤ y) I_3;                               ]

is surjective and, for all x,y,u,v∈𝕊^2,

    𝒫(u,v)=𝒫(x,y)   if and only if    (u,v) ∈{± (x,y) },

see Section <ref> and Theorem <ref> 
of Appendix <ref>. 

In order to specify the inverse map, note that 𝕊^2×𝕊^2± 1 
can be identified with
ℳ(3):={xy^⊤∈ℝ^3× 3 : x,y∈𝕊^2}. 
We now define ℒ:𝒢_2,4→ℳ(3),

    ℒ(P) :=     ([ 1/2(P_11+P_22-P_33-P_44)              P_23 - P_14              P_24 + P_13;              P_23 + P_14 1/2(P_11-P_22+P_33-P_44)               P_34- P_12;              P_24 - P_13              P_34 + P_12 1/2(P_11-P_22-P_33+P_44);                          ]),

and direct computations lead to

    ℒ(𝒫(x,y)) = xy^⊤,   x,y∈𝕊^2.

The right-hand side determines x and y up to the ambiguity (<ref>). Under the Frobenius norm, 𝒫 is distance preserving in the sense 


    𝒫(x,y)-𝒫(u,v)_ = xy^⊤ - uv^⊤_,   x,y,u,v∈𝕊^2.

The latter follows from (<ref>) in Lemma <ref> in
Section <ref> of Appendix <ref>. 



We shall now check how the spherical harmonics Y^m_l on ^2 relate to the eigenfunctions φ_λ,l∈ H_λ(_2,4) of the Laplace–Beltrami operator on _2,4, cf. (<ref>). The functions Y_k,l^m,n:𝒢_2,4→ℂ given by
  
    Y_k,l^m,n(𝒫(x,y)) := Y_k^m(x) · Y_l^n(y)

are well-defined for m+n ∈ 2ℕ, the latter taking into account the ambiguity (<ref>). 

For m_λ:=(λ_1+λ_2) and n_λ:=(λ_1-λ_2), we have

    H_λ(𝒢_2,4) = span{ Y_k,l^m_λ,n_λ, Y_l,k^n_λ,m_λ :
      k=-m_λ,…,m_λ,  l=-n_λ,…,n_λ}.


The proof is presented 
at the end of Section <ref> of Appendix <ref>. 

Note that the geodesic distance on _2,4 is __2,4(P,Q)=√(2)√(θ_1^2+θ_2^2), where θ_1,θ_2∈[0,π/2] are the principal angles determined by the two largest eigenvalues cos^2(θ_1) and cos^2(θ_2) of the matrix PQ. 
Aside from (<ref>), 𝒫 is also distance-preserving with respect to the respective geodesic distances, i.e.,[The geodesic distance on ^2 induces the geodesic distance on 𝕊^2×𝕊^2± 1 by

    ^2_𝕊^2×𝕊^2± 1(([ x; y ]),([ u; v ]))=min{^2_^2(x,u)+^2_^2(y,v),^2_^2(-x,u)+^2_^2(-y,v)}.

] 

    __2,4(𝒫(x,y),𝒫(u,v)) = _𝕊^2×𝕊^2± 1(([ x; y ]),([ u; v ])),   x,y,u,v∈^2.
 
This equality follows from (<ref>) in Lemma <ref> in the appendix 
via further direct calculations. 



The identity (<ref>) provides explicit expressions for the orthonormal basis {φ_λ,l}_l=1^(ℋ^4_2λ) of H_λ(_2,4) that is used to construct the reproducing kernel Q_λ in (<ref>). It also provides a fast Fourier transform on _2,4 from the respective transform on ^2×^2 that is developed in the subsequent section.












 §.§ Nonequispaced Fast Fourier Transform on _2,4

The nonequispaced fast (spherical) Fourier transform on ^2 has been developed in <cit.> under the acronym 
. Here, we shall derive the analogous transform on ^2×^2, which induces the nonequispaced fast Fourier transform on _2,4 via the mapping 𝒫 and (<ref>) with (<ref>).

For a given finite set of coefficients
f^m_1,m_2_k,l∈ℂ, m_1,m_2=0,…,M, k=-m_1,…,m_1, l=-m_2,…,m_2, 
we aim to evaluate

    F(x,y):=∑_m_1,m_2=0^M ∑_k=-m_1^m_1∑_l=-m_2^m_2 f^m_1,m_2_k,l Y^m_1_k(x)Y^m_2_l(y)

at n scattered locations (x_j,y_j)_j=1^n⊂^2×^2. 





Direct evaluation of (<ref>) leads to O(n M^4) operations. We shall now derive an approximative algorithm that is more efficient for n≫ M. 

By following the ideas in <cit.>, switching to spherical coordinates reveals that (<ref>) is a 4-dimensional trigonometric polynomial. This enables the use of the 4-dimensional nonequispaced fast Fourier transform  to significantly reduce the complexity. In spherical coordinates the spherical harmonics are trigonometric polynomials such that









    Y^m_k(z(θ,φ)) =e^ ikφ∑_k'=-m^m c^m_k,k'e^ ik'θ,     z(θ,φ) = ([ sin(θ)cos(φ); sin(θ)sin(φ);       cos(θ) ])
    ∈^2,

where 0≤θ≤π, 0≤φ≤ 2π, and c^m_k,k'∈ℂ are suitable coefficients that we assume to be given or precomputed. Thus, for
x = z(θ_1,φ_1) and y=z(θ_2,φ_2),  
there are coefficients b^k,l_k',l'∈ℂ such that

    F(x,y)=∑_k,l,k',l'=-M^M b^k,l_k',l' e^ ikφ_1e^ ilφ_2 e^ i k'θ_1e^ i l'θ_2.

We check 
in Section <ref> of Appendix <ref> that the set of coefficients b^k,l_k',l' can be computed by O(M^5) operations provided that the numbers  c^m_k,k' in (<ref>) are given. The expression (<ref>) can be evaluated at n scattered locations by the nonequispaced discrete Fourier transform  with O(n M^4) operations, cf. <cit.>. An efficient  approximative algorithm is the nonequispaced fast Fourier transform  that requires only O(M^4log(M)+n|log(ϵ)|^4) operations with accuracy ϵ, see <cit.> for details on accuracy. Thus, our algorithm for evaluating (<ref>) at n scattered locations requires O(M^5+n|log(ϵ)|^4) operations. For n≫ M, this is a significant reduction in complexity compared to the original O(nM^4) operations. We shall choose n∼ M^4 in the subsequent section, so that the complexity is reduced from O(M^8) to O(M^5+M^4|log(ϵ)|^4) operations. For potential further reduction, we refer to Remark <ref> in the appendix. 










 §.§ Numerical example on _2,4

By Theorem <ref>, we can calculate the coefficients of the kernel expansion

    K_β_16,√(2)(x,y) = ∑_λ_1≥λ_2≥ 0 a_λ∑_l=1^(ℋ^d_2λ)φ_λ,l(x)φ_λ,l(y),   x,y∈𝒢_2,4.

The eigenfunctions φ_λ,l are given by the tensor products of spherical harmonics in (<ref>), cf. Theorem <ref>. For (μ),(ν_n)⊂_2,4, the L_2-discrepancy (<ref>) of the kernel K_β_16,√(2)|__2,4×_2,4 is 

    𝒟_β_16,√(2)(μ,ν_n) = ∑_λ_1≥λ_2≥ 0 a_λ∑_l=1^(ℋ^d_2λ)|μ̂_λ,l- 1/n∑_j=1^n φ_λ,l(x_j)|^2,

where μ̂_λ,l is the Fourier coefficient of μ with respect to φ_λ,l, cf. (<ref>).  

Let us consider μ=σ__2,4. 
According to <cit.> (see also <cit.>), 
the lower bound 

    n^-5/4≲𝒟_β_16,√(2)(μ,ν_n)

holds for all n-point sets {x_1,…,x_n}⊂_2,4. We truncate the series (<ref>) and let ν^M_n=1/n∑_j=1^nδ_x^M_j denote a minimizer of 

    ∑_λ+λ_2≤ M a_λ∑_l=1^(ℋ^d_2λ)|μ̂_λ,l - 1/n∑_j=1^n φ_λ,l(x_j)|^2

among all n-point sets {x_1,…,x_n}⊂_2,4. A suitable choice n∼ M^4 leads to the optimal rate 

    𝒟_β_16,√(2)(μ,ν^M_n)∼  
    n^-5/4,

cf. <cit.>. 
Note that we can efficiently solve  the least squares minimization (<ref>) by using the nonequispaced fast Fourier transform on _2,4 derived from the nonequispaced fast Fourier transform on ^2×^2 of Section <ref> and applying Theorem  <ref>. Figure <ref> shows logarithmic plots of the number of points versus the L_2-discrepancy. We observe a line with slope -5/4 as predicted by (<ref>). 

























































































































































§ ACKNOWLEDGEMENT

ME and MG have been funded by the Vienna Science and Technology Fund (WWTF) through project VRG12-009. The research of CK has been partially supported by the Austrian Science Foundation (FWF) 
grant SFB F50 in the framework of the
Special Research Program
“Algorithmic and Enumerative Combinatorics". ME would like to thank Christine Bachoc, Karlheinz Gröchenig, Andreas Klotz, and Gerald Teschl for helpful discussions.




§ PROOFS FOR SECTION <REF>



Let h_d:[0,∞)→ denote Euclid's hat function given by 

    h_d(x)=1/(𝔹^d_1/2)(1_𝔹^d_1/2 * 1_𝔹^d_1/2)(x),    x∈^d,

where 1_𝔹^d_1/2 is the indicator function of 𝔹^d_1/2. 
For x,y∈^d and t=x-y, we derive

    K_d(x,y)     = ∫_0^∞1/(𝔹^d_r/2)∫_^d 1_𝔹^d_r/2(z)(x)1_𝔹^d_r/2(z)(y)  d z   dG_d(r)  =∫_t^∞ h_d(t/r) dG_d(r)
        = ∫_t^1 h_d(t/r)G_d'(r)  dr -h_d(t)G_d(1) =t∫_t^1 r^-2h_d'(t/r)G_d(r)  dr,
    where the last equality is due to partial integration. 
    An explicit expression 
    for h'_d is stated in <cit.>, so that we obtain
    K_d(x,y) 
        =d Γ(d/2)/√(π)Γ(d+1/2) t ∫_1^1/t (1-r^2t^2)^d-1/221-d+1/4,-d-1/4-d/2r^-2 dr 
        = d Γ(d/2)/√(π)Γ(d+1/2) t  ∑_j=0^d-1/2(-1)^jd-1/2jt^2j∑_l=0^d+1/4(-d+1/4)_l (-d-1/4)_l/(-d/2)_l l!∫_1^1/t r^2(j-l) dr 
    
    
        = d Γ(d/2)/√(π)Γ(d+1/2)∑_l=0^d+1/4 t^2l(-d+1/4)_l (-d-1/4)_l/(-d/2)_l l!∑_j=0^d-1/2(-1)^jd-1/2j/2j-2l+1
               - d Γ(d/2)/√(π)Γ(d+1/2)∑_j=0^d-1/2(-1)^jd-1/2jt^2j+1∑_l=0^d+1/4(-d+1/4)_l (-d-1/4)_l/(-d/2)_l l! (2j-2l+1)  .

We now compare coefficients of powers of t with those of the polynomial (1-t)^d+1/2. In order to check the coefficient of t^2l, we first observe 

    ∑_j=0^d-1/2(-1)^jd-1/2j/2j-2l+1 =
    21-d-1/2,1/2-l3/2-l11/1-2l
    
    =Γ(d+1/2)  Γ(1/2-l)/2Γ(d/2-l+1),

where the last equality makes use of Gauss' Theorem for the hypergeometric series
evaluated at 1 (cf. <cit.>). 
Direct computation yields

    (-d+1/4)_l (-d-1/4)_l/(-d/2)_l
    =
    
    d+1/22l(2l)!/(-4)^lΓ(d/2-l+1)/Γ(d/2+1),
 
so that the coefficient of t^2l in K_d(x,y) is 

    d+1/22l(2l)!/(-4)^l l!Γ(1/2-l) /Γ(d/2+1)  2d Γ(d/2)/√(π) = d+1/22l(2l)!/(-4)^l l!Γ(1/2-l) /√(π)
    =d+1/22l.

Hence, the coefficients for even powers of t match. To check the coefficient of t^2j+1, we first assume d-1/4∈. 
The Pfaff–Saalschütz Theorem (cf. <cit.>) yields

    ∑_l=0^d+1/4(-d+1/4)_l (-d-1/4)_l/(-d/2)_l l! (2j-2l+1)    = 32-d+1/4,-d-1/4,-1/2-j-d/2,1/2-j11/1+2j
       =
    (-d-1/4)_d-1/4(j-d-1/2)_d-1/4/(-d/2)_d-1/4(j-d-3/4)_d-1/41/1+2j.

Thus, the coefficient of t^2j+1 in K_d(x,y) is nonzero if and only if j≤d-1/4. Moreover, it is given by 

    -d Γ(d/2)(-1)^jd-1/2j(-d-1/4)_d-1/4(j-d-1/2)_d-1/4/√(π)Γ(d+1/2)(1+2j)(-d/2)_d-1/4(j-d-3/4)_d-1/4
    =-d/√(π) Γ(d/2)(-1)^j (d-1/4)!(d+3/4-j)_d-1/4/(d-1/2-j)! j!(1+2j)(d+1/4+1)_d-1/4(1/2-j)_d-1/4.

Further computations using the duplication formula eventually lead to -d+1/22l+1. The case d+1/4∈ is checked analogously, and we conclude the proof. 














§ PROOFS FOR SECTION <REF>


































By expressing the ^d-Fourier transform 1_𝔹^d_r(ξ) in terms of the Bessel function of the first kind of order d/2 and using its asymptotics, we deduce


    |1_𝔹^d_r(ξ)| ≲(1+ξ^2)^-(d+1)/4,   ξ∈^d.

Due to the 
zeros of the Bessel function, the respective lower bound cannot hold. This implies the embedding claims for ℍ^d+1/2(^d),  in particular, for d=3. 

To address the restricted kernel, we observe that the value K_d,r(x,y) coincides with the ^d-volume of the two intersecting balls 𝔹^d_r(x)∩𝔹^d_r(y). This volume has been explicitly computed in <cit.>. For d=3, we obtain

    K_3,r(x,y)=π(2r-x-y)^2(x-y+4r)/12,   x,y∈^2,

so that K_3,r|_^2×^2 is a polynomial of degree 3 in x-y. Its Fourier coefficients (a_m)_m∈ are linear combinations of the Fourier coefficients of the monomial terms, so that (<ref>) implies a_m≲ m^-3. 

We have checked that there are no cancelations in these linear combinations. Therefore, the asymptotics (<ref>) also imply the associated bound from below, which leads to a_m ∼ m^-3. Thus, K_3,r|_^2×^2 reproduces the Sobolev space ℍ^3/2(^2) with equivalent norms. 








§ PROOFS FOR SECTION <REF>



Let us consider s=1. The general case follows from rescaling. In order to determine the Fourier expansion of K_β_1,1|_[-1,1]×[-1,1], we need to determine the eigenfunctions with respect to the positive eigenvalues of the integral operator

    T : L_2([-1,1])→𝒞([-1,1]),   f↦∫_-1^1 f(t)(1-12|·-t|)  dt.

For λ>0, the equation Tϕ=λϕ is equivalent to the associated Sturm–Liouville eigenvalue problem. Indeed, 
differentiating twice on both sides and carrying out a short calculation,
we arrive at the second order homogeneous differential equation 

    ϕ”(x) +1/λϕ(x)=0,   x∈(-1,1),

with boundary conditions 
ϕ'(1) = - ϕ'(-1) and 

ϕ(1) + ϕ(-1) = -2ϕ'(1). 
The general solution of (<ref>) is 
ϕ(x)=c_1cos(x/√(λ))+c_2sin(x/√(λ)). 
Direct calculations using the boundary conditions determine c_1,c_2 and λ, so that Mercer's Theorem and normalization of the eigenfunctions provide the claimed Fourier expansion of the kernel. 

Computations analogous to <cit.> and <cit.> verify the claimed form of the reproducing kernel Hilbert space.
















§ PROOFS FOR SECTION <REF>

The two linearly independent eigenfunctions of the differential operator 

    D_m := ∂^2_r + d-1r∂_r - m(m+d-2)r^2

on [0,1] with respect to a possibly complex eigenvalue -ω^2 are

    𝒥^d,ω_m(r):=J_m+d/2-1(ω r )/r^d/2-1,  and   𝒴^d,ω_m(r):= Y_m+d/2-1(ω r )/r^d/2-1,

where J_ν and Y_ν are the Bessel functions of first and second type, respectively. 

For odd d≥ 3, odd p>1-d, and m∈, any eigenfunction of T_m^d,p with eigenvalue λ≠ 0 is a linear combination of {𝒥^d,ω_ℓ_m: ℓ=1,…,d+p/2}, where

    ω_ℓ =   |v|^1/d+p e^π i2ℓ/d+p,     (-1)^d+p/2v > 0,
     e^π i2ℓ+1/d+p,     (-1)^d+p/2v < 0,

and

    v=-λ^-1 2^d+p-2(d+2m-2)(2-d2)_d+p/2-1(d+p2-1)!.










Up to a constant depending on d and p, K^d,p is the Green's function of the polyharmonic equation Δ^d+p/2u=f on 𝔹^d with certain nonlocal boundary conditions, cf. <cit.>. In particular and by specifying the constant, one deduces that any eigenfunction of 

    T^d,p : L_2(𝔹^d) →𝒞(𝔹^d),   f↦∫_𝔹^d K^d,p(·,y)f(y) d y

with eigenvalue λ̃≠ 0 is an eigenfunction of Δ^d+p/2 with eigenvalue 

    - λ̃^-1 2^d+p-2(2-d2)_d+p/2-1(d+p2-1)! (d-2)(^d-1).



The Laplacian in polar coordinates is Δ = ∂^2_r + d-1/r∂_r + 1/r^2Δ_𝕊^d-1. The decomposition (<ref>) (see also (<ref>)) 

yields Δ^d+p/2φ^d,p_m,j,l=(D_m^d+p/2φ^d,p_m,j)Y^m_l on 𝔹^d, where D_m is as in (<ref>). Since λ̃= (d-2)vol(𝕊^d-1)/2m+d-2λ, any eigenfunction of T^d,p_m with eigenvalue λ≠ 0 is an eigenfunction of D^d+p/2_m with eigenvalue 

    v=-λ^-1 2^d+p-2(d+2m-2)(2-d2)_d+p/2-1(d+p2-1)!.

The linearly independent eigenfunctions of D_m^d+p/2 with respect to any eigenvalue v  0 are 

    𝒥^d,ω_ℓ_m,  and   𝒴^d,ω_ℓ_m,    ℓ=1,…,d+p2,

where

    (-ω_ℓ^2)^d+p/2 = (-1)^d+p/2ω_ℓ^d+p = v,     ℓ=1,…,d+p2,

and we take

    ω_ℓ := |v|^1/d+p e^2π iℓ/d+p,     (-1)^d+p/2v > 0
     e^π i2ℓ+1/d+p,     (-1)^d+p/2v < 0,    
          ℓ=1,…,d+p/2.

As an eigenfunction of a positive integer power of the Laplacian, ϕ^d,p_m,j,l is real analytic on 𝔹^d, cf. <cit.>. Hence, the radial part φ^d,p_m,j,l must be an analytic function on [0,1] with even or odd parity for m even or odd, respectively, cf. <cit.>. 

The functions 𝒴^d,ω_ℓ_m do not have matching parity, which concludes the proof. 


In order to identify the eigenvalues and the linear combination in Lemma <ref>, we check how T_m^d,p in (<ref>) acts on 𝒥^d,ω_m.

For d and p>1-d odd, we have

    (T_m^d,p𝒥^d,ω_m)(r) =  ω^-(d+p) 2^d+p-2(d+2m-2)(- p2)_d+p/2-1(d+p2-1)! 𝒥^d,ω_m(r) 
    
    -   (-p/2)_m/(d/2 -1)_m∑_i=1^d+p/2  2^i-1(d+p2-i+1)_i-1ω^-i𝒥^d,ω_m-i(1) r^m 21i-d+p2, m- p2 d2 + mr^2.



The idea of the proof is to use the series expansion of the Bessel functions, apply the integral operator to each term, and eventually recover the right-hand side of (<ref>). We shall provide the skeleton of the proof and omit some lengthy computations.

Let d ≥ 1, p ≥ 1-d and r≥ s. If d+p is even, direct computations yield

    𝒦_m^d,p(r,s)  = (-p/2)_m/(d/2 - 1)_m∑_k=0^d+p/2-1 ( -1)^k  d+p/2-1 k(-p/2+m)_k/(d/2 + m)_k s^m+2k r^p-m-2k.

We obtain that ( T_m^d,2l-d𝒥^d,ω_m)(r) equals
 
    (d/2 -l)_m/(d/2 -1)_m∑_k=0^l-1l-1 k(d/2-l+m)_k/( -1)^k(d/2 + m)_k[ ∫_0^r J_m+d/2-1(ω s ) s^d/2+m+2kd s/r^2l-d-m-2k+ ∫_r^1 J_m+d/2-1(ω s ) s^2l-d/2-m-2kds/r^m+2k].

    For α > 0, k ∈ℕ_0, ω∈ℂ and r > 0, 
integration of each term of the power series of the Bessel function eventually yields

    ∫_0^r J_α(ω s) s^α + 2k +1d s     =  ∑_i=0^k2^i (-k)_i/ω^i+1J_α + i+1(ω r)/r^i-α -2k-1,
    
    ∫_r^1 J_α(ω s) s^-α + 2k +1d s    =  ∑_i=0^k(-2)^i (-k)_i/ω^i+1(J_α-i-1(ω r)/r^i+α-2k-1 -J_α-i-1(ω )),

  which follows from direct computations and 
   
    ∑_j=0^∞(-k)_n+j-1/(α+n)_j+1 =   (-k)_n-1/α+n21n-k-1,1α+n+11=
    (-k)_n-1/α+k+1,

for α∈ℝ, k ∈ℕ_0, and n ∈ℕ with α∉{-(k+1),-k,…,-n }. 
By applying (<ref>) and (<ref>), 
we can express ( T_m^d,2l-d𝒥^d,ω_m)(r) as a sum of Bessel functions of various different orders. 
Straightforward but lengthy computations combined with the identity

    J_α(x) = ∑_i,k=0^l(x/2)^2l+1-i(-1)^k(-k)_i/(α-l)_l+1(l-k)! k![  (α-l)_kJ_α+i+1(x)/(α+1)_k
    
    
    +
           (α-l)_l-k J_α-i-1(x)/(-1)^l+i(α+1)_l-k]
 
eventually lead to the claimed equality (<ref>). 

In order to verify (<ref>),
according to the definition of the Bessel function J_α(x),
we have to show that the coefficient
of (x/2)^α+2m on the right-hand side of (<ref>) equals
(-1)^m/(m!Γ(m+α+1)), for m=0,1,…. Let first m≥ l.
Then this coefficient R(α,m) equals

    R(α,m)    :=
    ∑_i,k=0^l(-1)^k+m-l-1(-k)_i/(α+k-l)_l+1(l-k)! k!(m-l-1)!Γ(m-l+α+i+1)
       1cm
    +
    ∑_i,k=0^l(-1)^k+m
    (-k)_i/(α-k)_l+1(l-k)! k!(m-l+i)!Γ(m-l+α)
       =
    ∑_k=0^l(-1)^k+m-l-1211,-km-l+α+11/(α+k-l)_l+1(l-k)! k!(m-l-1)!Γ(m-l+α+1)
       1cm
    +
    ∑_k=0^l(-1)^k+m211,-km-l+11/(α-k)_l+1(l-k)! k!(m-l)!Γ(m-l+α) .
    
    Both hypergeometric series can be evaluated by means of Gauss'
    Theorem. Thus, we obtain
    R(α,m)   =∑_k=0^l(-1)^k+m+l+1/(α+k-l)_l+1(l-k)! k!(m-l-1)!Γ(m-l+α)(m-l+k+α)
       1cm
    +
    ∑_k=0^l(-1)^k+m/(α-k)_l+1(l-k)! k!(m-l-1)!Γ(m-l+α)(m-l+k).

Now we reverse the order of summation in the second sum, i.e.,
we replace k by l-k there. Then both sums can be conveniently put
together to yield

    R(α,m)   =
    ∑_k=0^l(-1)^k+l+m(α-l+2k)
    /(α-l+k)_l+1(l-k)!
      k!(m-l-1)!Γ(m-l+α)(m-l+k+α)(m-k)
       =(-1)^l+m54α-l,α/2-l/2+1,α-l+m,-m,-lα/2-l/2,1-m,α-l+m+1,α+11.
    /(α-l+1)_ll!
      (m-l-1)!Γ(m-l+α)(m-l+α)(m)

The hypergeometric function can be evaluated by means of the classical
terminating very-well-poised _5F_4-summation (cf. 
<cit.>). 
After some simplification one arrives at the desired expression
(-1)^m/(m!Γ(m+α+1)).

If m<l, then the first sum in (<ref>) does not
contribute anything because of the term (m-l-1)! in the denominator.
In the second sum, the summation over i may be started at i=l-m,
which can be evaluated by means of the binomial theorem. The result
is zero except if k=l. Again, in the end one obtains 
(-1)^m/(m!Γ(m+α+1)).
























 
We now combine Lemmata <ref> and <ref>. Let the linear combination
f=∑_ℓ=1^d+p/2 c_ℓ𝒥^d,ω_ℓ_m 
be an eigenfunction of T^d,p_m with eigenvalue λ and let ω_ℓ be as in (<ref>). We obtain for any ℓ = 1,…,d+p/2 that

    ω_ℓ^d+p    = (-1)^d+p/2 v 
        = λ^-1 2^d+p-2(d+2m-2)(-1)^d+p/2-1(2-d2)_d+p/2-1 (d+p2-1)! 
        = λ^-1 2^d+p-2(d+2m-2)(-p2)_d+p/2-1 (d+p2-1)!,

and thus

    T^d,p_m f - λ f =  (-p/2)_m/(d/2 -1)_m F(r)^⊤ A(ω)  c,

where A(ω) is as in (<ref>), c=(c_ℓ)_ℓ=1,…,d+p/2, and 

    F(r)  := ( 2^i-1(d+p2-i+1)_i-1 |ω|^-i r^m21i-d+p2,m-p2 d2 + mr^2)_i=1^d+p/2.

For i=1,…,d+p2, the hypergeometric functions in F are polynomials of exact degree d+p-2i and thus linearly independent. Hence, for c  0, the right-hand side of (<ref>)
  vanishes if and only if A(ω) is singular and c is in its nullspace. 



























§ PROOFS FOR SECTION <REF>



We shall derive the coefficients a_m(p,(3)) from the family of spherical coefficients a_m(p,^1). The half-angle identity sin(t/2)=√(1-cos(t)/2), for t∈[0,π], implies

    2^-p/2x-y^p = √(1-⟨ x,y⟩)^p = 2^p/2sin(s2)^p,     x,y∈^d-1,

where s=arccos(⟨ x,y⟩). For d=2, the addition theorem yields

    Y^m_1(x)Y^m_l(y)+Y^m_2(x)Y^m_2(y)=2𝒯_m(⟨ x,y⟩),     x,y∈^1,   m=1,2,…,

where 𝒯_m are the Chebyshev polynomials of the first kind, i.e., 𝒯_m(cos(s))=cos(ms), for s∈[0,π]. The relation (<ref>) for d=2 with (<ref>) and (<ref>) leads to 

    2^p/2sin(s2)^p = a_0(p,^1) + 2 ∑_m=1^∞ a_m(p,^1) 𝒯_m(cos(s)),     s∈[0,π],   p>0.
 

We now switch to (3). For x,y∈(3), the relation 

    x-y_=2^1/2√(3-(x^⊤ y)),
 
the choice s=arccos((x^⊤ y)-1/2), and (<ref>) imply

    2^-px-y_^p  = 2^-p/2√(3-(x^⊤ y))^p    = 2^p/2sin(s2)^p 
         = a_0(p,^1) + 2 ∑_m=1^∞ a_m(p,^1) 𝒯_m(cos(s))
        = a_0(p,^1) + 2∑_m=1^∞ a_m(p,^1) 𝒯_2m(cos(s2)).

The identity 𝒯_m = 1/2(𝒞^1_m-𝒞^1_m-2) with 𝒯_0 = 𝒞^1_0 and 𝒞^0_-1 = 0, for m ∈ℕ, leads to

    2^-px-y_^p  =∑_m=0^∞ (a_m(p,^1) - a_m+1(p,^1)) 𝒞^1_2m(cos(s2)).

By calculating the differences 2^p/2(a_m(p,^1) - a_m+1(p,^1)) in (<ref>) and applying the addition theorem of the Wigner 
𝒟-functions,

    ∑_k,l=-m^m 𝒟_k,l^m(x)𝒟_k,l^m(y) =
      (2m+1)𝒞^1_2m(cos(s2)),

we derive (<ref>). Analytic continuation and well-known relations for the gamma function cover the remaining values of p > -3.

Standard calculations yield Γ(m-p/2)/Γ(m+p/2+2)=m^-2-p(1+o(1)), which concludes the proof.




§ PROOFS FOR SECTION <REF>




 §.§ Proofs for Section <ref>



Proof of (<ref>): According to <cit.>, Q_λ is explicitly given in terms of Legendre polynomials by 

    Q_λ(x,y) =   c_λ/2(𝒞^1/2_λ_1+λ_2(ξ_+)𝒞^1/2_λ_1-λ_2(ξ_-) + 𝒞^1/2_λ_1+λ_2(ξ_-)𝒞^1/2_λ_1-λ_2(ξ_+)),

where c_λ:=(ℋ^4_2λ)=(2-δ_0,λ_2)((2λ_1+1)^2-4λ^2_2) and θ_1, θ_2 denote the principal angles between x and y and 

    ξ_+ = cos(θ_1+θ_2),   ξ_- = cos(θ_1-θ_2).

In order to write the integral (<ref>) in terms of the variables ξ_±, we first observe

    2^-p/2x-y^p_  = (2- (xy))^p/2    = (2-cos(θ_1)^2-cos(θ_2)^2)^p/2
        = (1-12cos(2θ_1)-12cos(2θ_2))^p/2
        = (1-ξ_+ξ_-)^p/2.

We set -q:=p/2 as well as m:=λ_1+λ_2 and n:=λ_1-λ_2. According to <cit.>, the measure μ__2,4⊗μ__2,4 in (<ref>) turns into dξ_+dξ_- for the variables ξ_± on |ξ_+| ≤ξ_-≤ 1, so that we obtain

    a_λ(p,_2,4)     = ∫_0^1∫_-ξ_-^ξ_- (1-ξ_+ξ_-)^-q1/2(𝒞^1/2_n(ξ_+)𝒞^1/2_m(ξ_-)+𝒞^1/2_n(ξ_-)C^1/2_m(ξ_+)) dξ_+dξ_-.
      Symmetry arguments yield
    a_λ(p,_2,4)   
         =  1/4∫_-1^1∫_-1^1 (1-xy)^-q1/2(𝒞^1/2_n(x)𝒞^1/2_m(y)+𝒞^1/2_n(y)𝒞^1/2_m(x))  dx  d y
       =  1/4∫_-1^1∫_-1^1 (1-xy)^-q𝒞^1/2_n(x)𝒞^1/2_m(y)  dx  dy.
        The series expansion of (1-xy)^-q converges absolutely for q<2, and the Legendre polynomials 𝒞^1/2_m are orthogonal to the monomials x^k for k<m or m ≢k  2. Therefore, the orthogonality relations yield 
     a_λ(p,_2,4)      = 1/4∑_k=0^∞(q)_k/k!∫_-1^1 𝒞^1/2_m(x)x^k dx ∫_-1^1  𝒞^1/2_n(y)y^k dy
       = 1/4∑_k=0^∞(q)_m+2k/(m+2k)!∫_-1^1 𝒞^1/2_m(x)x^2k+m dx ∫_-1^1 𝒞^1/2_n(y)y^2k+m dy
       =1/4∑_k=0^∞(q)_m+2k/(m+2k)! c_k^m c_m-n2+k^n,        where    c_k^m := ∫_-1^1 x^m+2k𝒞^1/2_m(x) dx.

The use of the generating function and further calculations lead to 

    c_k^m  
      = m! (m+1)_2k/2^m-1(32)_m 4^k  k!  (m+32)_k.

Application of (<ref>) and (q)_m+2k=(q)_m(m+q)_2k yields

    a_λ(p,_2,4)  
      =  1/4 (q)_m∑_k=0^∞(m+q)_2k/(m+2k)!·m! (m+1)_2k/2^m-1 (32)_m  4^k  k!  (m+32)_k
    ·n! (n+1)_m-n+2k/2^n-1(32)_n  4^m-n/2+k  (m-n2+k)!  (n+32)_m-n2+k.

By reordering and making use of n! (n+1)_m-n+2k = (m+2k)!, which
cancels the identical term in the denominator, 
we are led to

    a_λ(p,_2,4)  
         =   (q)_m  m!/4^m (32)_m∑_k=0^∞1/ (m+32)_k·(m+1)_2k/4^k·(m+q)_2k/4^k
       ·1/(m-n2+k)!·1/ (32)_n (n+32)_m-n2+k·1/k!
        =   (q)_m  m!/4^m (32)_m (m-n2)! (32)_n (n+32)_m-n2∑_k=0^∞ 
       (m+12)_k (m+22)_k (m+q2)_k (m+q+12)_k/ (m+32)_k (m-n+22)_k (m+n+32)_k·1/k!,

where the equality in the last line is due to the identities

    (m+1)_2k/4^k   =(m+12)_k (m+22)_k,    (m+q)_2k/ 4^k    = (m+q2)_k (m+q+12)_k  ,
    
       (m-n2+k)!   =(m-n2)! (m-n+22)_k,   
       (n+32)_m-n2+k   =(m+n+32)_k (n+32)_m-n2.

The relation (32)_n (n+32)_m-n2=(32)_m+n2 and 
our choices m=|λ|, n=λ_1-λ_2  yield (m-n2)!(32)_n (n+32)_m-n2=(32)_λ_1λ_2!=(32)_λ, so that q=-p/2 and the definition of 
the _4F_3 hypergeometric series conclude the proof of (<ref>).

Proof of (<ref>): 
The proof of the decay property for a_(p,_2,4) requires
some preparation and auxiliary results. For p∉2, direct calculations yield

    a_(p,_2,4)
    =2^-p/2-3/(-p/2)
    ∑_k=0^∞Γ(||/2+1/2+k) Γ(||/2+1+k) Γ(||/2+1/2-p/4+k) Γ(||/2-p/4+k)
    /Γ(||+3/2+k) Γ(_1+3/2+k) Γ(_2+1+k) Γ(k+1).

In the following, we treat the case 
_1=an and _2= (1-a)n
for n∈ and n→∞, with a fixed a∈[1/2,1], so that ||=n. Summarizing the proof of (<ref>), we shall first verify that the summand in (<ref>),

    S(n,k):=Γ(n/2+1/2+k) Γ(n/2+1+k) Γ(n/2+1/2-p/4+k) Γ(n/2-p/4+k)
    /Γ(n+3/2+k) Γ(an+3/2+k) Γ((1-a)n+1+k) Γ(k+1),

as a sequence in k, is unimodal, i.e., it first increases until it has reached its maximum and then 
decreases, see Lemma <ref>. Second, approximation of S(n,k) for k≥  n^2 with the help of Stirling's formula, where >0
but fixed, leads to an asymptotic formula for 
∑_k≥ n^2S(n,k); see Lemma <ref>. 
Third, we let →0 to obtain the asymptotic behavior of the full sum
∑_k≥0S(n,k) for n→∞; see (<ref>).
If the result is substituted in (<ref>), the claimed decay in
(<ref>) follows immediately upon observing ^2=(2a^2-2a+1)n^2.

To start with, we may consider S(n,k) as a function of real k. 
 
Given n∈, let k_0=k_0(n)∈ (0,∞) be such that ∂/∂
kS(n,k_0)=0. Then 

    k_0=k_0(n)=2a^2-2a+1/p+6n^2+O(n),
      for n→∞.


 
If n is large enough, S(n,k) has a unique — local
and global — maximum for k∈[0,∞).

 
Let >0 be fixed and let A:=2a^2-2a+1/p+6. Then

    ∑_k≥ n^2S(n,k)=A^-p+4/2
    n^-(p+4)(1+O(1n))
    ∫_-1+/A^∞
    (1+x)^-p+6/2exp(-p+62(1+x)) dx,
      for n→∞.


We postpone the proofs of Lemmata <ref>, <ref>, and <ref>, and discuss their consequences first. 
If we let →0 and substitute t=p+6/2(1+x) in the above integral, then the
integral definition of the gamma function yields 

    ∫_-1^∞(1+x)^-p+6/2exp(
    -p+6/2(1+x)) dx
    =(p+6/2)^-p+6/2 (p+4/2).

Lemma <ref> and (<ref>) provide the asymptotic
lower bound on ∑_k≥ 0S(n,k) of the following two-sided claim:

    ∑_k≥ 0S(n,k)=
    (2/2a^2-2a+1)^p+4/2(p+4/2)n^-p-4(1+o(1)),
      for n→∞.

To verify the asymptotic upper bound on ∑_k≥ 0S(n,k) in (<ref>), we observe 
    ∑_k≥0S(n,k)≤∑_k≥ n^2S(n,k)
    + n^2 S(n, n^2) ,  <A,

which is due to Lemma <ref>, saying that S(n,k) grows until its maximum at k=k_0∼ An^2. 
By Lemma <ref> and letting →0, we obtain the upper bound in (<ref>). By taking into account the additional factor 2^-p/2-3/(-p/2) in (<ref>) and the relation ^2=(2a^2-2a+1)n^2, we observe that (<ref>) provides our claim (<ref>). 

To complete the proof of (<ref>) in Theorem <ref>, it remains to prove Lemmata <ref>, <ref>, and <ref>.

The condition 0=∂/∂
kS(n,k_0) implies 0=∂/∂ k(log S(n,k_0)). By using the digamma function ψ(z), the logarithmic derivative of (<ref>)
can be written as

    ∂/∂ k(log S(n,k)) = ψ(n/2+1/2+k)
    +ψ(n/2+1+k)
    +ψ(n/2+1/2-p/4+k)
    +ψ(n/2-p/4+k)
    
     -ψ(n+3/2+k)
    -ψ(an+3/2+k)
    -ψ((1-a)n+1+k)
    -ψ(k+1).

For k=o(n), the above expression is certainly positive for large n,
hence nonzero. If the order of magnitude of k is at least the one of n
(in symbols, n=O(k)), then we may estimate the logarithmic derivative by 

    ∂/∂ k(log S(n,k)) =
    log(n/2+1/2+k)
    +log(n/2+1+k)
    +log(n/2+1/2-p/4+k)
    +log(n/2-p/4+k)
    
     -log(n+3/2+k)
    -log(an+3/2+k)
    -log((1-a)n+1+k)
    -log(k+1)
     +O(n^-1),


where the asymptotics ψ(x)=log x-1/2x+O(x^-2), for x→∞, 
cf. <cit.>), were used. 



































By applying the exponential function 
on both sides of ∂/∂ k(log S(n,k_0))=0, 
together with the above estimation for ∂/∂ k(log S(n,k)) we obtain 

    Q(n,k_0)=1+O(n^-1),

where 

    Q(n,k):=(n/2+1/2+k)
    (n/2+1+k)
    (n/2+1/2-p/4+k)
    (n/2-p/4+k)/(n+3/2+k)
    (an+3/2+k)
    ((1-a)n+1+k)
    (k+1).

For k∼ cn, we observe 

    lim_n→∞Q(n,k)=(c+1/2)^4/(c+1)(c+a)(c+1-a)c
    >1,

so that we may assume that k_0 is of larger asymptotic order of magnitude
than n. 

Define P(n,k) by 

    Q(n,k)=1-P(n,k)/(n+3/2+k)
    (an+3/2+k)
    ((1-a)n+1+k)
    (k+1).

For k of larger asymptotic order of magnitude than n, the asymptotics of the digamma function implies that the error term O(n^-1) in (<ref>) and (<ref>) may be replaced by O(k^-2) and hence by o(k^-1). The relations (<ref>) and (<ref>) lead to P(n,k_0)=o(k_0^3). Since direct computations yield

    P(n,k) =(p+6/2) k^3
    - (a^2 -a
       +1/2)(k^2n^2+k n^3) -1/16 n^4 +  lower order terms,

asymptotically leading terms in P(n,k) must cancel each other. We have already excluded k_0∼ cn, so that we now consider k_0∼ cn^2 for an appropriate constant c. The terms k^3 and k^2n^2 in P(n,k) must cancel each other, so that 

    1/2(p+6)c^3-1/2(2a^2-2a+1)c^2=0.

The solution c=1/p+6(2a^2-2a+1) yields the leading term in (<ref>). In order to derive the O(n) term in (<ref>), we have to perform “boot­strap”, i.e., we substitute k_0=cn^2+k_1 in (<ref>) and  apply analogous arguments to eventually conclude 
k_1=O(n).




We already saw in the previous proof that ∂/∂
klog S(n,0)>0, for sufficiently large n. Hence, ∂/∂
k S(n,0)>0 holds, so that S(n,k) does not have a local maximum in k=0. Convergence of the series (<ref>) implies S(n,k)→ 0 for integers k→∞. Thus, for sufficiently large n, S(n,k) attains a local maximum at some k_0∈(0,∞).  
Lemma <ref> implies k_0∼1/p+6(2a^2-2a+1)n^2. In order to investigate S(n,k) in a neighborhood of k_0, we compute ∂^2/∂
s^2log S(n,k_0+s), which is 

    ψ^(1)(n/2+1/2+k)
    +ψ^(1)(n/2+1+k)
    +ψ^(1)(n/2+1/2-p/4+k)
    +ψ^(1)(n/2-p/4+k)
    
    -ψ^(1)(n+3/2+k)
    -ψ^(1)(an+3/2+k)
    -ψ^(1)((1-a)n+1+k)
    -ψ^(1)(k+1),

where k=k_0+s and ψ^(1)(x) denotes the derivative of ψ(x). We claim that, for | s|=o( n^2) and sufficiently
large n, we have ∂^2/∂
s^2log S(n,k_0+s)<0. In order to establish this claim, we make use of ψ^(1)(x)=1/x+1/2x^2+O(x^-3), for x→∞, cf. <cit.>), to estimate the individual expressions in (<ref>). 
Using k=k_0+s=k_0+o(n^2)= k_0+o( k_0), we obtain

    ψ^(1)(n/2+1/2+k)
       =1/n/2+1/2+k
    +1/2(n/2+1/2+k)^2+O(k_0^-3)
       =1/k(1+n/2+1/2/k)
    +1/2 k_0^2+o(k_0^-2)
       =1/k(1-n/2+1/2/k
    +n^2/4k^2)
    +1/2 k_0^2+o(k_0^-2),

where we have used 1/1+x=1-x+x^2+O(x^3), for x→ 0. 
We treat the other terms in (<ref>) analogously. Note that k=k_0+o(n^2) also implies
k= 1/p+6(2a^2-2a+1)n^2+o(n^2), 
so that, by putting the individual estimates together, we arrive at

    ∂^2/∂ s^2log S(n,k_0+s)   =
    p/2+3/k^2
    +n^2/k^3(-a^2-(1-a)^2)+o(
    k_0^-2)
       =-p/2+3/k^2+o(k_0^-2)
    =-p/2+3/k_0^2+o(k_0^-2).

For sufficiently large n (and hence large k_0), this is
evidently negative, as claimed. Thus, S(n,k) has a strict local maximum in k_0. 

Moreover, say there are k_0, k'_0∈(0,∞), where S(n,k) has a local maximum, then Lemma <ref> implies that the magnitude of |k_0- k_0'| is of smaller order than n^2. The above considerations imply k_0=k_0', which completes the proof.



To prove Lemma <ref>, we shall approximate the summand S(n,k), given in (<ref>),
with the help of Stirling's formula (<ref>). 
It turns out that, under this approximation, the sum ∑_k≥0S(n,k) can then be interpreted as a Riemann integral.

Let k_0 be the unique location of the maximum of S(n,k). 
We recall that Lemma <ref> yields k_0=An^2+O(n). We consider log S(n,k), for k≥ n^2, and write k=k_0+s, so that s≥ -k_0+ n^2. 
Stirling's formula 

    log((z))=(z-1/2)log(z)-z+1/2log(2π)
    +O(z^-1)

leads to the estimates 

    log((n/2+1/2+k))
       =
    (n/2+k_0+s)
    log(n/2+1/2+k_0+s)
       2cm
    -(n/2+1/2+k_0+s)+1/2log(2π)
    +O(n^-2)
       =
    (n/2+k_0+s)[
    log(k_0+s)
    +
    log(1+n/2+1/2/k_0+s)]
       2cm
    -(n/2+1/2+k_0+s)+1/2log(2π)
    +O(n^-2),
    where we have factored out k_0+s. By applying log(1+x)=x-x^2/2+O(x^3), we obtainlog((n/2+1/2+k))   =
    (n/2+k_0+s)[
    log(k_0)
    +
    log(1+s/k_0)
    +n/2+1/2/k_0+s
    -1/2(n/2+1/2/k_0+s)^2]
       2cm
    -(n/2+1/2+k_0+s)+1/2log(2π)
    +O(n^-1).
    
    Here, terms such as n^3/(k_0+s)^2 or n/k_0+s 
    are of the order O(n^ -1) and can therefore be subsumed in the error term.
    Thus, we obtainlog((n/2+1/2+k))   =
    (n/2+k_0+s)[
    log(k_0)
    +
    log(1+s/k_0)] +1/2(n/2)^2/(k_0+s)
       3.7cm
    -(k_0+s)+1/2log(2π)
    +O(n^-1).

The reasoning for the O( . )-terms are based on our
restriction to k_0+s≥ n^2. However, the constants in these error
terms do contain . 

For the other gamma functions in (<ref>), we proceed similarly. If everything is put together, then we obtain

    log(S(n,k))
       =-p+6/2[log(k_0)+log(1+s/k_0)]
    -1/2n^2(2a^2-2a+1)/(k_0+s)
    +O(n^-1)
       =-p+6/2[log(k_0)+log(1+s/k_0)]
    -p+6/2(1+s/k_0)
    +O(n^-1).

For the sum of the S(n,k), we have

    ∑_k≥ n^2S(n,k)
    =∑_s≥ -k_0+ n^2
    k_0^-p+6/2(1+s/k_0)^-p+6/2exp(
    -(p+6)/2(1+sk_0))
    (1+O(n^-1)).

By understanding, the sum over s is taken over those s, for which
k_0+s is an integer. The error of the Riemann sum approximation 


    1/k_0∑_s≥-k_0+ n^2(1+s/k_0)^-p+6/2exp(
    -(p+6)/2(1+s/k_0))
    ≈∫_-1+/A^∞
    (1+x)^-p+6/2exp(-(p+6)/2(1+x)) dx

is of the order of magnitude O(k_0^-1). This follows from the 
fact that the summand in the sum attains a unique local and global
maximum — namely at s=0 — and therefore the error is bounded 
above by k_0^-1 times the absolute variation of the summand — 
which equals twice the maximum. Substitution of all this in (<ref>) concludes the proof. 

Our proof of (<ref>) in Theorem <ref> is now complete.





 §.§ Proofs for Section <ref>

Recall that the action of (3)×(3) on 𝕊^2×𝕊^2 by left multiplication is transitive, and  (4) acts transitively on _2,4 by conjugation. For each O∈(4), there are a,b∈𝕊^3 such that O=L_a R_b, where the left and right isoclinic rotations are

    L_a :=
    [  a_1 -a_2 -a_3 -a_4;  a_2  a_1 -a_4  a_3;  a_3  a_4  a_1 -a_2;  a_4 -a_3  a_2  a_1;      ]^⊤,   
    R_b :=
    [  b_1 -b_2 -b_3 -b_4;  b_2  b_1  b_4 -b_3;  b_3 -b_4  b_1  b_2;  b_4  b_3 -b_2  b_1;      ].

For a∈𝕊^3, the Euler–Rodrigues formula yields

    S_a:=
    [ a_1^2+a_2^2-a_3^2-a_4^2        2(a_2a_3-a_1a_4)        2(a_2a_4+a_1a_3);       2( a_2a_3+a_1a_4) a_1^2-a_2^2+a_3^2-a_4^2       2( a_3a_4-a_1a_2);       2( a_2a_4-a_1a_3)        2(a_3a_4+a_1a_2) a_1^2-a_2^2-a_3^2+a_4^2;                         ]^⊤∈(3).

We are looking for 𝒫:^2×^2→_2,4 satisfying

    (L_a R_b)·𝒫(x,y)·(L_a R_b)^⊤
      =𝒫(S_a x,S_b y),     a,b∈𝕊^3,   x,y∈𝕊^2.











There are exactly two mappings ^2×^2→_2,4 satisfying (<ref>). One is 𝒫 as in (<ref>), and the other is I_4-𝒫. In particular, 𝒫 is surjective and, for all x,y,u,v∈𝕊^2,

    𝒫(u,v)=𝒫(x,y)  if and only if   (u,v) ∈{± (x,y) }.



The identity (<ref>) for the specific choice of 𝒫 is verified by expanding both sides of the equality and comparing the polynomial expressions. We omit the straightforward but lengthy computation.

Since the conjugate action of (4) on _2,4 is transitive, the identity (<ref>) also implies surjectivity.
Since left and right eigenspaces of ℒ(𝒫(x,y)) in (<ref>) and (<ref>) are uniquely determined, we deduce that (<ref>) holds.

Let us now address the uniqueness statement. For a:=(α,0,0), b:=(β,0,0) with α,β∈𝕊^1, we obtain the isoclinic rotations

    L_a = [ A_α   0;   0 A_α ],
        
    R_b = [   B_β     0;     0 B^⊤_β ],

where A_α=([  α_1 -α_2;  α_2  α_1 ]), B_β = ([  β_1 -β_2;  β_2  β_1 ])∈(2). Since S_ae_1=S_b e_1=e_1, any mapping 𝒫̃:^2×^2→_2,4 satisfying (<ref>) must obey

    L_a R_b𝒫̃(e_1,e_2)(L_a R_b)^⊤ = 𝒫̃(e_1,e_2),  α,β∈𝕊^1.

Let V denote the range of 𝒫̃(e_1,e_2), which is a two-dimensional subspace of ℝ^4. The relation (<ref>) means that L_aR_b maps V into itself, i.e., L_aR_b V = V. For all U_1,U_2∈(2), there are α,β∈𝕊^1 such that

    L_aR_b = [   A_α B_β         0;         0 A_α B_β^⊤ ] =
    [ U_1   0;   0 U_2 ],

so that V must either coincide with {e_1,e_2} or with {e_3,e_4}. Hence, 𝒫̃(e_1,e_2) must coincide with either (1,1,0,0) or (0,0,1,1). It follows from (<ref>) and (3) acting transitively on 𝕊^2 that in the first case 𝒫̃=𝒫 and in the latter 𝒫̃=I_4-𝒫.


The Grassmannian _2,4 has been parametrized in <cit.> by means of angles but the explicit identity (<ref>) that steers our present approach has not been considered there. Our parametrization is compatible with the one used in <cit.>.



The following properties of 𝒫 are useful.

The probability measure σ__2,4, induced by the Haar measure on (4), is the push-forward measure of σ_𝕊^2⊗σ_𝕊^2 under 𝒫, and we have, for all x,y,u,v∈^2,

    ⟨𝒫(x,y),𝒫(u,v)⟩_   =1+⟨ xy^⊤,uv^⊤⟩_=1+ ⟨ x,u⟩⟨ y,v⟩, 
    ℐ_-·𝒫(x,y)·ℐ_-^⊤    = 𝒫(y,x),
    
    I_4-𝒫(x,y)     = 𝒫(-x,y) = 𝒫(x,-y),
    {⟨ x,u⟩,⟨ y,v⟩}   = ±{
    cos(θ_1+θ_2),cos(θ_1-θ_2)
    },

where ℐ_-:=(-1,1,1,1), and θ_1, θ_2 denote the principal angles between 𝒫(x,y) and 𝒫(u,v).







Note that the right-hand side of (<ref>) is ±{ξ_+,ξ_-} for ξ_+,ξ_- as in (<ref>). 

The product measure σ_𝕊^2⊗σ_𝕊^2 is (3)×(3) invariant. According to (<ref>), the pushforward measure of σ_𝕊^2⊗σ_𝕊^2 on 𝕊^2 ×𝕊^2 under 𝒫 is (4) invariant, so that the uniqueness of the Haar measure implies the first claim of the lemma.

Each of the remaining claims is first proved for 𝒫 as defined in Proposition <ref> and then argued that it also holds for I_4-𝒫. The identity (<ref>) is easily observed for u,v=e_1 first and then (<ref>) yields the general case. According to

    ⟨ I_4-𝒫(x,y),I_4-𝒫(u,v)⟩_ =⟨𝒫(x,y),𝒫(u,v)⟩_,

the identity (<ref>) also holds for I_4-𝒫. The statement (<ref>) follows from expanding both sides of the equality and comparing polynomial expressions in x and y. If it holds for 𝒫, then it must also hold for I_4-𝒫. One directly calculates (<ref>). 

In order to check (<ref>), we first recognize that principal angles between I-𝒫(x,y) and I-𝒫(u,v) coincide with the ones between 𝒫(x,y) and 𝒫(u,v). By ξ_± as in (<ref>) and as at the beginning of the proof of Theorem <ref>, we observe (𝒫(x,y)𝒫(u,v))=1+ξ_+ξ_-. Hence, (<ref>) leads to
ξ_+ξ_- = ⟨ x,u⟩⟨ y,v⟩. 
Theorem <ref> implies that Q_λ in (<ref>) satisfies

    Q_λ(𝒫(x,y),𝒫(u,v)) = c_λ(𝒞^1/2_λ_1+λ_2(⟨ x,u⟩)𝒞^1/2_λ_1-λ_2(⟨ y,v⟩) + 𝒞^1/2_λ_1+λ_2(⟨ y,v ⟩)𝒞^1/2_λ_1-λ_2(⟨ x,u ⟩)).

Comparison of this identity with (<ref>) and few further calculations eventually lead to (<ref>). 














The action of (3)×(3) leads to the irreducible decomposition 
  
    L_2(𝕊^2×𝕊^2± 1)    =⊕_m+n∈2ℕspan{ Y_k^m⊗ Y_l^n : k = -m,…,m,  l=-n,…,n}.

The unit quaternions provide a group structure on ^3, so that the mapping a↦ S_a between 𝕊^3± 1 and (3) as well as (a,b)↦ L_a R_b between
(𝕊^3×𝕊^3)± 1 and (4) become group isomorphisms. Their combination induces a group isomorphism between (3)×(3) and (4)± 1. 

Condition (<ref>) requires that the respective actions on ^2×^2 and _2,4 commute with 𝒫. Hence, the induced pullback

    𝒫^*:L_2(_2,4)→  L_2(𝕊^2×𝕊^2± 1),     f↦ f(𝒫(·,·))

 is an intertwining isomorphism. In particular, 𝒫^* maps one irreducible subspace into the other.

 Thus, the irreducible decomposition of L_2(_2,4) under the action of (4) is
  
    L_2(_2,4) =⊕_m+n∈2ℕspan{ Y_k,l^m,n: k = -m,…,m,  l=-n,…,n}.

In comparison to the irreducible decomposition of L_2(_2,4) for the action of (4) in (<ref>), the irreducible components H_λ(_2,4) with respect to (4) in (<ref>) are usually larger. 

Let us denote H^m,n:=span{ Y_k,l^m,n: k = -m,…,m,  l=-n,…,n}. Property (<ref>) yields that the
irreducible subspaces of L_2(_2,4) under the action of (4) are

    V^m,n:=
    
    H^m,n ,    m=n,
    
    H^m,n⊕ H^n,m,    m≠ n,

where m+n∈2ℕ and m≥ n, i.e.,

    L_2(_2,4) = ⊕_m+n∈2ℕ
     m≥ n V^m,n=⊕_λ_1≥λ_2≥ 0 V^m_λ,n_λ,

where m_λ=λ_1+λ_2 and n_λ=λ_1-λ_2. 

By considering ℒ(P)=xy^⊤, ℒ(P)·ℒ(P)^⊤=xx^⊤, and ℒ(P)^⊤·ℒ(P)=yy^⊤, we observe that the homogeneous polynomials
x_iy_j and x_ix_j, y_iy_j, for i,j = 1,…,3, 
can be written as homogeneous polynomials in the matrix entries of P ∈_2,4 of degree 1 and 2, respectively. The monomial x^α y^β with |α|=m_λ and |β|=n_λ is composed of n_λ factors of the form x_iy_j and (m_λ-n_λ)/2 factors of the form x_ix_j. Thus, x^α y^β is a homogeneous polynomial of degree n_λ+2((m_λ-n_λ)/2) = m_λ in the matrix entries of P. We deduce

    ⊕_λ_1+λ_2≤ t V^m_λ,n_λ⊂⊕_λ_1+λ_2≤ t H_λ(_k,d)

since the right-hand side of (<ref>) are the polynomials of degree at most t in the matrix entries of P∈_2,4, cf. <cit.>. 
By applying <cit.>, 
we see that the dimension of H_λ(𝒢_2,4) is

    dim(H_λ(𝒢_2,4)) =(2-δ_m_λ,n_λ)(2m_λ+1)(2n_λ+1),

which matches (V^m_λ,n_λ). Hence, there holds equality in (<ref>). For fixed t, the dimensions in (<ref>) are pairwise different, for all λ_1+λ_2=t. An induction over t leads to H_λ(𝒢_2,4)=V^m_λ,n_λ.





 §.§ Proofs for Section <ref>

The relation (<ref>) yields

    F(x,y) = ∑_k,l=-M^M e^ ikφ_1e^ ilφ_2 B_k,l(θ_1,θ_2),
 
where B_k,l(θ_1,θ_2) are given by

    B_k,l(θ_1,θ_2)     := ∑_m_1=|k|^M ∑_m_2=|l|^M f^m_1,m_2_k,l∑_k'=-m_1^m_1∑_l'=-m_2^m_2 c^m_1_k,k' c^m_2_l,l' e^ ik'θ_1 e^ il'θ_2
        = ∑_k',l'=-M^M e^ ik'θ_1e^ il'θ_2∑_m_1=max(|k|,N-k')^M   ∑_m_2=max(|l|,M-l')^M f^m_1,m_2_k,l  c^m_1_k,k' c^m_2_l,l'.

Hence, the coefficients from (<ref>) satisfy

    b^k,l_k',l' = ∑_m_1=max(|k|,M-k')^M c^m_1_k,k'∑_m_2=max(|l|,M-l')^M f^m_1,m_2_k,l   c^m_2_l,l'.

First evaluating the inner sum for k,l,l'=-M,…,M and m_1=max(|k|,M-k'),…,M, and afterwards the outer sum enables the computation of the coefficients b^k,l_k',l', for k,l,k',l'=-M,…,M, in O(M^5) operations provided that the numbers  c^m_k,k' in (<ref>) are given. 


The complexity for evaluating the two sums in (<ref>) can be further reduced to O(M^4log^2(M)) by using a fast polynomial transform, cf. <cit.>. 
In this way, the nonequispaced fast Fourier transform on ^2×^2 
takes O(M^4log^2(M)+n|log(ϵ)|^4) elementary operations. For n∼ M^4, this implies the reduction of complexity from O(M^8) to O(M^4log^2(M)+M^4|log(ϵ)|^4) operations.


The adjoint nonequispaced fast Fourier transform is about computing 

    ∑_i=1^M b_i Y^m_k(x_i)Y^n_l(y_i),   m,n=0,…,N,   k=-m,…,m,   l=-n,…,n,
 
for given coefficients (b_i)_i=1^M⊂ℂ and locations (x_i,y_i)_i=1^M⊂^2×^2. Similar arguments as above and in Section <ref> provide a fast evaluation of (<ref>) by using the adjoint nonequispaced fast Fourier transform , cf. <cit.>.













10

Adcock:2011by
B. Adcock, On the convergence of expansions in polyharmonic eigenfunctions, J. Approx. Theory
  163 (2011), no. 11, 1638–1674.


Alexander:1975rp
R. Alexander, Generalized sums of distances, Pacific J. Math.
  56 (1975), no. 2, 297–304.

Alpay:2015mb
D. Alpay and P. Jorgensen, Spectral theory for Gaussian processes: reproducing kernels, boundaries, and L^2-wavelet generators with fractional scales, Numer. Funct. Anal. Optim. 36 (2015), no. 10, 1239–1285.

Ambrosio:2005lh
L. Ambrosio, N. Gigli, and G. Savaré, Gradient Flows in Metric Spaces and in the Space of Probability Measures, Birkhäuser Verlag, 2005.

Aronszajn:1983mb
N. Aronszajn, T. M. Creese, and L. J. Lipkin, Polyharmonic Functions, Clarendon Press, 1983.

Bachoc:2006aa
C. Bachoc, Linear programming bounds for codes in Grassmannian spaces,
  IEEE Trans. Inf. Th. 52 (2006), no. 5, 2111–2125.

Bachoc:2004fk
C. Bachoc, E. Bannai, and R. Coulangeon, Codes and designs in
  Grassmannian spaces, Discrete Math. 277 (2004), 15–28.

Bachoc:2002aa
C. Bachoc, R. Coulangeon, and G. Nebe, Designs in Grassmannian spaces
  and lattices, J. Algebraic Combin. 16 (2002), 5–19.

Baouendi:1974pr
M. Baouendi, C. Goulaouic, and L. Lipkin, On the operator Δ r^2+μ(∂/∂ r)r+λ, J. Differ. Equations 15 (1974), 499–509.



Baxter:2001pb
B. J. C. Baxter and S. Hubbert, Radial basis functions for the sphere,
in:  Recent Progress in Multivariate Approximation (Witten-Bommerholz, 2000),
Internat. Ser. Numer. Math., vol. 137, Birkhäuser, Basel, 2001, 
pp. 33–47.

Brandolini:2014oz
L. Brandolini, C. Choirat, L. Colzani, G. Gigante, R. Seri, and G. Travaglini,
  Quadrature rules and distribution of points on manifolds, Annali della
  Scuola Normale Superiore di Pisa – Classe di Scienze XIII (2014),
  no. 4, 889–923.

Brauchart:2013il
J. S. Brauchart and J. Dick, A characterization of Sobolev spaces on
  the sphere and an extension of Stolarsky's invariance principle to
  arbitrary smoothness, Constr. Approx. 38 (2013), no. 3, 397–445.

Brauchart:2013pt
, A simple proof of Stolarsky's invariance principle, Proc.
  Amer. Math. Soc. 141 (2013), no. 6, 2085–2096.

Brauchart:fk
J. S. Brauchart, E. B. Saff, I. H. Sloan, and R. S. Womersley, QMC
  designs: Optimal order quasi Monte Carlo integration schemes on the
  sphere, Math. Comp. 83 (2014), 2821–2851.

Breger:2016vn
A. Breger, M. Ehler, and M. Gräf, Quasi Monte Carlo integration
  and kernel-based function approximation on Grassmannians,
in: Frames and Other
  Bases in Abstract and Function Spaces: Novel Methods in Harmonic Analysis,
  vol. 1, 
Birkhäuser/Springer, 2017, pp. 333–353.

Breger:2016rc
A. Breger, M. Ehler, and M. Gräf, Points on manifolds with
  asymptotically optimal covering radius, J. Complexity 48 (2018), 1–14.


Chauffert:2017eb
N. Chauffert, P. Ciuciu, J. Kahn, and P. Weiss, A projection method on measures sets, Constr. Approx. 45 (2017), no. 1, 83–111.







Cohl:2013ao
H. S. Cohl, On a generalization of the generating function for Gegenbauer polynomials, Integral Transforms and Special Functions 24 (2013), no. 10, 807–816.


Conway:1996aa
J. H. Conway, R. H. Hardin, and N. J. A. Sloane, Packing lines, planes, etc.: packings in Grassmannian space, Experimental Math. 5
  (1996), 139–159.


Davis:1999dn
A. W. Davis, Spherical functions on the Grassmann manifold and
  generalized Jacobi polynomials – Part 2, Lin. Alg. Appl. 289
  (1999), no. 1–3, 95–119.

Dick:2014bj
J. Dick and F. Pillichshammer, Discrepancy theory and quasi-Monte Carlo integration, in: A Panorama of Discrepancy Theory, W. Chen, A. Srivastav, G. Travaglini (eds.). Lecture Notes in Mathematics, vol. 2107. Springer, Cham, pp. 539–619.

ErdeAA
A. Erdélyi, V. Magnus, F. Oberhettinger and F. Tricomi, 
Higher Transcendental Functions, vol. 1, McGraw–Hill, New York,
1953.






Ehler:2014zl
M. Ehler and M. Gräf, Reproducing kernels for the irreducible
  components of polynomial spaces on unions of Grassmannians, 
  Constr. Approx. 49 (2018), no. 1, 29–58.

Etayo:2016qq
U. Etayo, J. Marzo, and J. Ortega-Cerdà, Asymptotically optimal
  designs on compact algebraic manifolds, 
Monatsh. Math. 186 (2018), 235–248.
  







Fulton:1991fk
W. Fulton and J. Harris, Representation Theory: a First Course, Springer, 1991.


Fuselier:2012jt
E. Fuselier and G. B. Wright, Scattered data interpolation on embedded
  submanifolds with restricted positive definite kernels: Sobolev error
  estimates, SIAM J. Numer. Anal. 50 (2012), no. 3, 1753–1776.

Gilbarg:2001vl
D. Gilbarg and N. S. Trudinger, Elliptic Partial Differential Equations of Second Order, Springer, 2001. 


Gneiting:1999ay
T. Gneiting, Radial positive definite functions generated by Euclid's
  hat, J. Multivariate Analysis 69 (1999), 88–119.

Gnewuch:2012jy
M. Gnewuch, Weighted geometric discrepancies and numerical integration on reproducing kernel Hilbert spaces, J. Complexity 28 (2012), 2–17.

Gournay:2019kt
F. de Gournay, J. Kahn, L. Lebrat, and P. Weiss, Optimal transport approximation of 2-dimensional measures, SIAM J. Imaging Sciences
12 (2019), 762–787.

Graf:2013zl
M. Gräf, Efficient Algorithms for the Computation of Optimal Quadrature
  Points on Riemannian Manifolds, Ph.D. thesis, TU Chemnitz,
Universitätsverlag Chemnitz, 2013.

Graf:2009ye
M. Gräf and D. Potts, Sampling sets and quadrature formulae on the
  rotation group, Numer. Funct. Anal. Optim. 30 (2009), 665–688.

Graf:2011lp
, On the computation of spherical designs by a new optimization
  approach based on fast spherical Fourier transforms, Numer. Math.
  119 (2011), 699–724.

Graf:2013fk
M. Gräf, M. Potts, and G. Steidl, Quadrature errors, discrepancies and
  their relations to halftoning on the torus and the sphere, SIAM J. Sci.
  Comput. 34 (2012), A2760–A2791.

James:1974aa
A. T. James and A. G. Constantine, Generalized Jacobi polynomials as
  spherical functions of the Grassmann manifold, Proc. London Math. Soc.
  29 (1974), no. 3, 174–192.
  
John:1950ek
F. John, The fundamental solution of a linear elliptic differential equations with analytic coefficients, Comm. Pure Appl. Math.
  3 (1950), no. 3, 273–304.
  

Kalmenov:2011to
T. Sh. Kal'menov and D. Suragan, A boundary condition and spectral problems for the Newton potential, Operator Theory: Advances and Applications 216 (2011), 187–210.
 
 
Kalmenov:2012fq
T. Sh. Kal'menov and D. Suragan, Boundary conditions for the volume potential for the polyharmonic equation, J. Differ. Equations 48 (2012), no. 4, 604–608.


Keiner:2009cv
J. Keiner and S. Kunis and D. Potts, Using NFFT3 – a software library for various nonequispaced fast Fourier transforms, ACM Trans. Math. Software
  36 (2009), 1–30.

Kuipers:1974la
L. Kuipers and H. Niederreiter, Uniform Distribution of Sequences,
  Wiley, 1974.
  
Kunis:2003pt
S. Kunis and D. Potts, Fast spherical Fourier algorithms, J. Comput. Appl. Math.
  161 (2003), 75–98.


Matousek:2010kb
J. Matoušek, Geometric Discrepancy, Algorithms and Combinatorics,
  vol. 18, Springer, 2010.

Nowak:2010rr
E. Novak and 
H. Woźniakowski, Tractability of Multivariate
  Problems. Volume II, EMS Tracts in Mathematics, vol. 12, EMS
  Publishing House, Zürich, 2010.

Potts:2009gb
D. Potts, J. Prestin, and A. Vollrath, A fast algorithm for nonequispaced
  Fourier transforms on the rotation group, Numer. Algorithms 52
  (2009), 355–384.

Roy:2010fk
A. Roy, Bounds for codes and designs in complex subspaces, 
J. Algebraic Combin. 31 (2010), no. 1, 1–32.


Sk20
M. M. Skriganov, Stolarsky's invariance principle for projective spaces, to appear in J. Complexity, 2020.

SlatAC 
L. J. Slater, 
Generalized Hypergeometric Functions,
Cambridge University Press, Cambridge, 1966.

Varshalovich:1988vh
D. Varshalovich, A. Moskalev, and V. Khersonskii, Quantum Theory of Angular Momentum, World Scientific, Singapore, 1988.

Wendland:2004wd
H. Wendland, Scattered Data Approximation, Cambridge Monographs on
  Applied and Computational Mathematics (17), Cambridge University Press, 2004.




