## Algorithm xxx: FaVeST - Fast Vector Spherical Harmonic Transforms

QUOC T. LE GIA, The University of New South Wales, Australia MING LI*, Zhejiang Normal University, China; La Trobe University, Australia

YU GUANG WANG*, The University of New South Wales, Australia

Vector spherical harmonics on the unit sphere of R 3 have broad applications in geophysics, quantum mechanics and astrophysics. In the representation of a tangent vector field, one needs to evaluate the expansion and the Fourier coefficients of vector spherical harmonics. In this paper, we develop fast algorithms (FaVeST) for vector spherical harmonic transforms on these evaluations. The forward FaVeST evaluates the Fourier coefficients and has a computational cost proportional to 𝑁 log √ 𝑁 for 𝑁 number of evaluation points. The adjoint FaVeST which evaluates a linear combination of vector spherical harmonics with a degree up to √ 𝑀 for 𝑀 evaluation points has cost proportional to 𝑀 log √ 𝑀 . Numerical examples of simulated tangent fields illustrate the accuracy, efficiency and stability of FaVeST.

CCS Concepts: · Mathematics of computing → Mathematical analysis ; Numerical analysis ; Computations of transforms;

Additional Key Words and Phrases: Vector spherical harmonics, tangent vector fields, FFT

## 1 INTRODUCTION

Vector spherical harmonics on the unit sphere S 2 in R 3 are widely used in many areas such as astrophysics [1-3], quantum mechanics [32, 38, 69], geophysics and geomagnetics [5, 15, 23, 41, 42], 3D fluid mechanics [63, 64], global atmospherical modelling [13, 19, 33] and climate change modelling [56-58]. For example, in constructing the numerical solution to the Navier-Stokes equations on the unit sphere [17], divergence-free vector spherical harmonics are used. In simulating scattering waves by single or multiple spherical scatterers [47, 65] modelled by the 3-dimensional Helmholtz equation, both divergence-free and curl-free vector spherical harmonics are used. In these problems, the solutions which are vector fields are represented by an expansion of vector spherical harmonics. One then needs to evaluate the coefficients of vector spherical harmonics for a vector field and also a linear combination of vector spherical harmonics. They can be evaluated by the (discrete) forward vector spherical harmonic transform (FwdVSHT) and adjoint vector spherical harmonic transform (AdjVSHT) respectively. In this paper, we develop fast algorithms for the forward and adjoint vector spherical harmonic transforms for tangent (vector) fields on S 2 and their software implementation.

Let {( y ℓ,𝑚 , z ℓ,𝑚 ) : ℓ = 1 2 , , . . . , 𝑚 = -ℓ, . . . , ℓ } be a set of pairs of the (complex-valued) divergence-free and curl-free vector spherical harmonics on S 2 . The coefficients for divergence-free and curl-free vector spherical harmonics are given by, for ℓ = 1 2 , , . . . , 𝑚 = -ℓ, . . . , ℓ ,

$$\widehat { T } _ { \ell, m } = \int _ { \mathbb { S } ^ { 2 } } T ( { \mathbf x } ) { \mathbf y } _ { \ell, m } ^ { * } ( { \mathbf x } ) { \mathrm d } \sigma ( { \mathbf x } ) \quad \widetilde { T } _ { \ell, m } = \int _ { \mathbb { S } ^ { 2 } } T ( { \mathbf x } ) { \mathbf z } _ { \ell, m } ^ { * } ( { \mathbf x } ) { \mathrm d } \sigma ( { \mathbf x } ),$$

where the 𝑉 ∗ is the complex conjugate transpose of the tangent field 𝑉 . The FwdVSHT for a spherical tangent field 𝑇 : S 2 → C 3 evaluates these coefficients by approximating the integrals of the coefficients with a quadrature rule which

* Corresponding authors.

Authors' addresses: Quoc T. Le Gia, The University of New South Wales, Sydney, NSW, Australia, qlegia@unsw.edu.au; Ming Li*, Zhejiang Normal University, Jinhua, Zhejiang, China; La Trobe University, Melbourne, VIC, Australia, mingli@zjnu.edu.cn; Yu Guang Wang*, The University of New South Wales, Sydney, NSW, Australia, yuguang.wang@unsw.edu.au.

is a set of 𝑁 pairs of weights 𝑤 𝑖 and points x 𝑖 on S 2 :

$$\widehat { T } _ { \ell, m } \approx \sum _ { l = 1 } ^ { N } w _ { i } T ( x _ { i } ) y _ { \ell, m } ^ { * } ( x _ { i } ) \quad \widetilde { T } _ { \ell, m } \approx \sum _ { l = 1 } ^ { N } w _ { i } T ( x _ { i } ) z _ { \ell, m } ^ { * } ( x _ { i } ).$$

The AdjVSHT evaluates the expansion of y ℓ,𝑚 , z ℓ,𝑚 with two complex sequences 𝑎 ℓ,𝑚 , 𝑏 ℓ,𝑚 , ℓ = 1 2 , , . . . , 𝑚 = -ℓ, . . . , ℓ as coefficients for a set of spherical points { x 𝑖 } 𝑀 𝑖 = 1 , 𝑀 ≥ 1:

$$\sum _ { \ell = 1 } ^ { \infty } \sum _ { m = - \ell } ^ { \ell } \left ( a _ { \ell, m } y _ { \ell, m } ( x _ { i } ) + b _ { \ell, m } z _ { \ell, m } ( x _ { i } ) \right ), \ \ i = 1, \dots, M.$$

To directly compute the summation of (1) for the coefficients with degree up to 𝐿 for 𝐿 ≥ 1, the computational cost for FwdVSHT is O GLYPH&lt;0&gt; 𝑁𝐿 2 GLYPH&lt;1&gt; . Here the quadrature rule should be properly chosen to minimize the approximation error. The optimal-order number of nodes for this purpose is 𝑁 = O GLYPH&lt;0&gt; 𝐿 2 GLYPH&lt;1&gt; . (See Section 4.3.) Thus, the cost of direct computation for forward vector spherical harmonic transform is O GLYPH&lt;0&gt; 𝑁 2 GLYPH&lt;1&gt; . On the other hand, to directly compute the expansion (2) with truncation degree 𝐿 for ℓ incurs O GLYPH&lt;0&gt; 𝑀𝐿 2 GLYPH&lt;1&gt; computational steps. Then, to evaluate 𝑀 = O GLYPH&lt;0&gt; 𝐿 2 GLYPH&lt;1&gt; points, the computational complexity for direct evaluation of AdjVSHT is O GLYPH&lt;0&gt; 𝑀 2 GLYPH&lt;1&gt; .

In this paper, we develop fast computational strategies for FwdVSHT and AdjVSHT by explicit representations for the divergence-free and curl-free vector spherical harmonics in terms of scalar spherical harmonics. Here the close representation formula exploits Clebsch-Gordan coefficients in quantum mechanics. By this way, the fast scalar spherical harmonic transforms can be applied to speed up computation. The resulting algorithm reduces the computational cost to O GLYPH&lt;16&gt; 𝑁 log √ 𝑁 GLYPH&lt;17&gt; and O GLYPH&lt;16&gt; 𝑀 log √ 𝑀 GLYPH&lt;17&gt; , both of which are nearly linear. We thus call the algorithms the Fa st Ve ctor S pherical Harmonic T ransform or FaVeST . A software package in Matlab is provided for implementing FaVeST . We then validate the FaVeST algorithm by the numerical examples of simulated tangent fields on the sphere. The algorithm, accompanied by its software implementation fills the blank of fast algorithms for vector spherical harmonic transforms.

The rest of this paper is organized as follows. In Section 2, we review the fast Fourier transforms for scalar spherical harmonics on S 2 and related works on the computation of vector spherical harmonic transforms. In Section 3, we introduce definitions and notation on scalar and vector spherical harmonics. In Sections 4.1 and 4.2, we give the representations for FwdVSHT and AdjVSHT in terms of forward and adjoint scalar spherical harmonic transforms. From these representations, in Section 4.3, we describe the fast algorithms ( FaVeST ) for the evaluation of FwdVSHT and AdjVSHT, and show that the computational complexity of the proposed FaVeST is nearly linear. We also estimate the approximation error for forward and adjoint FaVeST s. Section 4.4 describes the software package in Matlab developed for FaVeST . In Section 5, we give numerical examples of simulated spherical tangent fields to test FaVeST .

## 2 RELATED WORKS

Fast Fourier Transform (FFT) for R 𝑑 is one of the most influential algorithms in science and engineering [8, 48, 49, 59]. On the sphere, fast transforms for scalar spherical harmonics have been extensively studied by many researchers [4, 11, 20, 24, 25, 35, 36, 36, 43, 46, 51, 53-55, 60, 61]. In particular, Keiner et al. provide the software library NFFT [36] which implements the fast forward and adjoint FFT algorithms for scalar spherical harmonics based on the non-equispaced FFT [35, 37, 39]. Their package is easy to use in Matlab environment and has been applied in many areas. Suda and Takami in [55] propose a fast scalar spherical harmonic transform algorithm with computational complexity O( 𝑁 log √ 𝑁 ) based on the divide-and-conquer approach with split Legendre functions (where 𝑁 is the number of nodes in the discretization for integral on S 2 ), and the algorithm is used to solve the shallow water equation

[54]. Rokhlin and Tygert [51] develop the fast algorithms for scalar spherical harmonic expansion with computational time proportional to 𝑁 log √ 𝑁 log ( 1 / ) 𝜖 for a given precision 𝜖 &gt; 0. Later, Tygert [60, 61] improves their algorithm to achieve computational cost proportional to 𝑁 log √ 𝑁 at any given precision. Reinecke and Seljebotn and Gorski et al. [20, 50] develop the FFTs for spherical harmonics to evaluate Hierarchical Equal Area isoLatitude Pixelation (HEALPix) points [20] in various programming languages. Their algorithm can work on millions of evaluation points with spherical harmonic degree up to 6 143. ,

In contrast, fast transforms for vector spherical harmonics receive less attention. To the best of our knowledge, there are no existing fast algorithms for the forward and adjoint vector spherical harmonic transforms. Ganesh et al. [17] use FFTs to speed up their algorithms for solving Navier-Stokes PDEs on the unit sphere. However, their method is based on the idea that applies conventional FFTs to evaluate complex azimuthal exponential terms involved in the formulation of vector spherical harmonics. As fast Legendre transforms are not implemented, their method is not a fast transform. Wang et al. [65] evaluate vector spherical harmonic expansions via spectral element grids, which is neither fast computation.

## 3 VECTOR SPHERICAL HARMONICS

In this section, we present definitions and properties about spherical tangent fields, scalar and vector spherical harmonics and Clebsch-Gordan coefficients, see, e.g. [9, 12], which we will use in the representation of FwdVSHT and AdjVSHT in the next section. A tangent (vector) field 𝑇 is a mapping from S 2 to C 3 satisfying the normal component ( 𝑇 · x x ) of 𝑇 is zero, here 𝑇 · x : = ˝ 3 𝑖 = 1 𝑇 ( 𝑖 ) 𝑥 ( 𝑖 ) is the inner product of C 3 for column vectors 𝑇 : = GLYPH&lt;0&gt; 𝑇 ( 1 ) , 𝑇 ( 2 ) , 𝑇 ( 3 ) GLYPH&lt;1&gt; ′ and x : = ( 𝑥 ( 1 ) , 𝑥 ( 2 ) , 𝑥 ( 3 ) ) ′ , and the ′ denotes the transpose of a vector (or matrix). Let L 2 ( S 2 ) be 𝐿 2 space of tangent fields on the sphere S 2 with inner product

$$\langle T, V \rangle = \int _ { \mathbb { S } ^ { 2 } } T ^ { * } ( \mathbf x ) V ( \mathbf x ) d \sigma ( \mathbf x )$$

and 𝐿 2 norm ∥ 𝑇 ∥ 2 = √︁ ⟨ 𝑇,𝑇 ⟩ , where 𝑇 ∗ ( x ) is the complex conjugate transpose of 𝑇 ( x ) . Using spherical coordinates, the scalar spherical harmonics can be explicitly written as, for ℓ = 0 1 , , . . . ,

$$Y _ { \ell, m } ( { \mathbf x } ) \coloneqq Y _ { \ell, m } ( \theta, \varphi ) \coloneqq \sqrt { \frac { 2 \ell + 1 } { 4 \pi } \frac { ( \ell - m )! } { ( \ell + m )! } } P _ { \ell } ^ { ( m ) } ( \cos \theta ) \, e ^ { i m \varphi }, \ \ m = 0, 1, \dots, \ell,$$

$$Y _ { \ell, m } ( { \mathbf x } ) \coloneqq ( - 1 ) ^ { m } Y _ { \ell, - m } ( { \mathbf x } ), \ \ m = - \ell, \dots, - 1.$$

We would suppress the variable x in 𝑌 ℓ,𝑚 ( x ) if no confusion arises.

In the following, we introduce vector spherical harmonics, see, e.g. [12, 62]. For each ℓ = 1 2 , , . . . , 𝑚 = -ℓ, . . . , ℓ , and integers 𝑗 1 , 𝑗 2 , 𝑚 1 , 𝑚 2 satisfying 𝑗 2 ≥ 𝑗 1 ≥ 0, 𝑗 2 -𝑗 1 ≤ ℓ ≤ 𝑗 2 + 𝑗 1 and -𝑗 𝑖 ≤ 𝑚 𝑖 ≤ 𝑗 𝑖 , 𝑖 = 1 2, the , Clebsch-Gordan (CG) coefficients are

$$C _ { j _ { 1 }, m _ { 1 }, j _ { 2 }, m _ { 2 } } ^ { \ell, m } \coloneqq ( - 1 ) ^ { ( m + j _ { 1 } - j _ { 2 } ) } \sqrt { 2 \ell + 1 } \left ( \begin{array} { c c c } j _ { 1 } & j _ { 2 } & \ell \\ m _ { 1 } & m _ { 2 } & - m \end{array} \right ),$$

see e.g. [44, Chapter 3.5]. The covariant spherical basis vectors are

$$e _ { + 1 } = - \frac { 1 } { \sqrt { 2 } } \left ( [ 1, 0, 0 ] + i [ 0, 1, 0 ] \right ) ^ { T } \quad e _ { 0 } = [ 0, 0, 1 ] ^ { T } \quad e _ { - 1 } = \frac { 1 } { \sqrt { 2 } } \left ( [ 1, 0, 0 ] - i [ 0, 1, 0 ] \right ) ^ { T }.$$

Let

$$c _ { \ell } \coloneqq \sqrt { \frac { \ell + 1 } { 2 \ell + 1 } } \quad d _ { \ell } \coloneqq \sqrt { \frac { \ell } { 2 \ell + 1 } }.$$

We define the coefficients

