Spherical harmonic coefficients of isotropic polynomial functions with applications to gravity field modeling

Various aspects of gravity field modeling rely upon analytical mathematical functions for calculating spherical harmonic coefficients. Such functions allow quick and efficient evaluation of cumbersome convolution integrals defined on the sphere. In this work, we present a new analytical method for determining spherical harmonic coefficients of isotropic polynomial functions. This method in computationally flexible and efficient, since it makes use of recurrence relations. Also, its use is universal and could be extended to piecewise polynomials and polynomials with compact support. Our numerical investigation of the proposed method shows that certain recurrence relations lose accuracy as the order of implemented polynomials increases because of accumulation of numerical errors. Propagation of these errors could be mitigated by hybrid methods or using extended precision arithmetic. We demonstrate the relevance of our method in gravity field modeling and discuss two areas of application. The first one is the design of B-spline windows and filter kernels for the low-pass filtering of gravity field functionals (e.g., GRACE Follow-On monthly geopotential solutions). The second one is the calculation of spherical harmonic coefficients of isotropic polynomial covariance functions.


Introduction
Global-scale geophysical signals defined on the Earth's surface can be expanded in a series of spherical harmonic functions, assuming a spherical approximation for the Earth's figure.This type of representation is typical for the Earth's of sampled data exist (Colombo 1981;Sneeuw 1994;Blais 2011a, b;Rexer and Hirt 2015;Plonka et al. 2018;Sun 2021).The evaluation of convolution integrals in the spherical harmonic domain usually depends on whether methods for the calculation of spherical harmonic coefficients of integral kernels are available.The development of mathematical expressions that perform this task in an accurate and efficient manner is therefore of great importance.
A representative example of an integral kernel in physical geodesy is the Stokes kernel, which is used for the gravimetric determination of geoidal undulations.In other words, geoidal undulations can be derived via the spherical convolution of gravity anomalies with the Stokes kernel.Several other Stokes-like kernels exist that form corresponding integral relations among different functionals of the Earth's disturbing potential.An extensive list of such kernels is provided in Hirt et al. (2011).Mathematical expressions for the calculation of their spherical harmonic coefficients have already been developed, and their derivation is usually based on skillful algebraic manipulation of the generating function of Legendre polynomials (Hotine 1969).
In many practical applications, geophysical signals are observed in a limited region of the Earth's surface.The evaluation of convolution integrals for local-scale studies is done by truncating the kernel function to a radius that encloses the area of available observations (Vanńček et al. 2003).The corresponding truncation error needs to be calculated in these cases, in order to estimate the effect of neglecting the areas outside the evaluation radius.Methods for the calculation of spherical harmonic coefficients of truncated kernels (i.e., compactly supported kernels) are therefore also needed.This issue has been extensively studied for the Stokes and Vening-Meinesz kernels (de Witte 1967;Hagiwara 1972;Paul 1973;Hagiwara 1976;Paul 1983).
Filtering of a geophysical signal defined on a spherical surface is also a process that can be described with a convolution integral, denoting the convolution of the spatial signal with the corresponding filter kernel.In this case, the filter kernel is also defined by an analytical function.Methods for the calculation of spherical harmonic coefficients of filter kernels are therefore required to perform filtering in the spherical harmonic domain.Mathematical expressions in the form of recurrence relations for the spherical harmonic coefficients of the isotropic Pellinen, Gaussian, truncated Gaussian and Hanning filter kernels were derived by Jekeli (1981).
Filtering methods in the spherical harmonic domain have gained an increased interest since the development of the first Gravity Recovery and Climate Experiment (GRACE) and GRACE Follow-On (GRACE-FO) temporal geopotential solutions.The isotropic Gaussian filter has also been used extensively in GRACE-related applications.In addition to physical geodesy, the task of calculating the spherical harmonic coefficients of analytical functions is also relevant in other branches of geosciences, such as geostatistics, to assess the validity (i.e., positive definiteness) and examine the properties of covariance functions (Huang et al. 2011;Lantuéjoul et al. 2019).
In this paper, we develop expressions for the calculation of spherical harmonic coefficients of monomial and polynomial functions.We are mostly interested in isotropic, i.e., rotationally symmetric, functions that are still used primarily in gravity field modeling.In addition to purely polynomial functions, our expressions can be practical for the calculation of spherical harmonic coefficients of the polynomial part of already existing kernels.
The use of polynomial expressions is common in geostatistics, and specifically in the design of polynomial covariance functions, to model the stochastic characteristics of a random field.Apart from their mathematical simplicity, another property of polynomial functions that makes them appealing for such applications is their infinite differentiability.To the best of our knowledge, there is a lack of supporting literature on the use of polynomial kernels in gravity field modeling, and the calculation of their spherical harmonic coefficients.The only relevant study appears to be the paper of Ma (2016), where expressions for the Gegenbauer coefficients of isotropic polynomial functions are derived.These expressions are then used to characterize covariance matrices whose entries are polynomials of order up to four.Since the Gegenbauer coefficients are a generalization of spherical harmonic coefficients for higher-dimensional spheres (i.e., spheres of dimension higher than two), the expressions of Ma (2016) are also a generalization of our expressions.
The novelty of our work lies on the fact that the expressions derived here can be used for the calculation of spherical harmonic coefficients of compactly supported and piecewise polynomials, whereas the study of Ma (2016) examines the case of polynomials defined on the entire sphere.We further discuss the relevance of our methods in gravity field modeling and provide examples and applications that are motivated from the discussion of the present section.
The paper is organized as follows.In Sect. 2 we provide a brief mathematical overview of the representation of isotropic functions in the spherical harmonic domain by means of the discrete Legendre transform.We also derive a generalization, i.e., the compactly supported discrete Legendre transform, that permits the spherical harmonic representation of piecewise functions and functions of compact support.In Sect. 3 we develop a recurrent algorithm for the calculation of spherical harmonic coefficients of isotropic monomials and extend our mathematical formulation to isotropic polynomials.In Sect. 4 we design a numerical experiment to assess the relative accuracy of our recurrent algorithm.Sections 5 and 6 discuss two potential applications of our algorithm in gravity field modeling.In Sect.5, the concept of isotropic B-spline filtering on the sphere is introduced, where filter kernels are described by piecewise polynomial functions.Analytical expressions for the representation of B-spline filter kernels in the spatial and spherical harmonic domains are also provided.The third-order B-spline filter is examined in detail and used for the denoising of a GRACE-FO monthly gravity field solution.In Sect.6, we review some isotropic polynomial covariance models on the sphere and provide analytical expressions for the calculation of their spherical harmonic coefficients that are based on our algorithm.Finally, Sect.7 provides a summary of this work and highlights the main conclusions.