$$B _ { + 1, \ell, m } & = c _ { \ell } C _ { \ell - 1, m - 1, 1, 1 } ^ { \ell, m } Y _ { \ell - 1, m - 1 } + d _ { \ell } C _ { \ell + 1, m - 1, 1, 1 } ^ { \ell, m } Y _ { \ell + 1, m - 1 } \\ B _ { 0, \ell, m } & = c _ { \ell } C _ { \ell - 1, m, 1, 0 } ^ { \ell, m } Y _ { \ell - 1, m } + d _ { \ell } C _ { \ell + 1, m, 1, 0 } ^ { \ell, m } Y _ { \ell + 1, m } \\ B _ { - 1, \ell, m } & = c _ { \ell } C _ { \ell - 1, m + 1, 1, - 1 } ^ { \ell, m } Y _ { \ell - 1, m + 1 } + d _ { \ell } C _ { \ell + 1, m + 1, 1, - 1 } ^ { \ell, m } Y _ { \ell + 1, m + 1 }$$

and

$$D _ { + 1, \ell, m } = i C _ { \ell, m - 1, 1, 1 } ^ { \ell, m } Y _ { \ell, m - 1 } \quad D _ { 0, \ell, m } = i C _ { \ell, m, 1, 0 } ^ { \ell, m } Y _ { \ell, m } \quad D _ { - 1, \ell, m } = i C _ { \ell, m + 1, 1, - 1 } ^ { \ell, m } Y _ { \ell, m + 1 }.$$

Definition 3.1 (Vector Spherical Harmonics). For ℓ = 1 2 , , . . . , 𝑚 = -ℓ, . . . , ℓ , using the notation of (3), (5) and (6), the divergence-free and curl-free vector spherical harmonics are defined by

$$y _ { \ell, m } = B _ { + 1, \ell, m } e _ { + 1 } + B _ { 0, \ell, m } e _ { 0 } + B _ { - 1, \ell, m } e _ { - 1 } \quad z _ { \ell, m } = D _ { + 1, \ell, m } e _ { + 1 } + D _ { 0, \ell, m } e _ { 0 } + D _ { - 1, \ell, m } e _ { - 1 }.$$

Or equivalently, by (3) and (7),

$$y _ { \ell, m } = \begin{pmatrix} - \frac { 1 } { \sqrt { 2 } } \left ( B _ { + 1, \ell, m } - B _ { - 1, \ell, m } \right ) \\ - \frac { 1 } { \sqrt { 2 } } i \left ( B _ { + 1, \ell, m } + B _ { - 1, \ell, m } \right ) \\ B _ { 0, \ell, m } \end{pmatrix} \quad z _ { \ell, m } = \begin{pmatrix} - \frac { 1 } { \sqrt { 2 } } \left ( D _ { + 1, \ell, m } - D _ { - 1, \ell, m } \right ) \\ - \frac { 1 } { \sqrt { 2 } } i \left ( D _ { + 1, \ell, m } + D _ { - 1, \ell, m } \right ) \\ D _ { 0, \ell, m } \end{pmatrix}.$$

«

‹

«

‹

The set of vector spherical harmonics { y ℓ,𝑚 , z ℓ,𝑚 : ℓ = 1 2 , , . . . , 𝑚 = -ℓ, . . . , ℓ } in (7) or (8) forms an orthonormal basis for L 2 ( S 2 ) . Using the property of 3𝑗 symbols in [10], the CG coefficients in (5) and (6) have the following explicit formula.

$$\text{axis for } L _ { 2 } ( \mathbb { S } ^ { 2 } ). \, \text{Using the property of } 3 - j \text{symbols in } [ 1 0 ], \text{ the CG coefficients in } ( 5 ) \text{ and } ( 6 ) \text{ have the following explicit} \\ \text{ormula.} \\ C _ { t - 1, m } ^ { t, m } - 1, 1, 1 = \sqrt { \frac { ( \ell + m ) ( \ell + m - 1 ) } { ( 2 \ell ) ( 2 \ell - 1 ) } } } \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
 \text{$C_{t}, m } ^ { t, m } - 1, 1, 0 = \sqrt { \frac { ( \ell + m ) ( \ell - m ) } { \ell ( 2 \ell - 1 ) } } \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\text{$C_{t}, m } ^ { t, m } + 1, 1, - 1 = \sqrt { \frac { ( \ell + m ) ( \ell + m - 1 ) } { ( 2 \ell ) ( 2 \ell - 1 ) } } \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\text{$C_{t}, m }^ { t, m } + 1, 1, - 1 = \sqrt { \frac { ( \ell + m ) ( \ell + m - 1 ) } { ( 2 \ell + 2 ) } } \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
\text{$c_{t}, m } ^ { t, m } + 1, 1, - 1 = \sqrt { \frac { ( \ell + m ) ( \ell + m ) } { \ell ( 2 \ell + 2 ) } }. \\ C _ { t, m, 1, 0 = \sqrt { \ell ( \ell + 1 ) } }. \\ \text{The expression of these } C G \text{ coefficients in } ( 9 ) \text{ will simplify computations in } the proposed fast algorithms below, and} \\ \text{dso help with the interested readers follow the routines in our software.}$$

The expression of these CG coefficients in (9) will simplify computations in the proposed fast algorithms below, and also help with the interested readers follow the routines in our software.

## 4 FAST VECTOR SPHERICAL HARMONIC TRANSFORMS

In Subsections 4.1 and 4.2 below, we prove two theorems to represent the vector spherical harmonic transforms by scalar spherical harmonics and Clebsch-Gordan coefficients. From the representation formula, we obtain a computational strategy for fast evaluation of FwdVSHT and AdjVSHT. Subsection 4.3 shows the analysis of computational complexity and approximation error for the proposed FaVeST . In Subsection 4.4, we provide the 'user guide' of the software implementation for FaVeST in Matlab environment.

## 4.1 Fast Computation for FwdVSHT

In this section, we provide an efficient way to evaluate the Fourier coefficients for vector spherical harmonics. The evaluation exploits the connection between the coefficients for vector and scalar spherical harmonics. This connection together with FFTs for scalar spherical harmonics allows fast computation of FwdVSHT and AdjVSHT.

The divergence-free and curl-free coefficients of a tangent field 𝑇 on S 2 are, for ℓ = 1 2 , , . . . , 𝑚 = -ℓ, . . . , ℓ ,

$$\widehat { T } _ { \ell, m } \coloneqq \left \langle T, y _ { \ell, m } \right \rangle \quad \widetilde { T } _ { \ell, m } \coloneqq \left \langle T, z _ { \ell, m } \right \rangle.$$

To numerically evaluate them, one needs to discretize the integrals of the coefficients by a quadrature rule [29], which is a set Q 𝑁 : = {( 𝑤 , 𝑖 x 𝑖 )} 𝑁 𝑖 = 1 of 𝑁 𝑁 , ≥ 2, pairs of real numbers and points on S 2 .

Definition 4.1 (FwdVSHT). For a sequence of column vectors { 𝑇 𝑘 } 𝑁 𝑘 = 1 in R 3 , discrete forward divergence-free and curl-free transforms for { 𝑇 𝑘 } 𝑁 𝑘 = 1 associated with Q 𝑁 , or simply FwdVSHT , are the weighted sums

$$\widehat { F } _ { \ell, m } ( T. ) & \coloneqq \widehat { F } _ { \ell, m } ( \{ T _ { k } \} _ { k = 1 } ^ { N }, Q _ { N } ) \coloneqq \sum _ { k = 1 } ^ { N } w _ { k } y _ { \ell, m } ^ { * } ( x _ { k } ) T _ { k } \\ \widetilde { F } _ { \ell, m } ( T. ) & \coloneqq \widetilde { F } _ { \ell, m } ( \{ T _ { k } \} _ { k = 1 } ^ { N }, Q _ { N } ) \coloneqq \sum _ { k = 1 } ^ { N } w _ { k } z _ { \ell, m } ^ { * } ( x _ { k } ) T _ { k }. \\ \tilde { \Delta } \Delta T _ { k } \, \tilde { \Delta } \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \alpha \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \beta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Dambda \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Lambda \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Dambd \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Depsilon \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Persilon \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Depsilor \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \epsilon \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Damb \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \dambda \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Depsilion \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \depsilon \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Lambd \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Delta \Lambfrak { | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | c | }$$

For a tangent field 𝑇 , the divergence-free and curl-free coefficients can be approximated by FwdVSHT for the sequence of values { 𝑇 ( x 𝑖 )} 𝑁 𝑖 = 1 of the tangent field at quadrature nodes { x 𝑖 } 𝑁 𝑖 = 1 :

$$\cdots \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \cdot \colon \\ \widehat { T } _ { \ell, m } & \approx \widehat { F } _ { \ell, m } ( \{ T ( x _ { k } ) \} _ { k = 1 } ^ { N }, \mathcal { Q } _ { N } ) \coloneqq \sum _ { k = 1 } ^ { N } w _ { k } y _ { \ell, m } ^ { * } ( x _ { k } ) T ( x _ { k } ) \\ \widetilde { T } _ { \ell, m } & \approx \widetilde { F } _ { \ell, m } ( \{ T ( x _ { k } ) \} _ { k = 1 } ^ { N }, \mathcal { Q } _ { N } ) \coloneqq \sum _ { k = 1 } ^ { N } w _ { k } z _ { \ell, m } ^ { * } ( x _ { k } ) T ( x _ { k } ). \\ \text{ve call } ( 1 0 \partial \text{ forand vector theorical harmonic transform or found } \text{or found} V W T F \text{ for tangent field } T$$

As mentioned, we call (10) forward vector spherical harmonic transform , or FwdVSHT for tangent field 𝑇 .

Now, we design fast computation for FwdVSHT in Definition 4.1 using the scalar version of FwdVSHT. Let { 𝑓 𝑘 } 𝑁 𝑘 = 1 be a real sequence. For ℓ ≥ 1 , 𝑚 = -ℓ, . . . , ℓ , the discrete forward scalar spherical harmonic transform (FwdSHT) is

$$F _ { \ell, m } ( f. ) \coloneqq F _ { \ell, m } ( f _ { k } ) \coloneqq F _ { \ell, m } ( \{ f _ { k } \} _ { k = 1 } ^ { N }, \mathcal { Q } _ { N } ) \coloneqq \sum _ { k = 1 } ^ { N } w _ { k } f _ { k } Y _ { \ell, m } ^ { * } ( \mathbf x _ { k } ).$$

They are approximations of spherical harmonic coefficients ⟨ 𝑓 , 𝑌 ℓ,𝑚 ⟩ for 𝑓 by the quadrature rule Q 𝑁 . The following theorem shows a representation of the FwdVSHT by FwdSHT for a sequence { 𝑇 𝑘 } 𝑁 𝑘 = 1 , which would allow us to efficiently compute b F ℓ,𝑚 ( 𝑇 · ) and e F ℓ,𝑚 ( 𝑇 · ) . For each 𝑘 = 1 , . . . , 𝑁 , let GLYPH&lt;0&gt; 𝑇 ( 1 ) 𝑘 , 𝑇 ( 2 ) 𝑘 , 𝑇 ( 3 ) 𝑘 GLYPH&lt;1&gt; be the components of the vector 𝑇 𝑘 . For ℓ = 1 , . . . , 𝑚 = -ℓ, . . . , ℓ , using the notation of (4) and (9), we define the coefficients

$$\dots \, \cfrac { 1 } { 2 } \cfrac { 1 } { 2 } \cfrac { 1 } { 2 } \cdots \, \cfrac { 1 } { 2 } \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \, \cdots \
 \xi ^ { ( 1 ) } _ { \ell, m } \coloneqq c _ { \ell + 1 } C ^ { \ell + 1, m + 1 } _ { \ell, m, 1, 1 } \xi ^ { ( 2 ) } _ { \ell, m } \coloneqq d _ { \ell - 1 } C ^ { \ell - 1, m + 1 } _ { \ell, m, 1, 1 } \\ \xi ^ { ( 3 ) } _ { \ell, m } \coloneqq c _ { \ell + 1 } C ^ { \ell + 1, m - 1 } _ { \ell, m, 1, - 1 } \xi ^ { ( 4 ) } _ { \ell, m } \coloneqq d _ { \ell - 1 } C ^ { \ell - 1, m - 1 } _ { \ell, m, 1, - 1 } \\ \xi ^ { ( 5 ) } _ { \ell, m } \coloneqq c _ { \ell + 1 } C ^ { \ell + 1, m } _ { \ell, m, 1, 0 } \xi ^ { ( 6 ) } _ { \ell, m } \coloneqq d _ { \ell - 1 } C ^ { \ell - 1, m } _ { \ell, m, 1, 0 }$$

and

$$\mu _ { \ell, m } ^ { ( 1 ) } \coloneqq C _ { \ell, m, 1, 1 } ^ { \ell, m + 1 } \quad \mu _ { \ell, m } ^ { ( 2 ) } \coloneqq C _ { \ell, m, 1, 0 } ^ { \ell, m } \quad \mu _ { \ell, m } ^ { ( 3 ) } \coloneqq C _ { \ell, m, 1 - 1 } ^ { \ell, m - 1 }.$$

5

Theorem 4.2. Let { 𝑇 𝑘 } 𝑁 𝑘 = 1 be a sequence of column vectors in R 3 and Q 𝑁 : = {( 𝑤 , 𝑖 x 𝑖 )} 𝑁 𝑖 = 1 a quadrature rule on S 2 . Then, for ℓ = 1 2 , , . . . , 𝑚 = -ℓ, . . . , ℓ , the FwdVSHT can be represented by its scalar version FwdSHT, as follows.

$$& \cdots \kappa = 1 \cdots \ell \cdots \ell \cdots \kappa = 1 \cdots \ell \\ & \text{Then, for } \ell = 1, 2, \dots, m = - \ell, \dots, \ell, \text{ the FwdVSHT can be represented by its scalar version} \ FwdSHT, as follows. \\ & \widehat { F } _ { c, m } ( T. ) = \frac { 1 } { \sqrt { 2 } } \left \{ \xi ^ { ( 1 ) } _ { c ^ { - 1 } - 1, m - 1 } \left [ - \widetilde { F } _ { c ^ { - 1 } - 1, m - 1 } \left ( T ^ { ( 1 ) } _ { c ^ { - 1 } } \right ) + i \widetilde { F } _ { c ^ { - 1 } - 1, m - 1 } \left ( T ^ { ( 2 ) } _ { c ^ { - 1 } } \right ) \right ] + \xi ^ { ( 2 ) } _ { c ^ { - 1 } + 1, m - 1 } \left [ - \widetilde { F } _ { c ^ { - 1 } + 1, m - 1 } \left ( T ^ { ( 1 ) } _ { c ^ { - 1 } } \right ) + i \widetilde { F } _ { c ^ { - 1 } + 1, m - 1 } \left ( T ^ { ( 2 ) } _ { c ^ { - 1 } } \right ) \right ] \\ & \quad + \xi ^ { ( 3 ) } _ { c ^ { - 1 } - 1, m + 1 } \left [ \widetilde { F } _ { c ^ { - 1 } - 1, m + 1 } \left ( T ^ { ( 1 ) } _ { c ^ { - 1 } } \right ) + i \widetilde { F } _ { c ^ { - 1 } - 1, m + 1 } \left ( T ^ { ( 2 ) } _ { c ^ { - 1 } } \right ) \right ] \\ & \quad + \xi ^ { ( 5 ) } _ { c ^ { - 1 } - 1, m } \widetilde { F } _ { c ^ { - 1 } - 1, m } \left ( T ^ { ( 3 ) } _ { c ^ { - 1 } } \right ) + \xi ^ { ( 6 ) } _ { c ^ { - 1 } + 1, m } \widetilde { F } _ { c ^ { - 1 } + 1, m } \left ( T ^ { ( 3 ) } _ { c ^ { - 1 } } \right ) \\ & \widetilde { F } _ { c, m } ( T. ) = - \frac { 1 } { \sqrt { 2 } } \left \{ \mu ^ { ( 1 ) } _ { c ^ { - 1 }, m - 1 } \left ( - \widetilde { F } _ { c, m - 1 } ( T ^ { ( 1 ) } _ { c ^ { - 1 } } ) + i \widetilde { F } _ { c, m - 1 } ( T ^ { ( 2 ) } _ { c ^ { - 1 } } ) \right ) + \mu ^ { ( 3 ) } _ { c ^ { - 1 }, m + 1 } \left ( \widetilde { F } _ { c, m + 1 } ( T ^ { ( 1 ) } _ { c ^ { - 1 } } ) + i \widetilde { F } _ { c, m + 1 } ( T ^ { ( 2 ) } _ { c ^ { - 1 } } ) \right ) \right ] \\ & \quad - i \mu ^ { ( 2 ) } _ { c ^ { - 1 }, m } \widetilde { F } _ { c, m } ( T ^ { ( 3 ) } _ { c ^ { - 1 } } ), \\ & \text{where we use the notation of } ( 10 ), ( 12 ), ( 13 ) \text{ and } ( 14 ).$$

where we use the notation of (10) , (12) , (13) and (14) .

Proof. Let ℓ = 1 2 , , . . . , and 𝑚 = -ℓ, . . . , ℓ . By (5) and (13), the discrete forward divergence-free transform

- b F ℓ,𝑚 ( 𝑇 · ) = 𝑁 ∑︁ 𝑘 = 1 𝑤 𝑘 y ∗ ℓ,𝑚 ( x 𝑘 ) 𝑇 𝑘 = 𝑁 ∑︁ 𝑘 = 1 𝑤 𝑘 GLYPH&lt;18&gt; -√ 1 2 𝑇 ( 1 ) 𝑘 GLYPH&lt;0&gt; 𝐵 ∗ + 1 ,ℓ,𝑚 ( x 𝑘 ) -𝐵 ∗ -1 ,ℓ,𝑚 ( x 𝑘 ) GLYPH&lt;1&gt; + √ 1 2 𝑖𝑇 ( 2 ) 𝑘 GLYPH&lt;0&gt; 𝐵 ∗ + 1 ,ℓ,𝑚 ( x 𝑘 ) + 𝐵 ∗ -1 ,ℓ,𝑚 ( x 𝑘 ) GLYPH&lt;1&gt; + 𝑇 ( 3 ) 𝑘 𝐵 ∗ 0 ,ℓ,𝑚 ( x 𝑘 ) GLYPH&lt;19&gt; = 𝑁 ∑︁ 𝑘 = 1 𝑤 𝑘 GLYPH&lt;20&gt; GLYPH&lt;16&gt; 1 √ 2 GLYPH&lt;0&gt; -𝑇 ( 1 ) 𝑘 + 𝑖𝑇 ( 2 ) 𝑘 GLYPH&lt;1&gt; 𝑐 𝐶 ℓ ℓ,𝑚 ℓ -1 ,𝑚 -1 1 1 , , GLYPH&lt;17&gt; 𝑌 ∗ ℓ -1 ,𝑚 -1 ( x 𝑘 ) + GLYPH&lt;16&gt; 1 √ 2 GLYPH&lt;0&gt; -𝑇 ( 1 ) 𝑘 + 𝑖𝑇 ( 2 ) 𝑘 GLYPH&lt;1&gt; 𝑑 𝐶 ℓ ℓ,𝑚 ℓ + 1 ,𝑚 -1 1 1 , , GLYPH&lt;17&gt; 𝑌 ∗ ℓ + 1 ,𝑚 -1 ( x 𝑘 ) + GLYPH&lt;16&gt; 1 √ 2 GLYPH&lt;0&gt; 𝑇 ( 1 ) 𝑘 + 𝑖𝑇 ( 2 ) 𝑘 GLYPH&lt;1&gt; 𝑐 𝐶 ℓ ℓ,𝑚 ℓ -1 ,𝑚 + 1 1 , , -1 GLYPH&lt;17&gt; 𝑌 ∗ ℓ -1 ,𝑚 + 1 ( x 𝑘 ) + GLYPH&lt;16&gt; 1 √ 2 GLYPH&lt;0&gt; 𝑇 ( 1 ) 𝑘 + 𝑖𝑇 ( 2 ) 𝑘 GLYPH&lt;1&gt; 𝑑 𝐶 ℓ ℓ,𝑚 ℓ + 1 ,𝑚 + 1 1 , , -1 GLYPH&lt;17&gt; 𝑌 ∗ ℓ + 1 ,𝑚 + 1 ( x 𝑘 ) + GLYPH&lt;0&gt; 𝑇 ( 3 ) 𝑘 𝑐 𝐶 ℓ ℓ,𝑚 ℓ -1 ,𝑚, 1 0 , GLYPH&lt;1&gt; 𝑌 ∗ ℓ -1 ,𝑚 ( x 𝑘 ) + GLYPH&lt;0&gt; 𝑇 ( 3 ) 𝑘 𝑑 𝐶 ℓ ℓ,𝑚 ℓ + 1 ,𝑚, 1 0 , GLYPH&lt;1&gt; 𝑌 ∗ ℓ + 1 ,𝑚 ( x 𝑘 ) GLYPH&lt;21&gt; = 1 √ 2 GLYPH&lt;26&gt; 𝜉 ( 1 ) ℓ -1 ,𝑚 -1 h -F ℓ -1 ,𝑚 -1 GLYPH&lt;16&gt; 𝑇 ( 1 ) · GLYPH&lt;17&gt; + 𝑖 F ℓ -1 ,𝑚 -1 GLYPH&lt;16&gt; 𝑇 ( 2 ) · GLYPH&lt;17&gt; i + 𝜉 ( 2 ) ℓ + 1 ,𝑚 -1 h -F ℓ + 1 ,𝑚 -1 GLYPH&lt;16&gt; 𝑇 ( 1 ) · GLYPH&lt;17&gt; + 𝑖 F ℓ + 1 ,𝑚 -1 GLYPH&lt;16&gt; 𝑇 ( 2 ) · GLYPH&lt;17&gt; i + 𝜉 ( 3 ) ℓ -1 ,𝑚 + 1 h F ℓ -1 ,𝑚 + 1 GLYPH&lt;16&gt; 𝑇 ( 1 ) · GLYPH&lt;17&gt; + 𝑖 F ℓ -1 ,𝑚 + 1 GLYPH&lt;16&gt; 𝑇 ( 2 ) · GLYPH&lt;17&gt; i + 𝜉 ( 4 ) ℓ + 1 ,𝑚 + 1 h F ℓ + 1 ,𝑚 + 1 GLYPH&lt;16&gt; 𝑇 ( 1 ) · GLYPH&lt;17&gt; + 𝑖 F ℓ + 1 ,𝑚 + 1 GLYPH&lt;16&gt; 𝑇 ( 2 ) · GLYPH&lt;17&gt; i GLYPH&lt;27&gt; + 𝜉 ( 5 ) ℓ -1 ,𝑚 F ℓ -1 ,𝑚 GLYPH&lt;16&gt; 𝑇 ( 3 ) · GLYPH&lt;17&gt; + 𝜉 ( 6 ) ℓ + 1 ,𝑚 F ℓ + 1 ,𝑚 GLYPH&lt;16&gt; 𝑇 ( 3 ) · GLYPH&lt;17&gt; .

In a similar way, for the curl-free case, we use (6) and (14) to obtain

$$\widetilde { F } _ { \ell, m } ( T. ) \coloneqq \sum _ { k = 1 } ^ { N } w _ { k } z _ { \ell, m } ^ { * } ( { \mathbf x } _ { k } ) T ( { \mathbf x } _ { k } )$$

$$& = \sum _ { k = 1 } ^ { N } w _ { k } \left \{ - \frac { 1 } { \sqrt { 2 } } i \left [ \mu _ { \ell, m - 1 } ^ { ( 1 ) } Y _ { \ell, m - 1 } ^ { * } ( \mathbf x _ { k } ) \left ( - T _ { k } ^ { ( 1 ) } + i T _ { k } ^ { ( 2 ) } \right ) + \mu _ { \ell, m + 1 } ^ { ( 3 ) } Y _ { \ell, m + 1 } ^ { * } ( \mathbf x _ { k } ) \left ( T _ { k } ^ { ( 1 ) } + i T _ { k } ^ { ( 2 ) } \right ) \right ] \\ & \quad \ \ - i \mu _ { \ell, m } ^ { ( 2 ) } Y _ { \ell, m } ^ { * } ( \mathbf x _ { k } ) T _ { k } ^ { ( 3 ) } \right \} \\ & = - \frac { 1 } { \sqrt { 2 } } i \left [ \mu _ { \ell, m - 1 } ^ { ( 1 ) } \left ( - F _ { \ell, m - 1 } ( T _ { ^ { 1 } } ) + i \, F _ { \ell, m - 1 } ( T _ { ^ { 2 } } ) \right ) + \mu _ { \ell, m + 1 } ^ { ( 3 ) } \left ( F _ { \ell, m + 1 } ( T _ { ^ { 1 } } ) + i \, F _ { \ell, m + 1 } ( T _ { ^ { 2 } } ) \right ) \right ] \\ & \quad \ \ - i \mu _ { \ell, m } ^ { ( 2 ) } F _ { \ell, m } ( T _ { ^ { 1 } } ),$$

thus completing the proof.

## 4.2 Fast Computation for AdjVSHT

In this section, we study the adjoint vector spherical harmonic transforms for a tangent field on the sphere.

Definition 4.3 (AdjVSHT). For 𝐿 ≥ 1 and two complex sequences { 𝑎 ℓ,𝑚 , 𝑏 ℓ,𝑚 : ℓ = 1 2 , , . . . , 𝑚 = -ℓ, . . . , ℓ } , the adjoint vector spherical harmonic transform or AdjVSHT of degree 𝐿 is the Fourier partial sum

$$S _ { L } ( a _ { \ell, m }, b _ { \ell, m } ; x ) \coloneqq & \sum _ { \ell = 1 } ^ { L } \sum _ { m = - \ell } ^ { \ell } \left ( a _ { \ell, m } y _ { \ell, m } ( x ) + b _ { \ell, m } z _ { \ell, m } ( x ) \right ), \quad x \in \mathbb { S } ^ { 2 }.$$

Or equivalently, by (3) and (7),

$$S _ { L } ( a _ { \ell, m }, b _ { \ell, m } ) & = \sum _ { \ell = 1 } ^ { L } \sum _ { m = - \ell } ^ { t } \begin{pmatrix} - \frac { 1 } { \sqrt { 2 } } \left ( a _ { \ell, m } B _ { + 1, \ell, m } - a _ { \ell, m } B _ { - 1, \ell, m } \right ) \\ - \frac { 1 } { \sqrt { 2 } } \left ( a _ { \ell, m } B _ { + 1, \ell, m } + a _ { \ell, m } B _ { - 1, \ell, m } \right ) \\ a _ { \ell, m } B _ { 0, \ell, m } \end{pmatrix} + \sum _ { \ell = 1 } ^ { L } \sum _ { m = - \ell } ^ { t } \begin{pmatrix} - \frac { 1 } { \sqrt { 2 } } \left ( b _ { \ell, m } D _ { + 1, \ell, m } - b _ { \ell, m } D _ { - 1, \ell, m } \right ) \\ - \frac { 1 } { \sqrt { 2 } } \left ( b _ { \ell, m } D _ { + 1, \ell, m } + b _ { \ell, m } D _ { - 1, \ell, m } \right ) \\ b _ { \ell, m } D _ { 0, \ell, m }. \end{pmatrix}. \quad ( 1 7 ) \\ - \quad - \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdot \quad \cdots$$

For 𝐿 ≥ 0 and a finite complex sequence 𝑔 ℓ,𝑚 , ℓ = 0 1 , , . . . , 𝐿 , 𝑚 = -ℓ, . . . , ℓ , the adjoint scalar spherical harmonic transform or AdjSHT of degree 𝐿 is the Fourier partial sum of scalar spherical harmonics 𝑌 ℓ,𝑚 :

$$S _ { L } ( g _ { \ell, m } ) \coloneqq \sum _ { \ell = 0 } ^ { L } \sum _ { m = - \ell } ^ { \ell } g _ { \ell, m } Y _ { \ell, m }.$$

The AdjVSHT can be represented by its scalar version with CG coefficients, as shown by the following theorem. With the notation of (4) and (9), we define the coefficients

$$& \text{The AdJVSHT can be represented by its scalar version with CG coefficients, as shown by the following theorem.} \\ & \text{With the notation of } ( 4 ) \, \text{and} \, ( 9 ), \, \text{we define the coefficients} \\ & \text{$v_{t},m \colon = c_{t+1} \left ( a_{t+1},m+1 C_{t},m,1,1}^{c_{t+1},m-1}^{c_{t+1},m-1})$, \quad \, \text{$\ell=0,\dots,L-1,\ m=-\ell,\dots,\ell$,} \\ & \text{$v_{t},m \colon = \left \{ \, d_{t-1} \left ( a_{t-1},m+1 C_{t},m,1,1}^{c_{t-1},m-1}^{c_{t-1},m-1})$, \quad \, \text{$\ell=2,\dots,L+1,\ lm|=0,1,\dots,\ell-2,\dots,\ell$,} \\ & \text{$v_{t},m \colon = i c_{t+1} \left ( a_{t+1},m+1 C_{t},m,1,1}^{c_{t-1},m-1})$, \quad \, \text{$\ell=0,\dots,L-1,\ m=-\ell,\dots,\ell$,} \\ & \text{$v_{t},m \colon = \left \{ \, \text{id}_{t-1} \left ( a_{t-1},m+1 C_{t},m,1,1}^{c_{t-1},m-1})$, \quad \, \text{$\ell=2,\dots,L+1,\ m=-\ell,\dots,\ell$,} \\ & \text{$v_{t},m \colon = \left \{ \, \text{id}_{t-1} \left ( a_{t-1},m+1 C_{t},m,1,1}^{c_{t-1},m-1})$, \quad \, \text{$\ell=2,\dots,L+1,\ m=-\ell,\dots,\ell$,} \\ & \text{$v_{t},m \colon = a_{t+1},m c_{t+1}^{c_{t},m,1,0}^{c_{t-1},m-1}^{c_{t-1},m-1})$, \quad \, \text{$\ell=0,\dots,L-1,\ m=-\ell,\dots,\ell$,} \\ & \text{$v_{t},m \colon = \left \{ \, \text{id}_{t-1} \left ( a_{t-1},m+1 C_{t},m,1,0}^{c_{t-1},m-1}^{c_{t-1},m-1})$, \quad \, \text{$\ell=0,\dots,L-1,\ m=-\ell,\dots,\ell$,} \\ & \text{$v_{t},m \colon = \left \{ \, \text{id}_{t-1} \left ( a_{t-1},m+1 C_{t},m,1,0}^{c_{t-1},m-1})$, \quad \, \text{$\ell=0,\dots,L-1,\ m=-\ell,\dots,\ell$,}$$

□

and

$$\mathfrak { d } _ { \eta _ { \ell, m } ^ { ( 1 ) } } \coloneqq \begin{cases} \ i \left ( b _ { \ell, m + 1 } C _ { \ell, m, 1, 1 } ^ { \ell, m + 1 } - b _ { \ell, m - 1 } C _ { \ell, m, 1, - 1 } ^ { \ell, m - 1 } \right ), & \ell = 1, \dots, L, \ m = - \ell, \dots, \ell, \\ 0, & \ell = 0, \ m = 0, \end{cases} \\ \eta _ { \ell, m } ^ { ( 2 ) } \coloneqq \begin{cases} \ b _ { \ell, m + 1 } C _ { \ell, m, 1, 1 } ^ { \ell, m + 1 } + b _ { \ell, m - 1 } C _ { \ell, m, 1, - 1 } ^ { \ell, m - 1 }, & \ell = 1, \dots, L, \ m = - \ell, \dots, \ell, \\ 0, & \ell = 0, \ m = 0, \end{cases} \\ \mu _ { \ell, m } ^ { ( 3 ) } \coloneqq \begin{cases} \ i b _ { \ell, m } C _ { \ell, m, 1, 0 } ^ { \ell, m }, & \ell = 1, \dots, L, \ m = - \ell, \dots, \ell, \\ 0, & \ell = 0, \ m = 0. \end{cases} \\ \nu _ { \mu } ^ { \dots } \dots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots$$

Theorem 4.4. Let { 𝑎 ℓ,𝑚 , 𝑏 ℓ,𝑚 : ℓ = 1 2 , , . . . , 𝑚 = -ℓ, . . . , ℓ } be two complex sequences. For 𝐿 ≥ 1 , the AdjVSHT for 𝑎 ℓ,𝑚 , 𝑏 ℓ,𝑚 can be represented by its scalar version AdjSHT, as follows.

$$\dots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \cdots \\ S _ { L } ( a _ { \ell, m }, b _ { \ell, m } ) = \begin{pmatrix} - \frac { 1 } { \sqrt { 2 } } \left ( S _ { L - 1 } ( v _ { \ell, m } ^ { ( 1 ) } ) + S _ { L + 1 } ( v _ { \ell, m } ^ { ( 2 ) } ) + S _ { L } ( \eta _ { \ell, m } ^ { ( 1 ) } ) \right ) \\ - \frac { 1 } { \sqrt { 2 } } \left ( S _ { L - 1 } ( v _ { \ell, m } ^ { ( 3 ) } ) + S _ { L + 1 } ( v _ { \ell, m } ^ { ( 4 ) } ) - S _ { L } ( \eta _ { \ell, m } ^ { ( 2 ) } ) \right ) \\ S _ { L - 1 } ( v _ { \ell, m } ^ { ( 5 ) } ) + S _ { L + 1 } ( v _ { \ell, m } ^ { ( 6 ) } ) + S _ { L } ( \eta _ { \ell, m } ^ { ( 3 ) } ) \end{pmatrix}, \\ \ e t h e n o t a t i o n \, \text{of} \, ( 1 8 ), ( 1 9 ), ( 2 0 ) \, \text{and} \, ( 2 1 ).$$

where we use the notation of (18) , (19) , (20) and (21) .

Proof. By (5), we write each component of (17) by AdjSHT, as follows. For the divergence-free term in (17),

$$& \text{PROF.} \text{ By (5), we write each component of (17) by AdjSHT, as follows. For the divergence-free term in (17), \\ & \quad - \frac { 1 } { \sqrt { 2 } } \sum _ { t = 1 } ^ { L } \sum _ { m = - t } ^ { t } \left ( a _ { t, m } B + 1, t, m - a _ { t, m } B - 1, t, m \right ) \\ & \quad = - \frac { 1 } { \sqrt { 2 } } \left [ \sum _ { t = 1 } ^ { L } \sum _ { m = - t } ^ { t } a _ { t, m } \left ( c _ { t } C _ { t - 1, m - 1, 1, 1 } \right ) Y _ { t - 1, m - 1 } + \sum _ { t = 1 } ^ { L } \sum _ { m = - t } ^ { t } a _ { t, m } \left ( d _ { t } C _ { t + 1, m - 1, 1, 1 } \right ) Y _ { t + 1, m - 1 } \\ & \quad - \sum _ { t = 1 } ^ { L } \sum _ { m = - t } ^ { t } a _ { t, m } \left ( c _ { t } C _ { t - 1, m + 1, 1, - 1 } \right ) Y _ { t - 1, m + 1 } - \sum _ { t = 1 } ^ { L } \sum _ { m = - t } ^ { t } a _ { t, m } \left ( d _ { t } C _ { t + 1, m + 1, 1, - 1 } \right ) Y _ { t + 1, m + 1 } \right ] \\ & \quad = - \frac { 1 } { \sqrt { 2 } } \left [ \sum _ { t = 0 } ^ { L - 1 } \sum _ { m = - t } ^ { t } a _ { t + 1, m + 1 } \left ( c _ { t + 1 } C _ { t + 1, m + 1 } \right ) Y _ { t, m } + \sum _ { t = 2 } ^ { L + 1 } \sum _ { m = - t } ^ { t - 2 } a _ { t - 1, m + 1 } \left ( d _ { t } C _ { t - 1 } C _ { t, m, 1, 1 } \right ) Y _ { t, m } \\ & \quad - \sum _ { t = 0 } ^ { L - 1 } \sum _ { m = - t } ^ { t + 2 } a _ { t + 1, m - 1 } \left ( c _ { t + 1 } C _ { t, m, 1, - 1 } C _ { t, m, 1, - 1 } \right ) Y _ { t, m } - \sum _ { t = 2 } ^ { L + 1 } \sum _ { m = - t } ^ { t } a _ { t - 1, m - 1 } \left ( d _ { t - 1 } C _ { t, m, 1, - 1 } C _ { t, m, 1, - 1 } \right ) Y _ { t, m } \right ] \\ & \text{This and (9) give} \\ & \quad. \quad L \quad t$$

This and (9) give

$$& \text{and} \left ( 9 \right ) \text{give} \\ & \quad - \frac { 1 } { \sqrt { 2 } } \sum _ { t = 1 } ^ { L } \sum _ { m = - t } ^ { t } \left ( a _ { t, m } B _ { + 1, t, m } - a _ { t, m } B _ { - 1, t, m } \right ) \\ & \quad = - \frac { 1 } { \sqrt { 2 } } \left [ \sum _ { t = 0 } ^ { L - 1 } \sum _ { m = - t - 2 } ^ { t + 2 } c _ { t + 1 } \left ( a _ { t + 1, m + 1 } C _ { t, m, 1, 1 } ^ { t + 1, m + 1 } - a _ { t + 1, m - 1 } C _ { t, m, 1, - 1 } ^ { t + 1, m - 1 } \right ) Y _ { t, m } \\ & \quad + \sum _ { t = 2 } ^ { L + 1 } \sum _ { m = - t + 2 } ^ { t - 2 } d _ { t - 1 } \left ( a _ { t - 1, m + 1 } C _ { t, m, 1, 1 } ^ { t - 1, m + 1 } - a _ { t - 1, m - 1 } C _ { t, m, 1, - 1 } ^ { t - 1, m - 1 } \right ) Y _ { t, m } \right ], \\ & \text{for the} \, & \left | m \right | > t \left ( \text{which exceeds the range of $m$ for spherical harmonics} \right ), Y _ { t, m } = 0. \, \text{We then let}$$

where for the | 𝑚 | &gt; ℓ (which exceeds the range of 𝑚 for spherical harmonics), 𝑌 ℓ,𝑚 = 0. We then let

$$v _ { \ell, m } ^ { ( 1 ) } \coloneqq c _ { \ell + 1 } \left ( a _ { \ell + 1, m + 1 } C _ { \ell, m, 1, 1 } ^ { \ell + 1, m + 1 } - a _ { \ell + 1, m - 1 } C _ { \ell, m, 1, - 1 } ^ { \ell + 1, m - 1 } \right ), \quad \ell = 0, \dots, L - 1, \ m = - \ell, \dots, \ell,$$

Similarly,

$$y, \\ & - \frac { 1 } { \sqrt { 2 } } i \sum _ { t = 1 } ^ { L } \sum _ { m = - t } ^ { t } \left ( a _ { t, m } B _ { + 1, t, m } + a _ { t, m } B _ { - 1, t, m } \right ) \\ & = - \frac { 1 } { \sqrt { 2 } } \left [ \sum _ { t = 0 } ^ { L - 1 } \sum _ { m = - t - 2 } ^ { \ell + 2 } i c _ { t + 1 } \left ( a _ { t + 1, m + 1 } C _ { t, m, 1, 1 } ^ { \ell + 1, m + 1 } + a _ { t + 1, m - 1 } C _ { t, m, 1, - 1 } ^ { \ell + 1, m - 1 } \right ) Y _ { t, m } \\ & \quad + \sum _ { t = 2 } ^ { L + 1 } \sum _ { m = - t + 2 } ^ { \ell - 2 } i d _ { t - 1 } \left ( a _ { t - 1, m + 1 } C _ { t, m, 1, 1 } ^ { \ell - 1, m + 1 } + a _ { t - 1, m - 1 } C _ { t, m, 1, - 1 } ^ { \ell - 1, m - 1 } \right ) Y _ { t, m } \right ].$$

$$v _ { \ell, m } ^ { ( 3 ) } \coloneqq i c _ { \ell + 1 } \left ( a _ { \ell + 1, m + 1 } C _ { \ell, m, 1, 1 } ^ { \ell + 1, m + 1 } + a _ { \ell + 1, m - 1 } C _ { \ell, m, 1, - 1 } ^ { \ell + 1, m - 1 } \right ), \quad \ell = 0, \dots, L - 1, \ m = - \ell, \dots, \ell, \\ v _ { \ell, m } ^ { ( 4 ) } \coloneqq \begin{cases} \ i d _ { \ell - 1 } \left ( a _ { \ell - 1, m + 1 } C _ { \ell, m, 1, 1 } ^ { \ell - 1, m + 1 } + a _ { \ell - 1, m - 1 } C _ { \ell, m, 1, - 1 } ^ { \ell - 1, m - 1 } \right ), \quad \ell = 2, \dots, L + 1, \ m = - \ell, \dots, \ell, \\ 0, & \ell = 0, 1 \text{ or } | m | = \ell - 2, \ \ell - 1, \end{cases}$$

$$- \frac { 1 } { \sqrt { 2 } } i \sum _ { \ell = 1 } ^ { L } \sum _ { m = - \ell } ^ { \ell } \left ( a _ { \ell, m } B _ { + 1, \ell, m } + a _ { \ell, m } B _ { - 1, \ell, m } \right ) = - \frac { 1 } { \sqrt { 2 } } \left ( S _ { L - 1 } ( v _ { \ell, m } ^ { ( 3 ) } ) + S _ { L + 1 } ( v _ { \ell, m } ^ { ( 4 ) } ) \right ).$$

$$\ e s _ { \ell = 1 } ^ { L } \sum _ { m = - \ell } ^ { \ell } a _ { \ell, m } B _ { 0, \ell, m } = \sum _ { \ell = 0 } ^ { L - 1 } \sum _ { m = - \ell } ^ { \ell } a _ { \ell + 1, m } c _ { \ell + 1 } C _ { \ell, m, 1, 0 } ^ { \ell + 1, m } Y _ { \ell, m } + \sum _ { \ell = 2 } ^ { L + 1 } \sum _ { m = - \ell + 1 } ^ { \ell - 1 } a _ { \ell - 1, m } d _ { \ell - 1 } C _ { \ell, m, 1, 0 } ^ { \ell - 1, m } Y _ { \ell, m },$$

$$& \nu _ { \ell, m } ^ { ( 5 ) } \coloneqq a _ { \ell + 1, m } c _ { \ell + 1 } C _ { \ell, m, 1, 0 } ^ { \ell + 1, m }, \quad \ell = 0, \dots, L - 1, \ m = - \ell, \dots, \ell, \\ & \nu _ { \ell, m } ^ { ( 6 ) } \coloneqq \begin{cases} \ a _ { \ell - 1, m } d _ { \ell - 1 } C _ { \ell, m, 1, 0 } ^ { \ell - 1, m }, \quad \ell = 2, \dots, L + 1, \ | m | = 0, 1, \dots, \ell - 1, \\ 0, \quad \ell = 0, 1 \ \text{or} \ | m | = \ell, \end{cases}$$

then,

$$\sum _ { \ell = 1 } ^ { L } \sum _ { m = - \ell } ^ { \ell } a _ { \ell, m } B _ { 0, \ell, m } = S _ { L - 1 } ( \nu _ { \ell, m } ^ { ( 5 ) } ) + S _ { L + 1 } ( \nu _ { \ell, m } ^ { ( 6 ) } ).$$

For the curl-free term in (17),

$$\cdots \cdots & \cdots \\ - & \frac { 1 } { \sqrt { 2 } } \sum _ { \ell = 1 } ^ { L } \sum _ { m = - \ell } ^ { \ell } \left ( b _ { \ell, m } D _ { + 1, \ell, m } - b _ { \ell, m } D _ { - 1, \ell, m } \right ) \\ & = - \frac { 1 } { \sqrt { 2 } } \left ( \sum _ { \ell = 1 } ^ { L } \sum _ { m = - \ell } ^ { \ell } i b _ { \ell, m } C _ { \ell, m - 1, 1, 1 } ^ { \ell, m } Y _ { \ell, m - 1 } - \sum _ { \ell = 1 } ^ { L } \sum _ { m = - \ell } ^ { \ell } i b _ { \ell, m } C _ { \ell, m + 1, 1, - 1 } ^ { \ell, m } Y _ { \ell, m + 1 } \right )$$

Let then

As we let

$$v _ { \ell, m } ^ { ( 2 ) } \coloneqq \begin{cases} \ d _ { \ell - 1 } \left ( a _ { \ell - 1, m + 1 } C _ { \ell, m, 1, 1 } ^ { \ell - 1, m + 1 } - a _ { \ell - 1, m - 1 } C _ { \ell, m, 1, - 1 } ^ { \ell - 1, m - 1 } \right ), & \ell = 2, \dots, L + 1, \ | m | = 0, 1, \dots, \ell - 2, \\ 0, & \ell = 0, 1 \ \text{or} \ | m | = \ell - 1, \ \ell, \end{cases}$$

by which and (23),

$$- \frac { 1 } { \sqrt { 2 } } \sum _ { t = 1 } ^ { L } \sum _ { m = - t } ^ { t } \left ( a _ { t, m } B _ { + 1, t, m } - a _ { t, m } B _ { - 1, t, m } \right ) = - \frac { 1 } { \sqrt { 2 } } \left ( S _ { L - 1 } ( v _ { t, m } ^ { ( 1 ) } ) + S _ { L + 1 } ( v _ { t, m } ^ { ( 2 ) } ) \right ).$$

$$= - \frac { 1 } { \sqrt { 2 } } \left ( \sum _ { \ell = 1 } ^ { L } \sum _ { m = - \ell - 1 } ^ { \ell - 1 } i b _ { \ell, m + 1 } C _ { \ell, m, 1, 1 } ^ { \ell, m + 1 } Y _ { \ell, m } - \sum _ { \ell = 1 } ^ { L } \sum _ { m = - \ell + 1 } ^ { \ell + 1 } i b _ { \ell, m - 1 } C _ { \ell, m, 1, - 1 } ^ { \ell, m - 1 } Y _ { \ell, m } \right )$$

$$= - \frac { 1 } { \sqrt { 2 } } \sum _ { \ell = 1 } ^ { L } \sum _ { m = - \ell } ^ { \ell } i \left ( b _ { \ell, m + 1 } C _ { \ell, m, 1, 1 } ^ { \ell, m + 1 } - b _ { \ell, m - 1 } C _ { \ell, m, 1, - 1 } ^ { \ell, m - 1 } \right ) Y _ { \ell, m },$$

where we use (4) and (9). Let

$$\eta _ { \ell, m } ^ { ( 1 ) } \coloneqq \begin{cases} \ i \left ( b _ { \ell, m + 1 } C _ { \ell, m, 1, 1 } ^ { \ell, m + 1 } - b _ { \ell, m - 1 } C _ { \ell, m, 1, - 1 } ^ { \ell, m - 1 } \right ), & \ell = 1, \dots, L, \ m = - \ell, \dots, \ell, \\ 0, & \ell = 0, \ m = 0, \end{cases}$$

we then obtain the equation

$$- \frac { 1 } { \sqrt { 2 } } \sum _ { \ell = 1 } ^ { L } \sum _ { m = - \ell } ^ { \ell } \left ( b _ { \ell, m } D _ { + 1, \ell, m } - b _ { \ell, m } D _ { - 1, \ell, m } \right ) = - \frac { 1 } { \sqrt { 2 } } \, S _ { L } ( \eta _ { \ell, m } ^ { ( 1 ) } ).$$

In a similar way,

$$- \frac { 1 } { \sqrt { 2 } } i \sum _ { \ell = 1 } ^ { L } \sum _ { m = - \ell } ^ { t } \left ( b _ { \ell, m } D _ { + 1, \ell, m } + b _ { \ell, m } D _ { - 1, \ell, m } \right ) = \frac { 1 } { \sqrt { 2 } } \sum _ { \ell = 1 } ^ { L } \sum _ { m = - \ell } ^ { t } \left ( b _ { \ell, m + 1 } C _ { \ell, m, 1, 1 } ^ { \ell, m + 1 } + b _ { \ell, m - 1 } C _ { \ell, m, 1, - 1 } ^ { \ell, m - 1 } \right ) Y _ { \ell, m }.$$

If letting

$$\eta _ { \ell, m } ^ { ( 2 ) } \coloneqq \begin{cases} \ \left ( b _ { \ell, m + 1 } C _ { \ell, m, 1, 1 } ^ { \ell, m + 1 } + b _ { \ell, m - 1 } C _ { \ell, m, 1, - 1 } ^ { \ell, m - 1 } \right ), & \ell = 1, \dots, L, \ m = - \ell, \dots, \ell, \\ 0, & \ell = 0, \ m = 0, \end{cases}$$

we then obtain

As we let

$$\eta _ { \ell, m } ^ { ( 3 ) } \coloneqq \begin{cases} \ i b _ { \ell, m } C _ { \ell, m, 1, 0 } ^ { \ell, m }, & \ell = 1, \dots, L, \ m = - \ell, \dots, \ell, \\ 0, & \ell = 0, \ m = 0, \end{cases}$$

one then has the representation ˝ 𝐿 ℓ = 1 ˝ ℓ 𝑚 = -ℓ 𝑏 ℓ,𝑚 𝐷 0 ,ℓ,𝑚 = 𝑆 𝐿 ( 𝜂 ( 3 ) ℓ,𝑚 ) . Thus,

$$\text{the representation $\sum_{t=1}^{L} \sum_{m=-t}^{t} b_{t,m} D_{0,t,m} = S_{L}(\eta_{t,m})$. Thus,} \\ S _ { L } ( a _ { t, m }, b _ { t, m } ) = \begin{pmatrix} - \frac { 1 } { \sqrt { 2 } } \left ( S _ { L - 1 } ( v _ { t, m } ^ { ( 1 ) } ) + S _ { L + 1 } ( v _ { t, m } ^ { ( 2 ) } ) \right ) \\ - \frac { 1 } { \sqrt { 2 } } \left ( S _ { L - 1 } ( v _ { t, m } ^ { ( 3 ) } ) + S _ { L + 1 } ( v _ { t, m } ^ { ( 4 ) } ) \right ) \\ S _ { L - 1 } ( v _ { t, m } ^ { ( 5 ) } ) + S _ { L + 1 } ( v _ { t, m } ^ { ( 6 ) } ) \end{pmatrix} + \begin{pmatrix} - \frac { 1 } { \sqrt { 2 } } S _ { L } ( \eta _ { t, m } ^ { ( 1 ) } ) \\ \frac { 1 } { \sqrt { 2 } } S _ { L } ( \eta _ { t, m } ^ { ( 3 ) } ) \\ S _ { L } ( \eta _ { t, m } ^ { ( 3 ) } ) \end{pmatrix} \\ = \begin{pmatrix} - \frac { 1 } { \sqrt { 2 } } \left ( S _ { L - 1 } ( v _ { t, m } ^ { ( 1 ) } ) + S _ { L + 1 } ( v _ { t, m } ^ { ( 2 ) } ) + S _ { L } ( \eta _ { t, m } ^ { ( 1 ) } ) \right ) \\ - \frac { 1 } { \sqrt { 2 } } \left ( S _ { L - 1 } ( v _ { t, m } ^ { ( 3 ) } ) + S _ { L + 1 } ( v _ { t, m } ^ { ( 4 ) } ) - S _ { L } ( \eta _ { t, m } ^ { ( 2 ) } ) \right ) \\ S _ { L - 1 } ( v _ { t, m } ^ { ( 5 ) } ) + S _ { L + 1 } ( v _ { t, m } ^ { ( 6 ) } ) + S _ { L } ( \eta _ { t, m } ^ { ( 3 ) } ) \end{pmatrix}, \\ \text{mpletes the proof.}$$

which then completes the proof.

$$- \frac { 1 } { \sqrt { 2 } } i \sum _ { \ell = 1 } ^ { L } \sum _ { m = - \ell } ^ { \ell } ( b _ { \ell, m } D _ { + 1, \ell, m } + b _ { \ell, m } D _ { - 1, \ell, m } ) = \frac { 1 } { \sqrt { 2 } } \, S _ { L } ( \eta _ { \ell, m } ^ { ( 2 ) } ).$$

$$\sum _ { \ell = 1 } ^ { L } \sum _ { m = - \ell } ^ { \ell } b _ { \ell, m } D _ { 0, \ell, m } = \sum _ { \ell = 1 } ^ { L } \sum _ { m = - \ell } ^ { \ell } i b _ { \ell, m } C ^ { \ell, m } _ { \ell, m, 1, 0 } Y _ { \ell, m },$$

## 4.3 Algorithms and Errors

4.3.1 Fast algorithms. In Algorithms 1 and 2, we write down the algorithms for FwdVSHT and AdjVSHT from Theorems 4.2 and 4.4. They achieve fast computation by FFTs for scalar spherical harmonics. In this paper, we use the NFFT package [36] for fast scalar spherical harmonic transforms on S 2 which have the computational cost O GLYPH&lt;16&gt; 𝑁 log √ 𝑁 GLYPH&lt;17&gt; and O GLYPH&lt;16&gt; 𝑀 log √ 𝑀 GLYPH&lt;17&gt; for 𝑁 evaluation points and 𝑀 (Fourier) coefficients.

## ALGORITHM 1: Forward FaVeST

Input: A sequence { 𝑇 , . . . , 𝑇 1 𝑁 } ⊂ R 3 , 𝑁 ≥ 2 and a quadrature rule Q 𝑁 : = {( 𝑤 , 𝑖 x 𝑖 ) } 𝑁 𝑖 = 1 on S 2 , and maximal degree 𝐿 𝐿 , ≥ 1. Output: Complex sequences of AdjVSHT b F ℓ,𝑚 ( 𝑇 𝑘 ) and e F ℓ,𝑚 ( 𝑇 𝑘 ) in (10), ℓ = 1 , . . . , 𝐿, 𝑚 = - ℓ, . . . , ℓ .

Step 1 Compute the AdjSHT F ℓ,𝑚 GLYPH&lt;16&gt; - · 𝑇 ( 1 ) + 𝑖𝑇 ( 2 ) · GLYPH&lt;17&gt; , F ℓ,𝑚 GLYPH&lt;16&gt; 𝑇 ( 1 ) · + 𝑖𝑇 ( 2 ) · GLYPH&lt;17&gt; and F ℓ,𝑚 GLYPH&lt;16&gt; 𝑇 ( 3 ) · GLYPH&lt;17&gt; for ℓ = 0 , . . . , 𝐿 + 1, 𝑚 = - ℓ, . . . , ℓ , by forward FFT for scalar spherical harmonics.

Step 2 Compute 𝜉 ( 𝑖 ) ℓ,𝑚 , 𝑖 = 1 2 , , . . . , 6 and 𝜇 ( 𝑗 ) ℓ,𝑚 , 𝑗 = 1 2 3, for , , ℓ = 0 , . . . , 𝐿 + 1, 𝑚 = - ℓ, . . . , ℓ . Step 3 Compute b F ℓ,𝑚 ( 𝑇 · ) and e F ℓ,𝑚 ( 𝑇 · ) for ℓ = 1 , . . . , 𝐿, 𝑚 = - ℓ, . . . , ℓ , by (15).

## ALGORITHM 2: Adjoint FaVeST

Input: Two complex sequences of coefficients {( 𝑎 ℓ,𝑚 , 𝑏 ℓ,𝑚 ) : ℓ = 1 , . . . , 𝐿, 𝑚 = - ℓ, . . . , ℓ } , and 𝑀 evaluation points { x 𝑖 } 𝑀 𝑖 = 1 . Output: Complex sequence of AdjVSHT S 𝐿 ( 𝑎 ℓ,𝑚 , 𝑏 ℓ,𝑚 ; x 𝑖 ) in (16), 𝑖 = 1 , . . . , 𝑀 . For each point , 𝑖 1 , . . . , 𝑀 , do the following steps.

x 𝑖 =

Step 1 Compute 𝜈 ( 𝑖 ) ℓ,𝑚 , 𝜂 ( 𝑗 ) ℓ,𝑚 , 𝑖 = 1 , . . . , 6, 𝑗 = 1 2 3, by (19), (20) and (21). , ,

Step 2 Evaluate 𝑆 𝐿 ( 𝜈 ( 1 ) ℓ,𝑚 ) , 𝑆 𝐿 + 1 ( 𝜈 ( 2 ) ℓ,𝑚 ) , 𝑆 𝐿 - 1 ( 𝜈 ( 3 ) ℓ,𝑚 ) , 𝑆 𝐿 + 1 ( 𝜈 ( 4 ) ℓ,𝑚 ) , 𝑆 𝐿 - 1 ( 𝜈 ( 5 ) ℓ,𝑚 ) , 𝑆 𝐿 + 1 ( 𝜈 ( 6 ) ℓ,𝑚 ) , 𝑆 𝐿 ( 𝜂 ( 1 ) ℓ,𝑚 ) , 𝑆 𝐿 ( 𝜂 ( 2 ) ℓ,𝑚 ) and 𝑆 𝐿 ( 𝜂 ( 3 ) ℓ,𝑚 ) , by adjoint FFT for scalar spherical harmonics.

Step 3 Compute S 𝐿 ( 𝑎 ℓ,𝑚 , 𝑏 ℓ,𝑚 ; x 𝑖 ) , by (22).

In Algorithm 1, using the forward FFT for scalar spherical harmonics, Step 1 uses 3 times scalar FFTs with degree up to 𝐿 + 1 and then has cost O GLYPH&lt;16&gt; 𝑁 log √ 𝑁 GLYPH&lt;17&gt; (assuming using 𝑁 = O GLYPH&lt;0&gt; 𝐿 2 GLYPH&lt;1&gt; points); Step 2 is the direct computation of 𝜉 ( 𝑖 ) ℓ,𝑚 and 𝜇 ( 𝑗 ) ℓ,𝑚 by (13), (14) and (9) for degree ℓ at most 𝐿 + 1 and then has cost O GLYPH&lt;16&gt; √ 𝑁 GLYPH&lt;17&gt; ; Step 3 evaluates (15) using the results of Steps 1 and 2 and has cost O ( 1 . Thus, evaluating the ) b F ℓ,𝑚 ({ 𝑇 𝑘 } 𝑁 𝑘 = 1 , Q 𝑁 ) and e F ℓ,𝑚 ({ 𝑇 𝑘 } 𝑁 𝑘 = 1 , Q 𝑁 ) incurs computational cost proportional to 𝑁 log √ 𝑁 . In Algorithm 2, with 𝑀 = O GLYPH&lt;0&gt; 𝐿 2 GLYPH&lt;1&gt; evaluation points, Step 1 is computed by (19), (20), (21) and (4), (9) for degree up to 𝐿 + 1 has computational cost O GLYPH&lt;16&gt; √ 𝑀 GLYPH&lt;17&gt; ; Step 2 uses 8 times adjoint FFTs for scalar spherical harmonics and then has computational cost O GLYPH&lt;16&gt; 𝑀 log √ 𝑀 GLYPH&lt;17&gt; ; Step 3 computes (22) by the results of the previous steps has cost O ( 1 . Thus, the ) FaVeST for the adjoint case has computational cost proportional to 𝑀 log √ 𝑀 . These analyses show that the computational complexity of the proposed algorithms is nearly linear. We then call the algorithms Fa st Ve ctor S pherical Harmonic T ransform or FaVeST , and call Algorithms 1 and 2 forward FaVeST and adjoint FaVeST respectively.

4.3.2 Errors. Let 𝑇 be a tangent field in L 2 ( S 2 ) . The approximation error of FwdVSHT b F ℓ,𝑚 ({ 𝑇 ( x 𝑘 )} 𝑁 𝑘 = 1 , Q 𝑁 ) and e F ℓ,𝑚 ({ 𝑇 ( x 𝑘 )} 𝑁 𝑘 = 1 , Q 𝑁 ) for the divergence-free and curl-free coefficients b 𝑇 ℓ,𝑚 and e 𝑇 ℓ,𝑚 in (11) depends on the approximation quality of quadrature rule {( 𝑤 , 𝑘 x 𝑘 )} 𝑁 𝑘 = 1 for integrals on the sphere and the smoothness of the tangent field 𝑇 . Given a tangent field, the choice of quadrature rule is the key to reducing the approximation error. A quadrature rule {( 𝑤 , 𝑖 x 𝑖 )} 𝑁 𝑖 = 1 is called exact for polynomials of degree 𝐿 𝐿 , ≥ 0 if for all spherical polynomials 𝑝 of degree at most 𝐿 ,

$$\sum _ { i = 1 } ^ { N } w _ { i } p ( { \mathbf x } _ { i } ) = \underset { \Phi ^ { 2 } } { \int } p ( { \mathbf x } ) d \sigma ( { \mathbf x } ),$$

see e.g. [29]. Algorithm 1 using a quadrature rule that is exact for scalar spherical polynomial of degree 𝐿 has the approximation error 𝐶𝐿 -𝑠 for y ∗ ℓ,𝑚 𝑇 lying in Sobolev space 𝐻 𝑠 ( S 2 ) (of scalar spherical functions) and 𝑠 &gt; 1, where the constant 𝐶 depends only on the Sobolev norm of the function y ∗ ℓ,𝑚 𝑇 . The order 𝐿 -𝑠 is optimal, as a consequence of [26-28]. On the other hand, Algorithm 2 can evaluate FwdVSHT with zero-loss. Its approximation for the expansion (2) is equal to the truncation error for the vector spherical harmonic expansion, as determined by

$$\left \| \sum _ { \ell = L + 1 } ^ { \infty } \sum _ { m = - \ell } ^ { \ell } \left ( a _ { \ell, m } y _ { \ell, m } + b _ { \ell, m } z _ { \ell, m } \right ) \right \| _ { L _ { 2 } } ^ { 2 } = \sum _ { \ell = L + 1 } ^ { \infty } \sum _ { m = - \ell } ^ { \ell } \left ( | a _ { \ell, m } | ^ { 2 } + | b _ { \ell, m } | ^ { 2 } \right ),$$

where we have used the orthogonality of vector spherical harmonics.

## 4.4 Software Description

We provide the software package in Matlab for FaVeST which includes Matlab demo and routines for FaVeST and the routines for numerical examples in the next section. The FaVeST package can be downloaded from GitHub at https://github.com/mingli-ai/FaVeST. It has been tested in Matlab environment in operating systems including Ubuntu 16.04.6, macOS High Sierra, macOS Mojave, Windows 7, Windows 8 and Windows 10. In the Matlab library repository, the main m-files include FaVeST\_fwd.m and FaVeST\_adj.m , corresponding to Algorithms 1 and 2 respectively. Inside these two functions, the package NFFT 1 [36] is used to run the scalar FFTs on the sphere. (The NFFT is the only package requested by the FaVeST package.) There are four inputs for the function FaVeST\_fwd.m : 𝑇, 𝐿, 𝑋, 𝑤 , where 𝑇 is a tangent field sampled at the point set 𝑋 𝐿 , is the highest degree of spherical harmonics, 𝑋 is the set of quadrature nodes and 𝑤 is the set of quadrature weights. The output includes the two sequences of divergence-free and curl-free coefficients b F ℓ,𝑚 ( 𝑇 · ) and e F ℓ,𝑚 ( 𝑇 · ) for degree ℓ = 1 , . . . , 𝐿 and 𝑚 = -ℓ, . . . , ℓ as given in (10). The input of FaVeST\_adj.m contains two finite sequences 𝑎 ℓ,𝑚 and 𝑏 ℓ,𝑚 (with ℓ ≤ 𝐿 for some 𝐿 ≥ 1 and 𝑚 = -ℓ, . . . , ℓ ) as the coefficients for divergence-free and curl-free parts, and a point set 𝑋 for evaluation. The output is the AdjVSHT in (16) for 𝑎 ℓ,𝑚 and 𝑏 ℓ,𝑚 .

## 5 NUMERICAL EXAMPLES

In this section, we show numerical examples for verification of the performance of the proposed FaVeST algorithm. We start from the description of two types of polynomial-exact quadrature rules used in the experiments, and then present the examples of tangent fields on the sphere. We show the reconstruction and error in orthographic projection on the sphere for FaVeST using both kinds of point sets. We also show the CPU computational time of FaVeST which evaluates a tangent field for degree up to 2 250 and at 10 million spherical points. ,

## 5.1 Quadrature Rules

We use two types of point sets on S 2 in the experiments, as follows.

- (1) Gauss-Legendre tensor product rule (GL) [30]. The Gauss-Legendre tensor (product) rule is a (polynomial-exact but not equal area) quadrature rule Q 𝑁 : = {( 𝑤 , 𝑖 x 𝑖 )} 𝑁 𝑖 = 1 , 𝑖 = 0 , . . . , 𝑁 } on the sphere generated by the tensor product of the Gauss-Legendre nodes on the interval [-1 1 , ] and equi-spaced nodes on the longitude with non-equal weights. To be exact for polynomials of degree 𝐿 , one needs around 2 𝐿 2 GL points. Figure 1(a) shows 𝑁 = 512 GL points for degree 𝐿 = 16. An alternative GL was developed by [21, 22, 34] with smaller number of points.

1 https://www-user.tu-chemnitz.de/~potts/nfft

- (2) Symmetric spherical 𝑡 -designs (SD) [71]. The symmetric spherical 𝑡 -design is a (polynomial-exact) quadrature rule Q 𝑁 : = {( 𝑤 , 𝑖 x 𝑖 )} 𝑁 𝑖 = 1 , 𝑖 = 0 , . . . , 𝑁 } on the sphere S 2 with equal weights 𝑤 𝑖 = 1 / 𝑁 . The points are almost uniformly (or in formal definition, quasi-uniform 2 ) distributed on the sphere. To be exact for polynomials of degree 𝐿 , one needs to use the symmetric spherical 𝑡 -design for 𝑡 = 𝐿 with around 𝐿 2 / 2 points. Figure 1(b) shows 𝑁 = 498 SD points for degree 𝐿 = 31.

Fig. 1. Point sets S 2 . (a) Nodes of Gauss-Legendre tensor rule (GL). (b) Nodes of symmetric spherical designs (SD).

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-e65bf61b-0f40-443e-a0b0-ac1e0acebdef.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-49d64774-a6d6-4952-91e6-1c06be5120e4.jpg)

## 5.2 Tangent Fields

To verify our theoretical results in Section 3, we use three types of simulated tangent fields as provided in [16]. All these tangent fields are generated using stream function and velocity potential so that we can easily split the divergence-free and curl-free parts of the field. Let 𝑠 and 𝑣 be the stream function and velocity potential, then, each of the tangent fields can be represented by

$$T = \underbrace { L s } _ { f ^ { d i v } } + \underbrace { \nabla _ { * } v } _ { f ^ { c u r l } }.$$

Recall from Section 3 that L and ∇∗ denote the surface curl and surface gradient, and L 𝑠 and ∇∗ 𝑣 are divergence-free and curl-free. We define the three tangent fields, as follows.

Tangent Field A. The stream function and velocity potential for this field are linear combinations of spherical harmonics. They can be used to generate realistic synoptic-scale meteorological wind fields. The stream function is defined by √

$$s _ { 1 } ( { \boldsymbol x } ) = - \frac { 1 } { \sqrt { 3 } } Y _ { 1, 0 } ( { \boldsymbol x } ) + \frac { 8 \sqrt { 2 } } { 3 \sqrt { 3 8 5 } } Y _ { 5, 4 } ( { \boldsymbol x } ),$$

which is known as a Rosby-Haurwitz wave and is an analytic solution to the nonlinear barotropic vorticity equation on the sphere [33, pp. 453-454]. In [70], 𝑠 1 was used as the initial condition for the shallow water wave equations on the sphere. The velocity potential is given by

$$v _ { 1 } ( { \mathbf x } ) = \frac { 1 } { 2 5 } ( Y _ { 4, 0 } ( { \mathbf x } ) + Y _ { 6, - 3 } ( { \mathbf x } ) ).$$

Note that we can choose spherical harmonics with different degrees and different coefficients in the above formula. Here, we have used the same setting as [16].

2 The quasi-uniformity describes the distribution of a sequence of point sets. A sequence of point sets on the sphere is called quasi-uniform if the quotient of the covering and separation radii of each point set is bounded, see e.g. [45, 52, 71].

Tangent Field B. This field still uses the Rosby-Haurwitz wave (24) as the stream function. But the velocity potential is a linear combination of compactly supported functions:

$$v _ { 2 } ( x ) = \frac { 1 } { 8 } f ( x ; 5, \pi / 6, 0 ) - \frac { 1 } { 7 } f ( x ; 3, \pi / 5, - \pi / 7 ) + \frac { 1 } { 9 } f ( x ; 5, - \pi / 6, \pi / 2 ) - \frac { 1 } { 8 } f ( x ; 3, - \pi / 5, \pi / 3 ),$$

where

$$f ( { \mathbf x } ; \sigma, \theta _ { c }, \lambda _ { c } ) = \frac { \sigma ^ { 3 } } { 1 2 } \sum _ { j = 0 } ^ { 4 } ( - 1 ) ^ { j } \left ( \begin{array} { c } 4 \\ j \end{array} \right ) \left | r - \frac { ( j - 2 ) } { \sigma } \right | ^ { 3 }.$$

Tangent Field C . Let x 𝑐 ∈ 𝑆 2 in spherical coordinates ( 𝜃 𝑐 , 𝜆 𝑐 ) , and 𝑡 = x · x 𝑐 and 𝑎 = 1 -𝑡 . Define

$$g ( x ; \theta _ { c }, \lambda _ { c } ) = - \frac { 1 } { 2 } ( ( 3 t + 3 \sqrt { 2 } a ^ { 3 / 2 } - 4 ) + ( 3 t ^ { 2 } - 4 t + 1 ) \log ( a ) + ( 3 t - 1 ) a \log ( \sqrt { 2 a } + a ) ).$$

The stream function for this tangent field is given by

$$s _ { 3 } ( x ) = \int _ { - \pi / 2 } ^ { \theta } \sin ^ { 1 4 } ( 2 \xi ) d \xi - 3 g ( x ; \pi / 4, - \pi / 1 2 ),$$

where 𝜃 denotes the latitudinal coordinate of x . With 𝑔 , the velocity potential is given by

$$v _ { 3 } ( { \boldsymbol x } ) = \frac { 5 } { 2 } g ( { \boldsymbol x } ; \pi / 4, 0 ) - \frac { 7 } { 4 } g ( { \boldsymbol x } ; \pi / 6, \pi / 9 ) - \frac { 3 } { 2 } g ( { \boldsymbol x } ; 5 \pi / 1 6, \pi / 1 0 ).$$

## 5.3 Reconstruction and Errors

The left columns of Figures 2 and 3 present the tangent fields sampled at 𝑁 = 1922 GL points and 𝑁 = 1894 SD points, respectively. The middle columns of Figure 2 and 3 show the reconstructed field 𝑇 rec for the tangent field 𝑇 with 𝑁 = 1922 GL points and 𝑁 = 1849 SD points for evaluation. The corresponding point-wise error 𝑇 -𝑇 rec (at the sampling points) is displayed in the right columns. Here, the length and arrow indicate the scalar value and direction of the field at a spherical point. We observe that the reconstruction performance of the FaVeST is excellent: the relative error for the reconstruction is small compared with the magnitude of the original field in both GL and SD cases. The reconstructed field restores the direction and scalar value of the original field in high precision, especially when the field is sufficiently smooth. The error fields for Tangent Fields B and C show some local features of the evaluation of FaVeST : the smoother part has a smaller error, which may be partially interpreted by the localization approximation behaviour of spherical polynomials, see, e.g. [18, 40, 66-68].

Table 1 reports the relative 𝐿 2 -errors ∥ 𝑇 -𝑇 rec ∥ 2 /∥ 𝑇 ∥ 2 for the reconstruction with degree 𝐿 up to 150 in the cases of GL and SD. We observe that the FaVeST s with GL and SD exhibit similar error orders for each tangent field. These results illustrate that the FaVeST with a polynomial-exact quadrature rule is precise. In both GL and SD cases, for each degree, the error increases as the smoothness of the field decreases (from A to C). It means that the accuracy of the algorithm improves as the smoothness of the tangent field reduces. The error here is the superposition of the errors for forward FaVeST and adjoint FaVeST . They are determined by the approximation error of numerical integration by quadrature rule and the truncation error of vector spherical harmonic expansion, both of which increase as the smoothness of the tangent field reduces, see, e.g. [6, 7, 14].

We use the same example to show the error incurred by repeated use of FaVeST . In the experiment, we repeat the forward and adjoint transforms twice to obtain the sequence 𝑇 0 →{ b e 𝑇 ,𝑇 0 0 } → 𝑇 1 →{ b e 𝑇 ,𝑇 1 1 } → 𝑇 2 . Here 𝑇 0 is an original field on the sphere with the vector spherical harmonic coefficients { b e 𝑇 ,𝑇 0 0 } ; 𝑇 1 is the reconstructed field from the coefficients of 𝑇 0 via adjoint FaVeST . 𝑇 2 is the adjoint FaVeST for the coefficients { b e 𝑇 ,𝑇 1 1 } of the field 𝑇 1 . In Table 2, we evaluate the errors for the above transforms on the tangent fields A, B, C over GL and SD points for vector spherical

|                 | Quadrature   |     𝐿 = 10 |     𝐿 = 30 |     𝐿 = 50 |    𝐿 = 100 |    𝐿 = 120 |    𝐿 = 150 |
|-----------------|--------------|------------|------------|------------|------------|------------|------------|
| Tangent Field A | GL           | 8.6133e-12 | 4.3287e-12 | 3.1993e-12 | 2.6626e-12 | 2.5678e-12 | 2.4932e-12 |
|                 | SD           | 5.3367e-12 | 3.2721e-12 | 2.9385e-12 | 2.5873e-12 | 4.0959e-12 | 1.4713e-10 |
| Tangent Field B | GL           | 0.075102   | 0.0029206  | 0.00065037 | 0.00010389 | 7.1028e-05 | 3.5907e-05 |
|                 | SD           | 0.071155   | 0.0030262  | 0.00069277 | 0.0001128  | 7.6203e-05 | 4.0761e-05 |
| Tangent Field C | GL           | 0.26919    | 0.0057499  | 0.0019874  | 0.00047176 | 0.0003386  | 0.00025335 |
| Tangent Field C | SD           | 0.17225    | 0.0056764  | 0.0021336  | 0.00055568 | 0.00041357 | 0.00024919 |

Table 1. Relative 𝐿 2 -errors ∥ 𝑇 -𝑇 rec ∥ 2 /∥ 𝑇 ∥ 2 in GL and SD cases for various numbers of nodes.

|                 | Quadrature   | ∥ 𝑇 1 - 𝑇 0 ∥ ∞       | ∥ 𝑇 2 - 𝑇 0 ∥ ∞       | ∥ 𝑇 2 - 𝑇 1 ∥ ∞       | max GLYPH<8> ∥ b 𝑇 1 - b 𝑇 0 ∥ ∞ , ∥ e 𝑇 1 - e 𝑇 0 ∥ ∞ GLYPH<9>   |
|-----------------|--------------|-----------------------|-----------------------|-----------------------|-------------------------------------------------------------------|
| Tangent Field A | GL SD        | 2.0874e-12 4.0994e-12 | 4.0647e-12 8.2248e-12 | 2.0074e-12 4.1511e-12 | 2.9352e-12 2.9310e-13                                             |
| Tangent Field B | GL SD        | 1.4946e-03 2.0078e-03 | 1.4946e-03 2.0078e-03 | 2.8605e-12 4.5585e-12 | 2.9400e-12 3.0065e-13                                             |
| Tangent Field C | GL           | 1.9968e-02            | 1.9968e-02            |                       |                                                                   |
|                 |              |                       |                       | 1.2819e-11 5.6859e-12 | 3.0886e-12 6.1313e-13                                             |
|                 | SD           | 3.8045e-02            | 3.8045e-02            |                       |                                                                   |

Table 2. Errors when using forward and adjoint FaVeST repeatedly: 𝑇 0 → GLYPH&lt;8&gt; b e GLYPH&lt;9&gt; 𝑇 ,𝑇 0 0 → 𝑇 1 → GLYPH&lt;8&gt; b e GLYPH&lt;9&gt; 𝑇 ,𝑇 1 1 → 𝑇 2 , for the tangent fields A, B, C over GL and SD points, respectively.

harmonics up to degree 40 respectively, by computing the infimum norms ∥ 𝑇 1 -𝑇 0 ∥∞ , ∥ 𝑇 2 -𝑇 0 ∥∞ , ∥ 𝑇 2 -𝑇 1 ∥∞ , and max GLYPH&lt;8&gt; ∥ b 𝑇 1 -b 𝑇 0 ∥∞ ∥ , e 𝑇 1 -e 𝑇 0 ∥∞ GLYPH&lt;9&gt; . For all the three fields, ∥ 𝑇 1 -𝑇 0 ∥∞ and ∥ 𝑇 2 -𝑇 0 ∥∞ are at the same order, and the errors are bigger as the original tangent field becomes rougher. This is due to that they contain the projection error of the spherical harmonic approximation to the tangent field. ∥ 𝑇 2 -𝑇 1 ∥∞ is at the order 10 -12 which is the same to that of the scalar spherical harmonic transforms NFFT that we have used in FaVeST . This error indeed comes from the error of the NFFT since 𝑇 1 and 𝑇 2 are the reconstructed fields that are expanded by vector spherical harmonics at the same degree. The error for the two sets of coefficients max GLYPH&lt;8&gt; ∥ b 𝑇 1 -b 𝑇 0 ∥∞ ∥ , e 𝑇 1 -e 𝑇 0 ∥∞ GLYPH&lt;9&gt; is at the order equal or less than 10 -12 , which is the analogue to that of ∥ 𝑇 2 -𝑇 1 ∥∞ dependent upon the NFFT's precision.

## 5.4 Computational Complexity

To test the time complexity of FaVeST , we carry out experiments with different numbers of GL points: let 𝐿 = 250 𝑘 + 500, 𝑘 = 0 1 , , . . . , 8, which corresponds to 𝑁 𝑘 ≈ 2 𝐿 2 𝑘 GL nodes for FwdVSHT, and 𝑀 𝑘 = 𝐿 2 𝑘 + 𝐿 𝑘 coefficients for AdjVSHT. The CPU time consumption by FaVeST s is reported in Table 3. As indicated by the quotient ratios (in the round brackets) of the CPU times for degrees 𝐿 𝑘 and 𝐿 𝑘 -1 , the computational time is almost proportional to 𝑁 𝑘 and 𝑀 𝑘 for forward and adjoint FaVeST s. This is also confirmed by Figure 4. Figure 4(a) shows the trend of CPU time changing with the increase of the number of GL points for forward FaVeST . Figure 4(b) shows the trend of CPU time varying with the increase of coefficient number for adjoint FaVeST . The fitting curves for the forward and adjoint FaVeST cases have the powers (up to a constant) 𝑁 1 2 . and 𝑀 1 2 . , which means that both the forward and adjoint FaVeST s have near linear computational complexity. This is consistent with the analysis in Section 4.3.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-12a5cb12-de7b-43b6-9a8c-3a649c03ac2a.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-0508e629-cf19-4053-961e-9a6f213c4927.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-1865127d-5846-47a1-afdb-2bfbbe1b76f1.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-3a30d55b-e353-4f97-83e0-9578b3cb4532.jpg)

(a) Reconstruction of Tangent Field A with GL

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-578a7e61-85d6-4937-9d9d-d3b170d0a9b4.jpg)

(b) Reconstruction of Tangent Field B with GL

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-0fc551d0-f417-4d38-af5c-739cecdf0664.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-c7986322-601b-4786-8fbc-2d91145ce1fb.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-b6bc7023-deed-41ce-9581-d741c26fb42e.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-11523957-2b8b-4c6e-8d40-060929300cda.jpg)

(c) Reconstruction of Tangent Field C with GL

Fig. 2. Visualization of reconstructed tangent field and error by FaVeST with quadrature rule of Gauss-Legendre tensor (GL). The first and the second columns show the target tangent field 𝑇 and reconstructed field 𝑇 rec . The third column shows the error 𝑇 -𝑇 rec . All plots are the orthographic projections of the fields evaluated at 𝑁 = 1922 GL nodes. The normalized maximum norms for 𝑇 𝑇 , rec and 𝑇 -𝑇 rec are displayed in the titles.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-5be65e16-b297-43b8-b83b-60469df68fc0.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-c802b98e-df7f-407b-8575-c353d8ebacf3.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-e9d14441-6fa7-48bb-ae3d-6b32d1681767.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-cca071b3-a9f9-43a6-9558-fd341e9ab48b.jpg)

(a) Reconstruction of Tangent Field A with SD

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-b7408eec-85e3-497f-9afb-efe5196c11b7.jpg)