Spherical harmonic representation of isotropic functions
Two-point functions are typically used in geodesy to model the spatial covariance of a random field, and as kernel functions of integral transforms to describe the relation between the integration and calculation points (Devaraju and Sneeuw 2018).A two-point function f (φ 1 , λ 1 , φ 2 , λ 2 ) defined on a spherical surface is isotropic, i.e., invariant under rotation of the reference system, when it only depends on the spheri- where ψ is given by and with φ ∈ [−π/2, π/2] and λ ∈ [0, 2π) being the geocentric latitude and longitude, respectively.The dimensionless spherical harmonic coefficients F n of an isotropic function f are given by the discrete Legendre transform J , which is defined as (Jekeli 2017, p. 54) with P n denoting the Legendre polynomials of degree n ∈ N 0 .The isotropic function f can be reconstructed from the sequence of coefficients F n using the inverse Legendre transform J −1 as follows: In practice, the sum of Eq. ( 4) is evaluated up to a maximum degree n max , for which the coefficients F n are available.As in Fourier analysis, a greater n max results in a better reconstruction of f .
The operational calculus and properties of the discrete Legendre transform have been studied by mathematicians since the early fifties (Tranter 1950;Churchill 1954), and several generalizations have been developed since then.For n ∈ R, the discrete Legendre transform can be generalized to the continuous Legendre transform (Butzer et al. 1980).As for the kernel function, the Gegenbauer and Jacobi transforms (Conte 1955;Scott 1953) are also considered generalizations of the discrete Legendre transform because the Gegenbauer and Jacobi polynomials are both generalizations of the Legendre polynomials.
In this section, we generalize the transform J with respect to the integration interval.We therefore define the compactly supported discrete Legendre transform J I , for which the signal is analyzed in the subinterval of spherical distance Based on the definition of Eq. ( 5), the transition from J to J I is done by constraining f to be compactly supported in I , i.e., f (ψ) = 0, ∀ψ / ∈ I .The following relation can therefore be formulated: where 1 I is the indicator function (i.e., a unit boxcar function), defined as A relation between F n and F

(I )
n is also provided in Appendix A. This relation is based on an alternative formulation of the indicator function and is derived using some well-known methods (e.g., Wieczorek and Simons 2005;Gupta and Narasimhan 2007;Dahlen and Simons 2008).The rationale for constraining f using the indicator function comes from the fact that, in the present study, f is assumed to be an analytical function (i.e., a polynomial function).The indicator function corresponds to an ideal window and ensures that f is perfectly suppressed outside and perfectly preserved inside interval I .The use of any other window w in Eq. ( 6) will introduce large distortions in both the shape and properties of f .As a consequence, the resulting spherical harmonic coefficients will correspond to the distorted analytical function w(ψ) f (ψ), and not the original intended one.In addition, the compactly supported Legendre transform as defined in Eq. ( 6) appears naturally when the spectrum of a piecewise or compactly supported analytical function is required.This is demonstrated in subsequent sections.
The spatial confinement of f by a user-specified window w (e.g., a Gaussian or any other smooth window) provides a more general definition of the compactly supported discrete Legendre transform, which is also in better agreement with the definition of the short-time Fourier transform.The choice of a smooth window, instead of an ideal window, would be a more appropriate choice for localized spectral analysis of geophysical signals on the sphere.In this type of applications, the accurate estimation of periodic constituents in the frequency spectrum is ensured by the use of a smooth window that mitigates leakage effects.A relation between F n and F (I ) n using any window w can also be derived using the method presented in Appendix A, provided that the spectrum of w can be easily calculated.For simple cases where w is centered at ψ = 0, and a Gaussian or a Hanning window is selected, the spectrum of w can be calculated using the recurrence relations given in Jekeli (1981).
Since this study assumes f to be an analytical function, and not a geophysical signal, we only consider the case w(ψ) = 1 I (ψ).It can be easily deduced that Accounting for the properties of the indicator function, it is also evident by Eq. (6) that n in the interval of support I is possible by the application of J −1 as follows: The compactly supported discrete Legendre transform can be especially useful for the analysis of piecewise functions (e.g., covariance functions, window functions and filter kernels with compact support).We do not examine the mathematical properties of the compactly supported discrete Legendre transform, as it is outside the main scope of this study.We henceforth use the term "Legendre transforms" to refer to both the discrete Legendre transform and the compactly supported discrete Legendre transform.