(b) Reconstruction of Tangent Field B with SD

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-230e028b-a6b3-4cd7-91a7-82aa6dfc1c54.jpg)

(c) Reconstruction of Tangent Field C with SD

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-4abd6b99-ad7a-4eb7-b00b-6d9c60a6a5ef.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-46c84068-88dc-4b34-bd22-aa2e2b182149.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-855d56d3-9413-4443-96c6-19e02794c926.jpg)

Fig. 3. Visualization of reconstructed tangent field and error by FaVeST with quadrature rule of symmetric spherical 𝑡 -designs (SD). The first and the second columns show the target tangent field 𝑇 and reconstructed field 𝑇 rec . The third column shows the error 𝑇 -𝑇 rec . All plots are the orthographic projections of the fields evaluated at 𝑁 = 1 849 , SD nodes. The normalized maximum norms for 𝑇 𝑇 , rec and 𝑇 -𝑇 rec are displayed in the titles.

| 𝐿     | 250        | 500        | 750                  | 1000         | 1250        | 1500        | 1750        | 2000        | 2250             |
|-------|------------|------------|----------------------|--------------|-------------|-------------|-------------|-------------|------------------|
| 𝑁     | 126,002    | 502,002    | 1,128,002 4.22(1.73) | 2,004,002    | 3,130,002   | 4,506,002   | 6,132,002   | 8,008,002   | 10,134,002 62.74 |
| 𝑡 fwd | 0.40(4.56) | 1.81(2.33) |                      | 7.2903(2.08) | 15.16(1.47) | 22.30(1.47) | 32.87(1.26) | 41.51(1.51) |                  |
| 𝑀     | 63,000     | 251,000    | 564,000              | 1,002,000    | 1,565,000   | 2,253,000   | 3,066,000   | 4,004,000   | 5,067,000        |
| 𝑡 adj | 0.83(4.64) | 3.87(2.46) | 9.52(1.88)           | 17.89(1.90)  | 33.97(1.48) | 50.32(1.47) | 73.94(1.32) | 97.32(1.49) | 144.70           |

Table 3. Forward FaVeST CPU time 𝑡 fwd v.s. number of points and adjoint FaVeST CPU time 𝑡 adj v.s. number of coefficients. For 𝐿 = 𝐿 𝑘 = 250 𝑘 + 500 , 𝑘 = 0 1 , , . . . , 8 , forward FaVeST uses Gauss-Legendre tensor rule which has 𝑁 = 𝑁 𝑘 ≈ 2 𝐿 2 𝑘 nodes, and adjoint FaVeST uses 𝑀 = 𝑀 𝑘 = 𝐿 2 𝑘 + 𝐿 𝑘 coefficients. The numbers inside the brackets are the ratios 𝑡 fwd ( 𝑁𝑘 ) 𝑡 fwd ( 𝑁𝑘 -1 ) and 𝑡 adj ( 𝑀𝑘 ) 𝑡 adj ( 𝑀𝑘 -1 ) . The numerical test is run under Intel Core i7-6700 CPU @ 3.40GHz with 16GB RAM in Windows 10.

Fig. 4. CPU time of Forward and Adjoint FaVeST s.

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-52f16d85-8172-40a7-a20f-21993b4bd129.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-e364185c-3095-4cdc-a3fd-2a00adad7f04.jpg)

## 5.5 Numerical Stability