Spherical harmonic coefficients of isotropic polynomials
In this section, we provide expressions for the Legendre transforms of the monomial ψ m , with m ∈ N 0 , which represents the basis function of a polynomial vector space.We first treat the general case of the compactly supported discrete Legendre transform, and then corresponding expressions are provided for the discrete Legendre transform.These results are afterward used to derive expressions for the calculation of spherical harmonic coefficients of any polynomial, including piecewise polynomials, of any order m.

The compactly supported discrete Legendre transform of monomials
The compactly supported discrete Legendre transform of the m-th-order monomial ψ m , denoted as n,m , is given by The derivation of a closed-form expression for the integral of Eq. ( 9) can be a very demanding, if not an impossible, task.It is therefore customary in the geodetic practice, and also for reasons of computational efficiency, to evaluate integrals in the form of Eq. ( 9) using recurrence relations.These are usually derived by exploiting corresponding recurrence relations and other properties of the Legendre polynomials.
In the case of Eq. ( 9), we closely follow the procedure given in Jekeli (1981, Appendix B), which we also outline in Appendix B, and derive a recurrence relation for (I )  n,m that reads where and with the parameters B n,m (ψ) and D n,m (ψ) defined using the auxiliary variable k n (t Recurrence relations for the calculation of B n,m (ψ) and D n,m (ψ) are provided in Appendix C. For n = {0, 1}, the coefficients 0,m and 1,m are given by the definite integrals Equations ( 13a) and (13b) can be evaluated by the closedform expressions (Gradshteyn and Ryzhik 2014, §2.633, eq. 1) (I ) with (•) s denoting the falling factorial, i.e., (x The following recurrence relations can also be derived for 0,m and (I ) 1,m using the results in Gradshteyn and Ryzhik (2014, §2.631, eqs. 1 and 2) with initial conditions and with initial conditions The relations provided in the present section and in Appendix C allow the development of a framework for the calculation of n,m in an entirely recurrent fashion (e.g., Algorithm 1).As an example, Fig. 1a shows the spherical harmonic coefficients of the monomials ψ, ψ 2 and ψ 3 up to n max = 200 and for I = [80 • , 120 • ].The reconstructed monomials are presented in Fig. 1b, with the Gibbs effect being clearly visible at discontinuity points, e.g., endpoints of interval I , (Gelb 1997).

The discrete Legendre transform of monomials
The discrete Legendre transform of the m-th-order monomial ψ m is given by and can be calculated by setting ψ 1 = 0 and The recurrence relations of Eqs. ( 10), ( 15) and ( 17) are simplified as follows: with initial conditions Equation ( 21a) can be verified using the more general results given in Ma (2016).In his study, he provides recurrence relations for the discrete Gegenbauer transform of ψ m , defined as Algorithm 1: Calculation of for maximum harmonic degree N and maximum monomial order M for n, m = {0, 1} ; // Eqs.(16,18) {ψ 1 , ψ 2 } and n =  is the Gegenbauer polynomials.The Legendre polynomials are special cases of Gegenbauer polynomials, with The recurrence relation of Ma (2016) for n,1/2,m reads after some algebraic manipulation Substituting Eq. ( 24) into Eq.( 25) and solving with respect to n,m gives In order to show that Eqs. ( 21a) and ( 26) are equivalent, it is sufficient to prove that which is done in Appendix D, thus confirming the validity of Eq. (21a).

The Legendre transforms of polynomials
The compactly supported discrete Legendre transform of the M-th-order polynomial can be derived using the linearity property where α m , a and b are constants.Implementing the linearity property to Eq. ( 28) results in the expression The calculation of (I ) n,m is discussed in Sect.3.1 and offers a simple and elegant way of evaluating the spherical harmonic coefficients of any polynomial.The recurrent calculation of n,m is also especially convenient and efficient for polynomials with no missing monomial terms, i.e., with α m = 0, ∀m ∈ {x ∈ N 0 : x ≤ M}. Figure 2 shows the spherical harmonic coefficients and spatial reconstruction of a fifth-order polynomial s(ψ), as defined by Eq. ( 28 The results given so far in this section can be further generalized for piecewise polynomials.Let q be a piecewise function of K nonzero segments, with each segment q k to be defined in the interval and with I k being mutually disjoint, i.e., I i ∩ I j = ∅ for i = j.Equation ( 31) can also take the form Assuming that each segment q k is a polynomial of order M k with constants α m,k , Eq. ( 32) is rewritten as We note that Eqs. ( 32) and ( 33), unlike Eq. ( 31), can also be used for overlapping intervals I k .The compactly supported discrete Legendre transform of q in I is given by where 34) again indicates the importance and relevance of the mathematical developments of Sect.3.1 for the calculation of spherical harmonic coefficients of piecewise polynomials.Analogous expressions to the ones provided in the present section can be derived for the discrete Legendre transform by substituting

Numerical assessment
It is widely known that a sequence evaluated using a recurrence relation is susceptible to inaccuracies.These inaccuracies are due to the accumulation of truncation and round-off errors coming from each cycle of the recurrence procedure (Gautschi 1967).Since our method for the calculation of n,m is also based on recurrence relations, we design a numerical experiment to access its accuracy.We compare the recurrent calculation of (I ) n,m using doubleprecision floating-point number format, which is MATLAB's default numeric setup, against the recurrent calculation using extended precision floating-point number format.Extending the numerical precision (i.e., significant digits) of all computations, usually at the expense of computational efficiency, is a simple approach for reducing the propagation of truncation and round-off errors.This approach has been recently used by Bucha et al. (2019) and Piretzidis and Sideris (2019) to evaluate unstable recurrence relations of geodetic interest.We define the relative accuracy of where the subscripts "dp" and "ep" denote recurrent calculation using double precision and extended precision, respectively.All extended precision calculations are performed with 32 significant digits.In theory, the relative accuracy σ (I )  n,m depends on four variables at most, i.e., the spherical harmonic degree n, the polynomial order m and the two endpoints ψ 1 and ψ 2 of the interval I .
Extensive numerical experiments on the evaluation of σ (I ) n,m for different combinations of ψ 1 and ψ 2 revealed that, in practice, σ (I ) n,m is highly dependent on ψ 2 and not on ψ 1 .Therefore, to simplify the presentation of our results, we use the approximation σ and further assess the relative accuracy for different values of n, m and ψ 0 .The results are given in the upper panels of Fig. 3 is strongly influenced by both m and ψ 0 , with an increasing m and decreasing ψ 0 to result in a relative accuracy of larger magnitude.
The results of Fig. 3 also suggest that the relative accuracy slightly deteriorates as n increases, but to a far lesser extent than it does for varying m and ψ 0 .The larger magnitude of σ for increasing m indicates that the order-wise recurrence relations given by Eqs. ( 15) and ( 17) result in a large error growth.These errors are subsequently propagated by Eq. ( 10) to the rest of the coefficients.The calculation of (I )  0,m and (I ) 1,m using Eqs.( 14a) and (14b), instead of Eqs. ( 15) and ( 17), respectively, is also inaccurate after a certain order and does not provide any advantages.
In order to mitigate the propagation of large errors, we propose a hybrid method for the calculation of 1,m are evaluated using numerical integration, and the coefficients n,m for n > 1 are evaluated recurrently using Eq. ( 10).The numerical integration is performed here using the Gauss-Legendre quadrature method (Kythe and Schäferkotter 2005, §3.2.2).The relative accuracy of the hybrid method is again calculated using the extended precision arithmetic method as a reference.The results are provided in the middle row of Fig. 3, where an improvement of the relative accuracy of the hybrid method with respect to the double-precision recurrent method is evident.Additional tests showed that extending the degree n for which the coefficients n,m are evaluated numerically can further improve the relative accuracy of the hybrid method.
Finally, the last row of Fig. 3 presents the relative accuracy of is not systematically affected by m and ψ 0 in their examined range, whereas it shows a larger magnitude as n increases.Previous studies (e.g., Piretzidis and Sideris 2019) have demonstrated that numerical integration is not accurate for evaluating coefficients with a magnitude (in absolute value) close to or smaller than the relative accuracy level of the selected floating-point format.Evidence of such behavior is also provided in the numerical experiments of Sect.6.2.Therefore, deteriorating relative accuracy should also be expected for the calculation of (I ) n,m using numerical integration as the absolute value of (I ) n,m decreases.

Isotropic B-spline filtering on the sphere
The study of spectral-domain filtering techniques is an important area of research in the geodetic community since the development of the first GRACE-based temporal geopotential models.These models are affected by large correlated errors and their filtering prior to any geophysical interpretation is always necessary.A first possible application of the methods developed in Sect. 3 is the efficient calculation of the spherical harmonic coefficients of polynomial lowpass filters.Simple yet rigorous techniques of deriving the coefficients of these filters permit the study of their spectral characteristics.Such low-pass filters can subsequently be used to smooth functionals of the Earth's gravity field in the spherical harmonic domain.
In the present section, we introduce the concept of isotropic B-spline filtering on the sphere.It is a special case of isotropic filtering with the filter kernel to be described by B-spline curves.B-splines are continuous piecewise polynomial functions of finite support, with applications to interpolation, regression and parametric modeling of curves and surfaces.Due to their characteristics, they provide an excellent testbed for implementing the methods provided in Sect.3. B-splines have been previously used in geodetic studies for atmosphere modeling (Haji-Aghajany et al. 2020;Schmidt et al. 2020), topography/bathymetry modeling (Boergens et al. 2021;Zhang et al. 2022) and multiresolutional analysis (Mautzl et al. 2005).A review of the literature suggests that B-splines have not been systematically utilized in the context of filtering a spatial geophysical signal defined on the sphere.Therefore, apart from testing the methodology of Sect.3, we take the opportunity to also examine B-splines in this context.We first summarize the underlying mathematical theory of isotropic filtering on the sphere and its implementation in the spatial and spherical harmonic domains.We then provide relations for the calculation of B-spline windows and filter kernels.As an example, we derive the expressions for the third-order B-spline filter and demonstrate its applicability to the filtering of a GRACE-FO monthly solution.

Isotropic filtering on the sphere
The filtering of a spatial field g(φ, λ) with an isotropic filter kernel h(ψ) is expressed with the following convolution integral: where " * " denotes convolution, is the integration domain (surface of unit sphere) and d = cos φ dφ dλ is the area element of .The spherical distance ψ in Eq. ( 36) corresponds to the distance between the evaluation point (φ, λ) and the running integration point (φ , λ ).The filter kernel h(ψ) is usually a normalized window function w(ψ) (Jekeli 1981, eq. 51), i.e., where w represents the spatial average of w(ψ) on the unit sphere and is calculated by The evaluation of Eq. ( 36) is performed efficiently in the spherical harmonic domain using convolution theorems.
The set of spherical harmonic coefficients {G C n, , G S n, } of the (φ, λ)-dependent function g are given by the forward Fourier-Legendre transform J , defined as where Pn, are the fully normalized associated Legendre functions of harmonic degree n and order  (not to be confused with the polynomial order m defined in Sect.3).The inverse Fourier-Legendre transform J −1 reads For an isotropic function, it holds J ≡ J .The spherical harmonic coefficients of ḡ are derived using the convolution theorem 123 Fig. 3 Relative accuracy of for different values of ψ 0 using the recurrent approach (top), the hybrid approach (middle) and Gauss-Legendre quadrature (bottom) where H n denotes the Legendre transform of the filter kernel h(ψ).In the sequel, we discuss the use of B-spline windows for the construction of related isotropic kernels that can be used as spatial low-pass filters.

B-spline windows and filters
B-spline windows are a family of tapered window functions that can be used for truncating (i.e., space-limiting) a signal or for the design of low-pass filters.The K -th-order B-spline window function is a piecewise polynomial of order K − 1 that is obtained from the cascaded convolution of K rectangular window functions.B-spline windows were introduced and described mathematically by Toraichi et al. (1989).Adapting their closed-form expression of the unnormalized K -th-order B-spline window for the spherical case results in the function where the interval I k is given by the intersection and with ψ 0 ∈ (0, π] denoting the window length for which v(ψ) = 0, ∀ψ > ψ 0 .The window length ψ 0 can also be expressed in terms of the great-circle distance r = Rψ 0 on the surface of a sphere of radius R.
It is evident that Eq. ( 42) has the same form as Eq. ( 32).Using the binomial theorem, Eq. ( 42) can also be written in the form of Eq. ( 33) as follows: The K -th-order B-spline window function, normalized so that w(0) = 1, is then given by We note that in the work of Toraichi et al. (1989), Eq. ( 42) is multiplied by the constant K K .Since w is normalized, this constant is disregarded here as trivial.An analytical expression for the spatial average w is derived by substituting Eqs. ( 44) and ( 45) into Eq.( 38) and reads as follows: The filter kernel h(ψ) is then derived using Eq. ( 37), and the filter coefficients H n are given by the relation where it is evident by the last two equations that H 0 = 1.The first-order B-spline window corresponds to a rectangular (or spherical cap) window, the second order is a Bartlett (or triangular) window and the fourth order is a Parzen window (Toraichi et al. 1989).In the next subsection we derive expressions for the third-order B-spline window and filter kernel.

The third-order B-spline filter
Substituting K = 3 in Eq. ( 42), we derive the following expression for the unnormalized third-order B-spline window: which, upon expanding the sum, yields with Equation ( 49) consists of two nonzero segments and can be written in symbolic form as Expanding the terms 1 through 3 of Eq. ( 51), and taking into account that all indicator functions equal 1 inside their interval of support, we derive the expression with v(0) = ψ 2 0 /9.The normalized third-order B-spline window is given by 123 The spatial average w has the analytical form and the isotropic filter kernel h(ψ) is given by Eq. ( 37).Based on Eq. ( 53) and the theory given in Sect.3, the spherical harmonic coefficients H n of the filter h(ψ) can be evaluated by Figure 4 shows the spatial and spherical harmonic representation of the isotropic third-order B-spline filter kernel for different window lengths.

Application to GRACE-FO filtering
We demonstrate the use of the third-order B-spline filter in the post-processing of GRACE-FO monthly solutions.We use the Level 2 RL06 products calculated by the Center for Space Research (CSR), The University of Texas at Austin.Level 2 products correspond to monthly sets of spherical harmonic coefficients {G C n, , G S n, } of the Earth's gravity field and have been available since June 2018.Time-varying geopotential models derived from GRACE and GRACE-FO contain correlated errors that dominate the higher frequencies (i.e., higher degrees and orders) of the gravity field spectrum.Due to GRACE-FO's orbital and sampling characteristics (Peidou and Pagiatakis 2020;Pagiatakis and Peidou 2021), the spatial-domain geometry of these errors follows a characteristic north-south distribution commonly referred to as striping.The correlated errors are usually reduced using lowpass filters or specialized decorrelation techniques.
Since the results of this section are not used for geophysical interpretation, no corrections (e.g., addition of degree-one coefficients, replacement of C 2,0 coefficients, etc.) or prior decorrelation are applied to Level 2 products.We use the standard procedure of Wahr et al. (1998) and calculate filtered surface mass variations σ in terms of equivalent water height (EWH) using the equation where ρ ave is the Earth's average density, ρ w is the water density, κ n are the load-deformation coefficients of degree n and { ḠC n, , ḠS n, } are the filtered spherical harmonic coefficient variations.
The coefficients { ḠC n, , ḠS n, } are calculated in two steps.In the first step, the contribution of the GRACEonly static geopotential model ITSG-Grace2018s (Kvas et al. 2019;Mayer-Gürr et al. 2018) In the second step, the coefficients { G C n, , G S n, } are filtered using Eq. ( 41) resulting in the filtered coefficient variations { ḠC n, , ḠS n, }. Figure 5 shows the filtered maps of EWH variations, corresponding to the monthly solution of November 2019, using the third-order B-spline filter with different window lengths.The dominant striping effect is clearly evident in the unfiltered EWH variations of Fig. 5a.The use of a larger window length reduces the striping effect and produces smoother EWH variations.This behavior is typical and resembles other low-pass filters used for denoising GRACE and GRACE-FO data, such as the isotropic Gaussian filter (Jekeli 1981;Piretzidis and Sideris 2022).

Isotropic polynomial covariance functions
As a second application of the methods given in Sect.3, we discuss the calculation and study of polynomial covariance functions in the spherical harmonic domain.In the first part of the present section, we outline some aspects of covariance functions used in gravity field modeling.In the second part, we examine some isotropic polynomial covariance functions and provide expressions for the calculation of their spherical harmonic coefficients.

Covariance functions in gravity field modeling
Covariance functions are important mathematical components in geostatistics, since they provide essential information for the spatial or spatiotemporal correlation of a random field.In physical geodesy and gravity field modeling, covariance functions are primarily used to model the stochastic behavior of the Earth's disturbing potential and its related functionals, as well as other signals of geophysical interest.Covariance functions are also essential in the implementation of statistical methods (e.g., least-squares collocation) for the optimal estimation, interpolation and filtering of spatial fields.The use of such statistical tools and techniques is motivated by the fact that certain geophysical signals exhibit a random behavior in the lateral direction (Jekeli 2010).
Since the true covariance function of a spatial field is never known, an empirical covariance function is usually estimated The isotropic covariance function C(ψ) of a spatial random field g(φ, λ) defined on the spherical surface is given by (Hofmann-Wellenhof and Moritz 2006, eq.9-32) where ψ and α are the polar coordinates (spherical distance and azimuth, respectively) of point A(φ, λ) and A (φ , λ ).
The spherical harmonic representation of C(ψ), denoted as n = J {C(ψ)}, can be derived from Eq. ( 3).As already discussed, only discrete values of g(φ, λ) are usually available in practice, which allow the empirical estimation of C(ψ) in one step and the fitting of an analytical model to the estimated empirical values in a subsequent step.Although a large variety of analytical covariance models are available in the literature, the ones that should be accounted for are those that are valid (i.e., positive definite) in the examined space.The criteria that ensure the positive definiteness of a covariance model on the sphere, based on the characterization of n , are derived by Schoenberg (1942) as follows: In addition to Schoenberg's criteria, the analytical covariance model to be used should capture the properties of the examined spatial field.For example, a covariance function for the gravitational potential should be harmonically continued into free space (i.e., outside the Earth's masses).
Another important attribute of covariance functions is their support length.Compactly supported covariance functions, also known as finite covariance functions, play an important role in stochastic modeling.This type of covariance functions assumes zero auto-correlation after a specific distance and offers certain computational advantages.For example, in the context of optimal statistical estimation, the use of a compactly supported covariance function produces sparse covariance matrices that are faster to be calculated and stored, and results in a linear system of equations that is easier to be solved.Such covariance functions have been previously studied in Sansò and Schuh (1987), Arabelos and Tscherning (1998), Moreaux et al. (1999), Moreaux (2008).In the following subsection, we discuss some examples of compactly supported polynomial covariance functions and provide expressions for their spectral representation.

Expressions for spherical harmonic coefficients of some isotropic polynomial covariance models on the sphere
A review of the recent literature on a wide variety of topics (e.g., geostatistics, multivariate analysis, stochastic analysis, etc.) reveals that isotropic polynomial covariance models exist and are a subject of ongoing research (Schubert and Schuh 2023).The validity of these models was firstly examined on Euclidean spaces and later on spherical surfaces.The transition to the sphere is done by a simple substitution of the Euclidean distance by the great-circle distance (Gneiting 2013).Early investigations showed that some of these models meet the positive definiteness criteria for the twodimensional sphere S 2 = {x ∈ R 3 : x = R} of arbitrary radius R and therefore allow the design of rigorous covariance functions on such surfaces.A list of commonly studied isotropic polynomial covariance models is provided in Table 1 along with their parameter range that guarantees positive definiteness in S 2 .For more information regarding the positive definiteness of these models, the reader is referred to the work of Gneiting (2013).Some examples on the study and use of these models in gravity field modeling, in the context of optimal estimation of gravity field functionals, are given in Moreaux (2008) and Klees and Prutkin (2010).It is evident by their expressions that all models listed in Table 1 are compactly supported.The ψ 0 parameter denotes the range of the covariance function, for which C(ψ) = 0 for ψ > ψ 0 , and controls the extent of spatial correlation (Christakos 1992, p. 76); (Wackernagel 2003, p. 58).The range ψ 0 can also be expressed as a great-circle distance r = Rψ 0 on the Earth's surface, as already done for the window length in Sect.5.2.The τ parameter controls the shape of the covariance function.
The spherical covariance model is given by a third-order polynomial expression.From a geometric perspective, this expression describes the normalized volume of the intersection of two spheres with diameter ψ 0 (in linear units) and with their centers separated by a vector of magnitude ψ (also in linear units) (Wackernagel 2003, p. 60).In other words, the spherical model is derived by the normalized selfconvolution of the indicator function of a ball with diameter ψ 0 in R 3 .Using the theory of Sect.3, the spherical harmonic coefficients of the spherical model are given by . ( The self-convolution of a ball in higher-dimensional spaces gives rise to the so-called Euclid's hat functions.Gneiting (1999) showed that particular representations of the Euclid's hat functions correspond to the Askey covariance model, which was introduced by Askey (1973) and is a polynomial function of order τ .Using the binomial theorem, the Askey covariance model is written as and its spherical harmonic coefficients are calculated by One common drawback of the spherical and Askey covariance models is their non-differentiability (Wackernagel 2003, p. 58); (Schaback 2011).A family of covariance models with no such limitation is the Wendland models (Wendland 1995).These models represent polynomials of arbitrary order and smoothness.They are developed via the repeated application of the "Montée" integral operator of  Matheron (1965) to the Askey covariance model (Gneiting 2002).The most frequently used Wendland models are the ones with smoothness degree two and four, termed C 2 -and C 4 -Wendland models, and correspond to polynomial functions of order τ + 1 and τ + 2, respectively.The binomial expansion of the C 2 -and C 4 -Wendland models reads as follows and their spherical harmonic coefficients are given by the expressions respectively.Figures 6 and 7 show some examples of the covariance models of Table 1 in the spatial and spherical harmonic domain for different values of ψ 0 and τ .The all-positive values of n in Fig. 7 are expected since the covariance models of Table 1 are positive definite; hence, the Schoenberg's criteria apply.
The method for the evaluation of n,m (e.g., recurrent, hybrid, etc.) plays a crucial role in the calculation of spherical harmonic coefficients of positive definite covariance functions as the propagation of numerical errors in higher degrees can result in negative n values.For example, the calculation of n for the C 4 -Wendland model with parameter pairs {r , τ } = {600 km, 7} and {1400 km, 6} using the recurrent method yields n values with increasing oscillations up to n ≈ 2300 and 1800, respectively, and negative values afterward (Fig. 8a).For the same model and parameter pairs, the hybrid method shows a behavior similar to the recurrent method, with the oscillatory and negative n values to become evident in higher degrees (Fig. 8b).The use of Gauss-Legendre quadrature (GLQ) numerical integration leads to inaccurate results when n values are close to the relative accuracy level of the implemented floating-point format.Since the GLQ numerical integration is performed using a double-precision floating-point format, inaccurate results start to manifest when n < 2.2 × 10 −16 , which corresponds to n > 2050 for the second parameter pair, and with negative n values appearing in higher degrees (Fig. 8c).The numerical integration method is also important.For comparison, Fig. 8c also shows the resulting n values using MATLAB's function integral that implements the adaptive quadrature (ADQ) method of Shampine (2008) with a built-in control of the absolute and relative error.The relative error tolerance is chosen to be equal to 2.2×10 −16 .The ADQ numerical integration results in negative and divergent values when n ≈ 10 −12 and appears to perform worse than the GLQ for this application.Finally, the recurrent calculation of n using extended precision arithmetic results in smoothly decreasing and positive values, corroborating the superiority of the method (Fig. 8d).We also note that any deviation of the covariance model parameters from their corresponding ranges given in the third column of Table 1 results in a nonpositive definite covariance function, i.e., some n values might be negative (Fig. 9).
As a final remark, we briefly mention the applicability of the methods developed in Sect. 3 for the analysis of piecewise polynomial covariance functions.Such functions are studied by Sünkel (1978) as non-positive definite approximations of exact covariance functions.In his work, he examined simple approximations using piecewise step functions, linear func-

Summary and conclusions
Despite their simple mathematical form, functions still provide viable modeling options in applied and natural sciences, including Earth sciences.In an effort to better facilitate the use of such functions in geodetic applications, and specifically in gravity field modeling where certain tasks are performed in the spherical harmonic domain, we developed a recurrent method for the calculation of spherical harmonic coefficients of monomials and polynomials defined on the two-dimensional sphere.We limit our analysis to isotropic polynomial functions; therefore, their spherical harmonic coefficients were evaluated using the discrete Legendre transform.The transition from the discrete Legendre transform to the more generalized compactly supported discrete Legendre transform enabled the extension of our recurrent method to piecewise polynomials and polynomials with compact support.
A numerical assessment of the method was performed to evaluate the accuracy of all recurrence relations.Results indicate that the recurrence relations used for the calculation of coefficients of fixed degree n = {0, 1} and increasing polynomial order m quickly become inaccurate due to the propagation of large numerical errors.Since these relations are used for the initialization of the recurrent method, their numerical errors also propagate to the rest of the coeffi- cients.The polynomial order for which the coefficients can be evaluated accurately depends on the interval of support and the overall accuracy requirements of the application under study.The numerical errors of the recurrent method were easily reduced using extended precision arithmetic techniques.Alternatively, a hybrid method was developed for the mitigation of numerical errors, where the coefficients for n = {0, 1} were evaluated using numerical integration and the coefficients for n > 1 were evaluated recurrently using standard double precision.Studying the relative accuracy of the hybrid method showed that it allows the accurate calculation of coefficients for higher monomial orders than the double-precision recurrent method.For application where accuracy is important, we recommend the implementation of the recurrent method using extended precision arithmetic.
Based on our developments, we constructed a framework for isotropic B-spline filtering on the sphere and provided explicit expressions for the calculation of B-spline weights and filter kernels in the spatial and spherical harmonic domains.We demonstrated the use of the third-order B-spline filter for reducing the correlated errors in a GRACE-FO monthly gravity field solution.Results showed that the thirdorder B-spline filter follows the same behavior as other low-pass filters, i.e., a larger window length results in a smoother gravity field solution.Additional analysis of the spectral characteristics of B-spline filters, and comparison with other low-pass filtering techniques, could shed more light on their performance.
We finally examine some positive definite polynomial covariance models in the spatial and spherical harmonic domains.The calculation of their spherical harmonic coeffi-Fig.8 Calculation of C 4 -Wendland model coefficients using (a) recurrence method with double precision, (b) hybrid method, (c) GLQ (solid lines) and ADQ (dotted lines) numerical integration, and (d) recurrence method with extended precision.Negative coefficients are ignored cients is carried out by implementing our recurrent method to analytical expressions that were also developed in this study.The first Schoenberg criterion, stating that a positive definite covariance function has nonnegative coefficients, was experimentally verified for all models.A decay of the coefficients toward zero is also observed for all covariance models.This decay can be attributed to the second Schoenberg criterion that guarantees the convergence of the infinite series of coefficients of positive definite functions.Our numerical investigations showed that the method of calculating spherical harmonic coefficients of positive definite covariance functions is extremely important, especially for higher degrees, since the propagation of even small numerical errors can result in negative coefficients.The recommended approach for calculating spherical harmonic coefficients of polynomial functions for such applications is again the recurrent method using extended precision arithmetic.
Using the orthogonality relation of the Legendre polynomials in the right-hand side of Eq. (A.6) and substituting W j with the expression in Eq. (A.3), we finally derive the following relation between F n and F (I )  n : where the last term in the parenthesis denotes the Wigner 3-j symbol that can be evaluated recurrently (Varshalovich et al. 1988).The k n parameter is also calculated recurrently using the relation (Piretzidis and Sideris 2021, eq. 32) with initial values k 0 (cos ψ) = 1 − cos ψ and k 1 (cos ψ) = 3 sin 2 ψ/2.Although Eq. (A.7) is a rigorous way of calculating n from F n , there are some practical limitations.For example, the double summation is evaluated up to the maximum degrees i max and j max for which F i and k j are available.This truncation results in the inaccurate calculation of F (I ) n in the vicinity of n max .From a computational viewpoint, the double summation, the calculation of the Wigner-3j symbol and the need for prior calculation of F n are all issues that add to the overall complexity of the method.It is therefore advisable that Eq. (A.7) be used in cases where it is difficult to derive an analytical formula or a recurrence relation for F (I ) n (e.g., in case of a non-homogeneous, non-isotropic function f ).A different relation between F n and F (I )  n can also be derived using the recent approach of Djellouli et al. (2021).
(I )n,m as the relative difference

Fig. 2
Fig. 2 (a) Absolute magnitude of spherical harmonic coefficients S (I ) n of a fifth-order polynomial s(ψ) as function of degree n.The evaluation is performed in the interval I = [52 • , 120 • ].(b) Spatial representation

Fig. 4
Fig. 4 (a) Kernel function h(ψ) as function of distance on the Earth's surface Rψ (R = 6378136.3m) and (b) absolute magnitude of spherical harmonic coefficients H n as function of degree n for the isotropic third-order B-spline filter