We test the numerical stability of FaVeST . The stability can be measured by the max ratio of the discrete Fourier transform (10) and the sum of absolute values of individual terms in (15). See for example [31]. For example, for forward FaVeST , the ratio is simplified as

$$\widehat { r } _ { L, N } = \max _ { \substack { \tau \leq L, | m | = 1, \dots, \tau \\ k = 1, \dots, N } } \frac { | y ^ { * } _ { \ell, m } ( { \mathbf x } _ { k } ) | } { \widehat { V } _ { \ell, m, k } } \quad \widetilde { r } _ { L, N } = \max _ { \substack { \tau \leq L, | m | = 1, \dots, \tau \\ k = 1, \dots, N } } \frac { | y ^ { * } _ { \ell, m } ( { \mathbf x } _ { k } ) | } { \widetilde { V } _ { \ell, m, k } },$$

where

$$\text{where} \\ \widehat { V } _ { \ell, m, k } & = \sqrt { 2 } \left ( \xi _ { \ell - 1, m - 1 } ^ { ( 1 ) } | Y _ { \ell - 1, m - 1 } ^ { * } ( x _ { k } ) | + \xi _ { \ell + 1, m - 1 } ^ { ( 2 ) } | Y _ { \ell + 1, m - 1 } ^ { * } ( x _ { k } ) | + \xi _ { \ell - 1, m + 1 } ^ { ( 3 ) } | Y _ { \ell - 1, m + 1 } ^ { * } ( x _ { k } ) | + \xi _ { \ell + 1, m + 1 } ^ { ( 4 ) } | Y _ { \ell + 1, m + 1 } ^ { * } ( x _ { k } ) | } \\ & \quad + \xi _ { \ell - 1, m } ^ { ( 5 ) } | Y _ { \ell - 1, m } ^ { * } ( x _ { k } ) | + \xi _ { \ell + 1, m } ^ { ( 6 ) } | Y _ { \ell + 1, m } ^ { * } ( x _ { k } ) | \right ) \\ \widetilde { V } _ { \ell, m, k } & = \sqrt { 2 } \left ( \mu _ { \ell, m - 1 } ^ { ( 1 ) } | Y _ { \ell, m - 1 } ^ { * } ( x _ { k } ) | + \mu _ { \ell, m + 1 } ^ { ( 3 ) } | Y _ { \ell, m + 1 } ^ { * } ( x _ { k } ) | + \mu _ { \ell, m } ^ { ( 2 ) } | Y _ { \ell, m } ^ { * } ( x _ { k } ) | \right ).$$

$$\widetilde { V } _ { \ell, m, k } = \sqrt { 2 } \left ( \mu _ { \ell, m - 1 } ^ { ( 1 ) } | Y _ { \ell, m - 1 } ^ { * } ( { \mathbf x } _ { k } ) | + \mu _ { \ell, m + 1 } ^ { ( 3 ) } | Y _ { \ell, m + 1 } ^ { * } ( { \mathbf x } _ { k } ) | + \mu _ { \ell, m } ^ { ( 2 ) } | Y _ { \ell, m } ^ { * } ( { \mathbf x } _ { k } ) | \right ).$$

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-e847d907-9b97-478e-923e-2fc556039418.jpg)

(a) b 𝑟 5 ,𝑁

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-bec2ec31-8bf2-4873-a783-c7758acdaf68.jpg)

(b) e 𝑟 5 ,𝑁

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-4586a271-3564-4b10-a59a-249ac942227c.jpg)

b

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-7a9b9528-0693-4f07-aa0a-f49af3a2893f.jpg)

e

Fig. 5. Stability test for forward FaVeST . The top row shows that the quantities b 𝑟 5 ,𝑁 and e 𝑟 5 ,𝑁 have an envelope with positive power rate, which indicates their increasing trend as number of points 𝑁 increases. The bottom row shows that the quantities b 𝑟 5 ,𝑁 / 𝑁 and e 𝑟 5 ,𝑁 / 𝑁 both have an envelope with negative power rate, which indicates their decreasing trend as 𝑁 increases.

The top row of Figure 5 shows the b 𝑟 5 ,𝑁 and e 𝑟 5 ,𝑁 for 𝑁 from 100 up to 10,000. The errors in two pictures in the top panel both have an envelope with an increasing trend as the number of points 𝑁 increases. The bottom row of Figure 5 shows the b 𝑟 5 ,𝑁 / 𝑁 and e 𝑟 5 ,𝑁 / 𝑁 for 𝑁 from 100 up to 10,000. The errors in two pictures in the bottom panel both have an envelope with a decreasing trend as 𝑁 increases. Similar trending results with respect to 𝑁 are illustrated in Figure 6, where the top row shows the b 𝑟 10 ,𝑁 and e 𝑟 10 ,𝑁 , and the bottom row shows the the b 𝑟 10 ,𝑁 / 𝑁 and e 𝑟 10 ,𝑁 / 𝑁 , for 𝑁 from 100 up to 100,000. They illustrate that forward FaVeST has good stability. Computing similar quantities as (25) we can show the adjoint FaVeST also has good stability.

## 6 CONCLUSIONS AND DISCUSSION

This work proposes the first concrete fast algorithm which evaluates the forward and adjoint transforms of vector spherical harmonics for tangent fields. The fast algorithm (which we call FaVeST ) is made possible from the representation of FwdVSHT and AdjVSHT by scalar spherical harmonics and Clebsch-Gordan coefficients. By scalar FFTs on the sphere,

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-d5a58bbe-caed-4223-a172-f56e79805cb8.jpg)

(a) b 𝑟 10 ,𝑁

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-d7eaa552-16b8-46a6-9a9d-0e6f67299023.jpg)

(b) e 𝑟 10 ,𝑁

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-c110d31d-dd83-46d2-8b22-aedd87ffce24.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/1908.00041v3-figure-e93270ad-a222-4630-859b-5bad10d1d9c9.jpg)

b

e

Fig. 6. Stability test for forward FaVeST . The top row shows that the quantities b 𝑟 10 ,𝑁 and e 𝑟 10 ,𝑁 have an envelope with positive power rate, which indicates their increasing trend as number of points 𝑁 increases. The bottom row shows that the quantities b 𝑟 10 ,𝑁 / 𝑁 and e 𝑟 10 ,𝑁 / 𝑁 both have an envelope with negative power rate, which indicates their decreasing trend as 𝑁 increases.

the proposed algorithms for FwdVSHT and AdjVSHT both achieve the near-linear computational complexity. The accuracy, computational speed and numerical stability of FaVeST are validated by the numerical examples of simulated tangent fields. We develop a software package in Matlab for FaVeST , which works for polynomial-exact quadrature rules. This package has been used in the fast computation for tensor needlet transform - a fast multiresolution analysis -for tangent field on the sphere in [41]. In future work, we will use the package for solving partial differential equations on the sphere such as Stokes or Navier-Stokes equations on S 2 . As full skype maps in cosmology are usually evaluated at HEALPix points [20], we will develop a software package of FaVeST for HEALPix in Python in future.

## ACKNOWLEDGMENTS

The authors thank the anonymous referees' helpful comments and E. J. Fuselier and G. B. Wright for providing their MATLAB program which generates simulated tangent fields. The authors also thank P. Broadbridge and A. Olenko for their helpful comments during the initial stage of this work. M. Li acknowledges support from the National Natural

Science Foundation of China under Grant 61802132, and the Australian Research Council under Discovery Project DP160101366 when he worked with P. Broadbridge and A. Olenko at La Trobe University. Q. T. Le Gia and Y. G. Wang acknowledges support from the Australian Research Council under Discovery Project DP180100506.

## REFERENCES

- [1] P. A. Ade, N. Aghanim, M. Arnaud, M. Ashdown, J. Aumont, C. Baccigalupi, A. Banday, R. Barreiro, J. Bartlett, N. Bartolo, et al. Planck 2015 results-XV. Gravitational lensing. Astronomy &amp; Astrophysics , 594:A15, 2016.
- [2] N. Aghanim, Y. Akrami, F. Arroja, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. Barreiro, N. Bartolo, et al. Planck 2018 results-I. Overview and the cosmological legacy of Planck. Astronomy &amp; Astrophysics , 641:A1, 2020.

[3] N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, M. Ballardini, A. J. Banday, R. Barreiro, N. Bartolo, S. Basak, et al. Planck 2018 results-VIII. Gravitational lensing. Astronomy &amp; Astrophysics , 641:A8, 2020.

[4] A. H. Barnett, J. Magland, and L. af Klinteberg. A parallel nonuniform fast fourier transform library based on an 'exponential of semicircle' kernel. SIAM Journal on Scientific Computing , 41(5):C479-C504, 2019.

[5] R. G. Barrera, G. A. Estevez, and J. Giraldo. Vector spherical harmonics and their application to magnetostatics. European Journal of Physics , 6(4):287-294, oct 1985.

[6] J. Brauchart, J. Dick, E. Saff, I. Sloan, Y. Wang, and R. Womersley. Covering of spheres by spherical caps and worst-case error for equal weight cubature in Sobolev spaces. Journal of Mathematical Analysis and Applications , 431(2):782-811, 2015.

[7] 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. Mathematics of Computation , 83(290):2821-2851, 2014.

[8] E. O. Brigham and E. O. Brigham. The Fast Fourier Transform and its Applications , volume 448. Prentice Hall Englewood Cliffs, NJ, 1988.

[9] F. Dai and Y. Xu. Approximation Theory and Harmonic Analysis on Spheres and Balls . Springer Monographs in Mathematics. Springer, New York, 2013.

[10]

NIST Digital Library of Mathematical Functions. http://dlmf.nist.gov, Release 1.0.9 of 2014-08-29, 2014. Online companion to [

?

].

[11] J. B. Drake, P. Worley, and E. D Azevedo. Algorithm 888: Spherical harmonic transform algorithms. ′ ACM Transactions on Mathematical Software , 35(3):23, 2008.

[12] A. R. Edmonds. Angular Momentum in Quantum Mechanics . Princeton University Press, 2016.

[13] M. Fan, D. Paul, T. C. M. Lee, and T. Matsuo. Modeling tangential vector fields on a sphere. Journal of the American Statistical Association , 113(524):1625-1636, 2018.

[14] W. Freeden, T. Gervens, and M. Schreiner. Constructive Approximation on the Sphere . Numerical Mathematics and Scientific Computation. The Clarendon Press, Oxford University Press, New York, 1998.

[15] W. Freeden and M. Schreiner. Spherical Functions of Mathematical Geosciences: A Scalar, Vectorial, and Tensorial Setup . Springer-Verlag, 2009.

[16] E. J. Fuselier and G. B. Wright. Stability and error estimates for vector field interpolation and decomposition on the sphere with RBFs. SIAM Journal on Numerical Analysis , 47(5):3213-3239, 2009.

[17] M. Ganesh, Q. T. Le Gia, and I. H. Sloan. A pseudospectral quadrature method for Navier-Stokes equations on rotating spheres. Mathematics of Computation , 80(275):1397-1430, 2011.

[18]

Q. T. L. Gia, I. H. Sloan, Y. G. Wang, and R. S. Womersley. Needlet approximation for isotropic random fields on the sphere.

Journal of Approximation

Theory

,

216:86-116, 2017.

[19] F. Giraldo and T. Rosmond. A scalable spectral element Eulerian atmospheric model (SEE-AM) for NWP: Dynamical core tests. Monthly Weather Review , 132:133-153, 2004.

[20] K. M. Gorski, E. Hivon, A. J. Banday, B. D. Wandelt, F. K. Hansen, M. Reinecke, and M. Bartelmann. HEALPix: A framework for high-resolution discretization and fast analysis of data distributed on the sphere. The Astrophysical Journal , 622(2):759-771, apr 2005.

[21] D.-M. M. Gräf. Efficient algorithms for the computation of optimal quadrature points on Riemannian manifolds . PhD thesis, Universitätsverlag Chemnitz, 2013.

[22] M. Gräf and D. Potts. On the computation of spherical designs by a new optimization approach based on fast spherical fourier transforms. Numerische Mathematik , 119(4):699-724, 2011.

[23] E. W. Grafarend. Three-dimensional deformation analysis: Global vector spherical harmonic and local finite element representation. Tectonophysics , 130(1):337-359, 1986.

[24] D. Healy, P. J. Kostelec, and D. Rockmore. Towards safe and effective high-order Legendre transforms with applications to FFTs for the 2-sphere. Advances in Computational Mathematics , 21(1-2):59-105, 2004.

[25] D. M. Healy, D. N. Rockmore, P. J. Kostelec, and S. Moore. FFTs for the 2-sphere-improvements and variations. Journal of Fourier Analysis and Applications , 9(4):341-385, 2003.

[26] K. Hesse and I. H. Sloan. Optimal lower bounds for cubature error on the sphere 𝑆 2 . Journal of Complexity , 21(6):790-803, 2005.

[27] K. Hesse and I. H. Sloan. Worst-case errors in a Sobolev space setting for cubature over the sphere 𝑆 2 . Bulletin of the Australian Mathematical Society , 71(1):81-105, 2005.

- [28] K. Hesse and I. H. Sloan. Cubature over the sphere 𝑆 2 in Sobolev spaces of arbitrary order. Journal of Approximation Theory , 141(2):118-133, 2006.
- [29] K. Hesse, I. H. Sloan, and R. S. Womersley. Numerical Integration on the Sphere , pages 2671-2710. Springer Berlin Heidelberg, Berlin, Heidelberg, 2015.
- [30] K. Hesse and R. S. Womersley. Numerical integration with polynomial exactness over a spherical cap. Advances in Computational Mathematics , 36(3):451-483, 2012.
- [31] N. J. Higham. Accuracy and Stability of Numerical Algorithms . SIAM, 2002.
- [32] E. L. Hill. The theory of vector spherical harmonics. American Journal of Physics , 22(4):211-214, 1954.
- [33] J. R. Holton. An introduction to dynamic meteorology. American Journal of Physics , 41(5):752-754, 1973.
- [34] J. Keiner, S. Kunis, and D. Potts. NFFT 3.5, C subroutine library. Contributors: F. Bartel, M. Fenn, T. Gorner, M. Kircheis, T. Knopp, M. Quellmalz, M. Schmischke, T. Volkmer, A. Vollrath.
- [35] J. Keiner, S. Kunis, and D. Potts. Efficient reconstruction of functions on the sphere from scattered data. Journal of Fourier Analysis and Applications , 13(4):435-458, Aug 2007.
- [36] J. Keiner, S. Kunis, and D. Potts. Using NFFT 3 - a software library for various nonequispaced fast Fourier transforms. ACM Transactions on Mathematical Software , 36(4):19, 2009.
- [37] J. Keiner and D. Potts. Fast evaluation of quadrature formulae on the sphere. Mathematics of Computation , 77(261):397-419, 2008.
- [38] S. Kim. The complete superconformal index for 𝑁 = 6 Chern-Simons theory. Nuclear Physics B ,
- 821(1):241-284, 2009.
- [39] S. Kunis and D. Potts. Fast spherical Fourier algorithms. Journal of Computational and Applied Mathematics , 161(1):75-98, 2003.
- [40] Q. T. Le Gia and H. N. Mhaskar. Localized linear polynomial operators and quadrature formulas on the sphere. SIAM Journal on Numerical Analysis , 47(1):440-466, 2008/09.
- [41] M. Li, P. Broadbridge, A. Olenko, and Y. G. Wang. Fast tensor needlet transforms for tangent vector fields on the sphere. arXiv preprint arXiv:1907.13339 , 2019.
- [42] F. J. Lowes. Mean-square values on sphere of spherical harmonic vector fields. Journal of Geophysical Research , 71(8):2179-2179, 1966.
- [43] F.-s. Lu, J.-q. Song, W.-q. Lin, Y.-f. Pang, K.-j. Ren, and P.-c. Shi. Efficient utilization of launched threads on GPUs: The spherical harmonic transform as a case study. Computer Physics Communications , 184(11):2494-2502, 2013.
- [44] D. Marinucci and G. Peccati. Random Fields on the Sphere , volume 389 of London Mathematical Society Lecture Note Series . Cambridge University Press, Cambridge, 2011.
- [45] H. N. Mhaskar, F. J. Narcowich, and J. D. Ward. Spherical Marcinkiewicz-Zygmund inequalities and positive quadrature. Mathematics of Computation , 70(235):1113-1130, 2001.
- [46] M. J. Mohlenkamp. A fast transform for spherical harmonics. Journal of Fourier Analysis and Applications , 5(2-3):159-184, 1999.
- [47] J. C. N'edélec. Acoustic and Electromagnetic Equations , volume 144 of Applied Mathematical Sciences . Springer-Verlag, New York, 2001.
- [48] D. Pekurovsky. P3DFFT: A framework for parallel computations of Fourier transforms in three dimensions. SIAM Journal on Scientific Computing , 34(4):C192-C209, 2012.
- [49] W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. Numerical Recipes 3rd Edition: The Art of Scientific Computing . Cambridge University Press, 2007.
- [50] Reinecke, M. and Seljebotn, D. S. Libsharp - spherical harmonic transforms revisited. Astronomy &amp; Astrophysics , 554:A112, 2013.
- [51] V. Rokhlin and M. Tygert. Fast algorithms for spherical harmonic expansions. SIAM Journal on Scientific Computing , 27(6):1903-1928, 2006.
- [52] M. Sourisseau, Y. G. Wang, H.-T. Wu, and W.-H. Yu. Optimization-based quasi-uniform spherical t-design and generalized multitaper for complex physiological time series. arXiv preprint arXiv:1907.13493 , 2019.
- [53] R. Suda. Stability analysis of the fast Legendre transform algorithm based on the fast multipole method. In Proceedings of the Estonian Academy of Sciences , volume 53, pages 107-115. Estonian Academy Publishers; 1999, 2004.
- [54] R. Suda. Fast spherical harmonic transform routine FLTSS applied to the shallow water test set. Monthly Weather Review , 133(3):634-648, 2005.
- [55] R. Suda and M. Takami. A fast spherical harmonics transform algorithm. Mathematics of Computation , 71(238):703-715, 2002.
- [56] P. N. Swarztrauber. Spectral transform methods for solving the shallow-water equations on the sphere. Monthly Weather Review , 124(4):730-744, 1996.

[57] P. N. Swarztrauber. Shallow water flow on the sphere. Monthly Weather Review , 132(12):3010-3018, 2004.

[58] P. N. Swarztrauber and W. F. Spotz. Generalized discrete spherical harmonic transforms. Journal of Computational Physics , 159(2):213-230, 2000.

[59] P. T. P. Tang. DFTI - a new interface for Fast Fourier Transform libraries. ACM Transactions on Mathematical Software , 31(4):475-507, 2005.

[60] M. Tygert. Fast algorithms for spherical harmonic expansions, II. Journal of Computational Physics , 227(8):4260-4279, 2008.

[61] M. Tygert. Fast algorithms for spherical harmonic expansions, III. Journal of Computational Physics , 229(18):6181-6192, 2010.

[62] D. A. Varshalovich, A. N. Moskalev, and V. K. Khersonskii. Quantum Theory of Angular Momentum . World Scientific, 1988.

[63] S. K. Veerapaneni, D. Gueyffier, G. Biros, and D. Zorin. A numerical method for simulating the dynamics of 3D axisymmetric vesicles suspended in viscous flows. Journal of Computational Physics , 228(19):7233 - 7249, 2009.

[64] S. K. Veerapaneni, A. Rahimian, G. Biros, and D. Zorin. A fast algorithm for simulating vesicle flows in three dimensions. Journal of Computational Physics , 230(14):5610 - 5634, 2011.

[65] B. Wang, L.-L. Wang, and Z. Xie. Accurate calculation of spherical and vector spherical harmonic expansions via spectral element grids. Advances in Computational Mathematics , 44(3):951-985, Jun 2018.

- [66] Y. Wang. Filtered polynomial approximation on the sphere. Bulletin of the Australian Mathematical Society , 93(1):162-163, 2016.
- [67] Y. G. Wang, Q. T. L. Gia, I. H. Sloan, and R. S. Womersley. Fully discrete needlet approximation on the sphere. Applied and Computational Harmonic Analysis , 43(2):292-316, 2017.

[68] Y. G. Wang, I. H. Sloan, and R. S. Womersley. Riemann localisation on the sphere. Journal of Fourier Analysis and Applications , 24(1):141-183, Feb 2018.

[69] E. J. Weinberg. Monopole vector spherical harmonics. Physical Review D , 49:1086-1092, Jan 1994.

[70] D. L. Williamson, J. B. Drake, J. J. Hack, R. Jakob, and P. N. Swarztrauber. A standard test set for numerical approximations to the shallow water equations in spherical geometry. Journal of Computational Physics , 102(1):211-224, 1992.

[71] R. S. Womersley. Efficient spherical designs with good geometric properties. In Contemporary Computational Mathematics - A Celebration of the 80th Birthday of Ian Sloan , pages 1243-1285. Springer. URL: https://web.maths.unsw.edu.au/~rsw/Sphere/EffSphDes/, 2018.
Quoc T. Le Gia Ming Li Yu Guang Wang
FaVeST: Fast Vector Spherical Harmonic Transforms
Vector spherical harmonics on the unit sphere of $\mathbb{R}^3$ have broad
applications in geophysics, quantum mechanics and astrophysics. In the
representation of a tangent vector field, one needs to evaluate the expansion
and the Fourier coefficients of vector spherical harmonics. In this paper, we
develop fast algorithms (FaVeST) for vector spherical harmonic transforms on
these evaluations. The forward FaVeST evaluates the Fourier coefficients and
has a computational cost proportional to $N\log \sqrt{N}$ for $N$ number of
evaluation points. The adjoint FaVeST which evaluates a linear combination of
vector spherical harmonics with a degree up to $\sqrt{M}$ for $M$ evaluation
points has cost proportional to $M\log\sqrt{M}$. Numerical examples of
simulated tangent fields illustrate the accuracy, efficiency and stability of
FaVeST.