Analytical Solutions for Gravitational Potential up to Its Third-order Derivatives of a Tesseroid, Spherical Zonal Band, and Spherical Shell

The spherical shell and spherical zonal band are two elemental geometries that are often used as benchmarks for gravity field modeling. When applying the spherical shell and spherical zonal band discretized into tesseroids, the errors may be reduced or cancelled for the superposition of the tesseroids due to the spherical symmetry of the spherical shell and spherical zonal band. In previous studies, this superposition error elimination effect (SEEE) of the spherical shell and spherical zonal band has not been taken seriously, and it needs to be investigated carefully. In this contribution, the analytical formulas of the signal of derivatives of the gravitational potential up to third order (e.g., V, Vz\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{z}$$\end{document}, Vzz\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{zz}$$\end{document}, Vxx\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{xx}$$\end{document}, Vyy\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{yy}$$\end{document}, Vzzz\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{zzz}$$\end{document}, Vxxz\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{xxz}$$\end{document}, and Vyyz\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{yyz}$$\end{document}) of a tesseroid are derived when the computation point is situated on the polar axis. In comparison with prior research, simpler analytical expressions of the gravitational effects of a spherical zonal band are derived from these novel expressions of a tesseroid. In the numerical experiments, the relative errors of the gravitational effects of the individual tesseroid are compared to those of the spherical zonal band and spherical shell not only with different 3D Gauss–Legendre quadrature orders ranging from (1,1,1) to (7,7,7) but also with different grid sizes (i.e., 5∘×5∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$5^{\circ }\times 5^{\circ }$$\end{document}, 2∘×2∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2^{\circ }\times 2^{\circ }$$\end{document}, 1∘×1∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1^{\circ }\times 1^{\circ }$$\end{document}, 30′×30′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$30^{\prime }\times 30^{\prime }$$\end{document}, and 15′×15′\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$15^{\prime }\times 15^{\prime }$$\end{document}) at a satellite altitude of 260 km. Numerical results reveal that the SEEE does not occur for the gravitational components V, Vz\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{z}$$\end{document}, Vzz\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{zz}$$\end{document}, and Vzzz\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{zzz}$$\end{document} of a spherical zonal band discretized into tesseroids. The SEEE can be found for the Vxx\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{xx}$$\end{document} and Vyy\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{yy}$$\end{document}, whereas the superposition error effect exists for the Vxxz\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{xxz}$$\end{document} and Vyyz\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V_{yyz}$$\end{document} of a spherical zonal band discretized into tesseroids on the overall average. In most instances, the SEEE occurs for a spherical shell discretized into tesseroids. In summary, numerical experiments demonstrate the existence of the SEEE of a spherical zonal band and a spherical shell, and the analytical solutions for a tesseroid can benefit the investigation of the SEEE. The single tesseroid benchmark can be proposed in comparison to the spherical shell and spherical zonal band benchmarks in gravity field modeling based on these new analytical formulas of a tesseroid.

Among several schemes for a spherical shell benchmark, a discretization using tesseroids has been a common strategy to investigate the numerical performance of the algorithms to compute the tesseroid's gravitational and magnetic effects in the gravity and magnetic fields. A whole spherical shell can be discretized into tesseroids along the coordinate lines, i.e., in longitude, latitude (or colatitude), and radial directions. The relative errors between the analytical gravitational or magnetic effects and the sum of the calculated gravitational or magnetic effects of the discretized tesseroids forming the entire spherical shell are investigated by using different numerical approaches, e.g., Taylor series expansion with different orders (Grombein et al. 2013;Shen and Deng 2016;Deng and Shen 2018a, b), Taylor series expansion and kernel matrix equivalence method (Zeng et al. 2022), 3D or 2D Gauss-Legendre quadrature (GLQ) (Deng and Shen 2018b;Zhong et al. 2019;Deng et al. 2020Deng et al. , 2022, 3D or 2D adaptive discretization stack-based GLQ approach (Uieda et al. 2016;Baykiev et al. 2016;Deng and Shen 2019;Lin and Denker 2019), conditional split and double exponential quadrature method (Fukushima 2018), 2D density-based adaptive discretization GLQ method (Soler et al. 2019), variable-order GLQ approach (Qiu and Chen 2020), combination of Taylor series expansion and GLQ method (Lin et al. 2020), and combination of GLQ method with 2D adaptive discretization and discrete convolution algorithm based on 1D fast Fourier transform (Ouyang et al. 2022). Note that most of the above-mentioned studies focused on global errors of a whole spherical shell rather than the local errors of the tesseroid for this commonly used strategy.
Moreover, the spherical zonal band has often been used in the literature as a geometrical element to test and validate (or to benchmark) gravity field modeling algorithms. In general, a spherical zonal band is the difference between two concentric spherical caps with different spherical distances i and i+1 (Deng 2022), which is a volumetric element. A spherical zonal band is a component of a spherical shell, i.e., it becomes a spherical shell when its first spherical latitude equals zero and its second spherical latitude equals . When the computation point is positioned on the polar axis, the analytical expressions for a spherical zonal band can be derived for the gravitational effects, e.g., GP (Heck and Seitz 2007;Yang et al. 2022), GV (Heck and Seitz 2007;Lin and Denker 2019;Deng 2022), GGT (Lin et al. 2020;Deng 2022), and GC (Deng 2022).
Analogously, a spherical zonal band can be discretized into tesseroids with respect to the longitude and radial coordinates. In other words, a tesseroid can be treated as a sector of a spherical zonal band. Recently, it was numerically confirmed that the computation efficiency of this strategy was higher and its computation errors were smaller than those of a spherical shell discretized using tesseroid (Deng 2022). The main reason is that the shell consists of many zonal bands. This strategy can also be applied to investigate the performance of different numerical methods for calculating the gravitational effects of the tesseroid, e.g., Taylor series expansion (Heck and Seitz 2007;Wild-Pfeiffer 2008;Lin and Denker 2019;Yang et al. 2022) and GLQ (Wild-Pfeiffer 2008; Lin and Denker 2019; Deng 2022) approaches. Generally, in previous numerical experiments, the global errors of the total spherical zonal band were also evaluated rather than the local errors of the tesseroid.
The strategies of a spherical zonal band and spherical shell discretized into tesseroids effectively calculate the global errors of the gravitational effects of the total spherical zonal band and spherical shell. The numerical results from the two strategies have constraints when employing the tesseroid with different numerical methods in practical applications. In particular, the numerical results of the error analysis from the aforementioned different numerical approaches by using the strategies of the spherical zonal band and spherical shell discretized into tesseroids can only be applied to global situations in the discretized tesseroids forming the entire spherical zonal band and spherical shell (Heck and Seitz 2007;Wild-Pfeiffer 2008;Grombein et al. 2013;Shen and Deng 2016;Uieda et al. 2016;Baykiev et al. 2016;Fukushima 2018;Shen 2018a, b, 2019;Zhong et al. 2019;Lin and Denker 2019;Soler et al. 2019;Lin et al. 2020;Qiu and Chen 2020;Deng et al. 2020Deng et al. , 2022Zeng et al. 2022;Yang et al. 2022;Deng 2022;Ouyang et al. 2022). However, for local practical calculations, the numerical results from these global error assessments do not apply. For example, in Section 'Evaluation of the accuracy' of Uieda et al. (2016), the spherical shell was adopted to evaluate the computation accuracy of the GP, GV, and GGT of the discretized tesseroids that comprise the entire spherical shell. The threshold error 0.1% was only for the global error assessments of the whole spherical shell. When using the tesseroids with the algorithm in Uieda et al. (2016) for local practical calculations, the conclusion derived from this global threshold error 0.1% cannot be applied. It should be re-evaluated for the specific local error assessments. Usually, they tend to underestimate error levels due to the superposition error elimination effect (SEEE) (Deng et al. 2022). The definition of the SEEE is the possibility that the symmetry of the global spherical shell with the discretized tesseroids might cancel out some errors from individual tesseroids (Deng et al. 2022). This SEEE was not noticed in the previous works. Careful investigation of the SEEE will help to make better use of tesseroids in the practical applications of gravity field modeling. In reality, the spherical zonal band may also have the SEEE. The SEEE is merely a hypothesis proposed in Deng et al. (2022) that has not yet been validated by numerical experiments. Concerning the spherical zonal band and spherical shell discretized into tesseroids, it is necessary to analyze the SEEE by comparing the global errors of the spherical zonal band and spherical shell with the local errors of the individual tesseroid. Consequently, it is vital to examine the variation in the local errors of the individual tesseroid in detail.
The challenge of evaluating the local errors of an individual tesseroid is to derive analytical expressions for the gravitational effects of a tesseroid. This is difficult because since the concept of the tesseroid was proposed in Anderson (1976), subsequent studies have used different numerical methods to calculate the gravitational effects of a tesseroid (Kuhn 2003;Heck and Seitz 2007;Asgharzadeh et al. 2007;Wild-Pfeiffer 2008;Tsoulis et al. 2009;Grombein et al. 2013;Shen and Deng 2016;Uieda et al. 2016;Deng and Shen 2018b;Zhao et al. 2019;Soler et al. 2019;Lin et al. 2020;Zeng et al. 2022;Ouyang et al. 2022). The primary reason is that in the general formulas of the gravitational effects of a tesseroid, the double integrals of the longitude and colatitude (or latitude) coordinates were considered incapable of analytical integration in previous studies. However, some research progress has brought a turning point to this challenge and provided the research idea for the solution of this difficult problem. For example, the gravitational effects of a spherical zonal band can be solved analytically if the computation point is located on the polar axis (Heck and Seitz 2007;Marotta and Barzaghi 2017;Lin and Denker 2019;Deng 2022). Marotta and Barzaghi (2017) derived the analytical gravitational acceleration of a tesseroid when the computation point is located on the polar axis based on Sect. 5.2 of Turcotte and Schubert (2002), see also Lin and Li (2022). Thus, the key idea is the special condition of the computation point located on the polar axis. This paper goes beyond the previous studies in that the analytical formulas of the gravitational effects including the GP, GV, GGT, and GC of a tesseroid are derived when the computation point is located on the polar axis. The formulas of the gravitational effects of a spherical zonal band are also analytically derived. These new formulas of the gravitational effects of a spherical zonal band are simpler in forms and more intuitive in methods than those in Deng (2022). Moreover, the relative errors of the gravitational effects of a single tesseroid are obtained to investigate the SEEE of the spherical zonal band and spherical shell discretized into tesseroids.
The paper is organized as follows. Section 2 offers the theoretical aspects of this paper, in which the analytical formulas of the GP, GV, GGT, and GC of a tesseroid are derived in Sects. 2.1, 2.2, 2.3, and 2.4, respectively. Sections 2.5 and 2.6 derive the analytical formulas of these gravitational effects of a spherical zonal band and spherical shell, respectively. Numerical experiments are performed in Sect. 3, in which Sect. 3.1 presents the basic setup. The influences of the 3D Gauss-Legendre quadrature orders and grid sizes on the gravitational effects of the tesseroid, spherical zonal band, and spherical shell are investigated in Sects. 3.2 and 3.3, respectively. Section 3.4 compares the gravitational effects of every individual tesseroid forming the whole spherical zonal band and spherical shell. Finally, the main conclusions and outlook for future studies are summarized in Sect. 4.

Analytical Formula of the GP of a Tesseroid
A tesseroid (see Fig. 1) is the mass body in the spherical or ellipsoidal coordinate system defined by three pairs of longitudes, colatitudes (or latitudes), and radial distances (Anderson 1976). In previous studies (Kuhn 2003;Heck and Seitz 2007;Asgharzadeh et al. 2007;Wild-Pfeiffer 2008;Grombein et al. 2013;Shen and Deng 2016;Uieda et al. 2016;Shen 2018a, b, 2019), the gravitational effects of a tesseroid were often calculated numerically because the integrals with respect to the longitude and colatitude (or latitude) in their expressions are not considered to be integrated analytically. In this section, the analytical formula of the GP of a tesseroid is derived when the computation point is located on the polar axis.
Generally, the formula of the GP (V) of a tesseroid is expressed in the spherical coordinate system as (Anderson 1976;Blakely 1995;Kuhn 2003;Heck and Seitz 2007;Asgharzadeh et al. 2007;Wild-Pfeiffer 2008;Shen and Deng 2016;Uieda et al. 2016;Deng and Shen 2019): Illustration of a tesseroid in the spherical coordinate system revised from Fig. 3.8 of Anderson (1976). Δ = 2 − 1 , Δ = 2 − 1 (or Δ = 2 − 1 ), and Δr = r 2 − r 1 are the longitude, colatitude (or latitude), and radial dimensions of the tesseroid, in which Δ = Δ . P( , , r) (or P( , , r) ) and Q( � , � , r � ) (or Q( � , � , r � ) ) are the computation and integration points, respectively. There is the relation between the colatitude (or ′ ) and latitude (or ′ ) of the computation point P (or integration point Q) as + = 90 • (or � + � = 90 • ). P (x, y, z) is in the local north-oriented frame system, where x, y, and z point to the north, east, and radial up directions, respectively 1 3 where G is the gravitational constant, is the constant density of a homogeneous tesseroid, and l is the Euclidean distance between the computation point P( , , r) and integration point Q( � , � , r � ) with their longitudes ( and ′ ), colatitudes ( and ′ ), and geocentric distances (r and r ′ ). [ 1 , 2 ], [ 1 , 2 ], and [ r 1 , r 2 ] are the longitude, colatitude, and radial boundaries of a tesseroid, i.e., the integration boundaries of the tesseroid are � = [ 1 , 2 ] , � = [ 1 , 2 ] , and r � = [r 1 , r 2 ]. When the computation point P is located on the polar axis, the colatitude of the computation point P is = 0 or = . To simplify the derivation process, = 0 is adopted as the example for the situation of the computation point located on the polar axis. The derivation progress with = is similar to that of = 0 and there are small changes in the solution for = . Under this circumstance, the expression for the GP of a tesseroid with the triple integral form in Eq. (1) has the analytical solution. The detailed derivation is presented in Appendix A. The analytical formula of the GP of a tesseroid when the computation point is located on the positive part of the polar axis (i.e., = 0) is obtained as: where the expressions for the A 22 , A 12 , A 21 , A 11 , B 2 , B 1 , C 22 , C 12 , C 21 , C 11 , l 22 , l 12 , l 21 , and l 11 are listed in Table 8 in Appendix A. Note that the presented formula in Eq. (4) is singular for the r = 0 , see the more general solution in Karcol (2021).

Analytical Formula of the GV of a Tesseroid
Differentiating Eq. (1) with respect to the geocentric distance r of the computation point P, the formula of the radial GV of a tesseroid is presented as (Heck and Seitz 2007;Asgharzadeh et al. 2007;Uieda et al. 2016;Deng and Shen 2019): Similarly, when the computation point is located on the polar axis, there exists the analytical solution for Eq. (5). The analytical expression for the GV (i.e., V r , V x , and V y ) of a tesseroid is derived in Appendix B. The analytical solution for the radial GV of a tesseroid with the colatitude of the computation point = 0 is obtained as: (3) cos = cos cos � + sin sin � cos( � − ) + 3B 2 ln C 12 + l 12 C 22 + l 22 + 3B 1 ln C 21 + l 21 C 11 + l 11 where the expressions for the B 2 , B 1 , C 22 , C 12 , C 21 , C 11 , l 22 , l 12 , l 21 , and l 11 are referred to in Table 8 and the expressions for the D 22 , D 12 , D 21 , and D 11 are listed in Table 9 in Appendix B.
Note that Marotta and Barzaghi (2017) derived the analytical formula of the gravitational acceleration Δg P in Appendix 2 when the computation point is located on the polar axis based on Sect. 5.2 of Turcotte and Schubert (2002). The differences between the V r in Eq. (7) in this paper and Δg P in Eq. (2) of Marotta and Barzaghi (2017) lie in that: (1) the usage of the r 2 , r 1 , 2 , and 1 in this paper corresponds to that of the R + h , R, ′ 2 , and ′ 1 in Marotta and Barzaghi (2017); and (2) the Δg P adds the minus symbol (i.e., −), whereas the V r has no minus symbol. Despite these differences, the formula of the V r in Eq. (7) of this paper is the same as that of Δg P in Marotta and Barzaghi (2017) after the parameter exchange and absolute operation.

Analytical Formulas of the GGT of a Tesseroid
Differentiating Eq. (1) with respect to the geocentric distance r of the computation point twice, the formula of the radial-radial GGT of a tesseroid is obtained as (Asgharzadeh et al. 2007;Grombein et al. 2013;Uieda et al. 2016;Deng and Shen 2019): Analogously, Eq. (8) has the analytical solution when the computation point is located on the polar axis. The derivation of the analytical solutions for the GGT of a tesseroid is presented in Appendix C. The analytical expression for the radial-radial GGT of a tesseroid with = 0 is obtained as: where the expressions for the B 2 , B 1 , C 22 , C 12 , C 21 , C 11 , l 22 , l 12 , l 21 , and l 11 are referred to in Table 8 and the expressions for the E 22 , E 12 , E 21 , and E 11 are listed in Table 10 in Appendix C. Due to the common usage in the previous works (Uieda et al. 2016;Lin and Denker 2019;Lin et al. 2020), V rr is denoted as V zz in the following part.
Regarding the other components of the GGT, the components V xx and V yy are chosen to derive their analytical expressions when the computation point is located on the polar axis, because these three components V xx , V yy , and V zz of the GGT satisfy Laplace's equation as V xx + V yy + V zz = 0 in theory. Laplace's equation will be applied to confirm the correctness of the derived analytical expressions for the V xx , V yy , and V zz of a tesseroid.
The general formulas of the V xx and V yy of a tesseroid in Cartesian integral kernels are expressed as (Grombein et al. 2013;Uieda et al. 2016;Deng and Shen 2019): + 6B 2 ln C 22 + l 22 C 12 + l 12 + 6B 1 ln C 11 + l 11 C 21 + l 21 When the computation point is located on the polar axis, Eqs. (10) and (12) have the analytical solutions. The analytical expressions for the V xx and V yy are derived in Appendix C. The analytical expressions for the V xx and V yy with = 0 are obtained as: where the expressions for the C 22 , C 12 , C 21 , C 11 , l 22 , l 12 , l 21 , and l 11 are referred to in Table 8. The expressions for the F 2 , F 1 , H 2 , H 1 , I 22 , I 12 , I 21 , I 11 , J 2 , J 1 , K 22 , K 12 , K 21 , and K 11 are listed in Table 11 in Appendix C. The expressions for the L 2 , L 1 , M 22 , M 12 , M 21 , M 11 , N 22 , N 12 , N 21 , and N 11 are presented in Table 12 in Appendix C. The differences of the sign between F 2 and L 2 ; F 1 and L 1 ; I 22 and M 22 ; I 12 and M 12 ; I 21 and M 21 ; I 11 and M 11 ; K 22 and N 22 ; K 12 and N 12 ; K 21 and N 21 ; K 11 and N 11 should be noted.

Analytical Formulas of the GC of a Tesseroid
Differentiating Eq. (1) with respect to the geocentric distance r of the computation point at three times, the formula of the radial-radial-radial GC of a tesseroid is obtained as Shen 2018a, b, 2019): Similarly, the analytical solution exists for Eq. (16) when the computation point is located on the polar axis. The derivation of the analytical solution for the GC of a tesseroid is (11) Δ x = r � sin cos � − cos sin � cos( � − ) The analytical formula of the radial-radial-radial GC of a tesseroid with = 0 is obtained as: where the expressions for the O 22 , O 12 , O 21 , and O 11 are listed in Table 13 in Appendix D, and the expressions for the l 22 , l 12 , l 21 and l 11 are referred to in Table 8. Similarly, V rrr is represented as V zzz due to the common utilization in Deng and Shen (2018a, 2018b. Analogously, these three components V xxz , V yyz , and V zzz of the GC satisfy Laplace's equation as V xxz + V yyz + V zzz = 0 . Thus, the analytical expressions for the V xxz and V yyz of a tesseroid are derived when the computation point is located on the polar axis. Generally, the expressions for the V xxz and V yyz of a tesseroid in Cartesian integral kernels are presented as Shen 2018b, 2019): When the colatitude of the computation point is = 0 , there exist the analytical solutions for Eqs. (18) and (19). The detailed derivation process of the analytical expressions for the V xxz and V yyz of a tesseroid is presented in Appendix D. The analytical expressions for the V xxz and V yyz with = 0 are obtained as: where the expressions for the l 22 , l 12 , l 21 , and l 11 are referred to in Table 8, the expressions for the P 22 , P 12 , P 21 , P 11 , Q 22 , Q 12 , Q 21 , Q 11 , R 22 , R 12 , R 21 , R 11 , S 22 , S 12 , S 21 , S 11 , T 2 , and T 1 Table 14 in Appendix D, and the expressions for the V 22 , V 12 , V 21 , V 11 , W 22 , W 12 , W 21 , and W 11 are presented in Table 15 in Appendix D. Note that there are the differences of the sign between P 22 and V 22 ; P 12 and V 12 ; P 21 and V 21 ; P 11 and V 11 ; R 22 and W 22 ; R 12 and W 12 ; R 21 and W 21 ; R 11 and W 11 . Regarding these expressions in Eqs. (4), (7), (9), (14), (15), (17), (20), and (21), the situation = 0 is considered. Other situations = and r = 0 are not considered, which could be investigated based on the work of Karcol (2021).

Analytical Formulas of Gravitational Effects of a Spherical Zonal Band
Generally, a tesseroid can be treated as a sector of a spherical zonal band (see Fig. 2). In other words, when the longitude ′ of the integration point Q integrates from 1 = 0 to 2 = 2 for a tesseroid, a spherical zonal band is obtained. In this study, the analytical solutions for the GP, GV, GGT, and GC of a spherical zonal band are derived based on these analytical expressions for a tesseroid.
In Appendix E, Laplace's equation is applied to confirm the correctness of the derived analytical expressions for the GGT and GC of a tesseroid and spherical zonal band. The results reveal that these derived analytical expressions for a tesseroid and a spherical zonal band are correct.
In the previous studies, the formulas of the spherical zonal band were derived by the difference between two concentric spherical caps for the GP ( V zb ) ( Visualization of a tesseroid (i.e., shadow part) and a spherical zonal band revised from Fig. 6 of Heck and Seitz (2007). The computation point P( , = 0, r) is located on the polar axis. Q( � , � , r � ) is the integrating point. O is the center of the sphere. Δ , Δ = 2 − 1 , and Δr = r 2 − r 1 are the longitude, colatitude, and radial dimensions of a tesseroid 1 3 and Denker 2019; Deng 2022), radial-radial GGT ( V zz zb ) (Lin et al. 2020;Deng 2022), and radial-radial-radial GC ( V zzz zb ) (Deng 2022). In theory, the numerical values of these previous formulas should be correspondingly equal to those of Eqs. (22), (23), (24), and (26), respectively. The consistency of the analytical expressions for the gravitational effects of a spherical zonal band between Deng (2022) and this paper is presented in Appendix F.
This study goes beyond the previous studies in that the analytical expressions for a spherical zonal band are efficiently obtained from the analytical formulas of a tesseroid by integrating the longitude of the integration point from 0 to 2 when the colatitude of the computation point is = 0 . Meanwhile, the analytical expressions for the V xx zb , V yy zb , V xxz zb , and V yyz zb are derived, in which V xx zb = V yy zb and V xxz zb = V yyz zb .

Analytical Formulas of Gravitational Effects of a Spherical Shell
In general, a spherical shell can be discretized into tesseroids with respect to the longitude, colatitude (or latitude), and radial directions. When the longitude of the integration point ′ integrates from 1 = 0 to 2 = 2 and its colatitude ′ integrates from 1 = 0 to 2 = for a tesseroid, a spherical shell is obtained. In other words, a spherical shell can be treated as a hollow sphere. In this study, the analytical expressions for the GP, GV, GGT, and GC of a spherical shell are derived based on these analytical expressions for a tesseroid when the computation point is located on the polar axis. Letting 1 = 0 , 2 = 2 , 1 = 0 , and 2 = in Eqs. (4), (7), (9), (14), (15), (17), (20) and (21), the analytical expressions for the GP ( sh , and V yyz sh ) of a spherical shell with = 0 are obtained as: Note that these formulas in Eqs. (28)-(33) can be applied to a wider range for the computation point due to the symmetry of the spherical shell. It is obvious that the expressions in Eqs. (28)-(33) are the same as the previous formulas of the spherical shell for the GP, radial GV (Tsoulis 1999;Vaníček et al. 2001Vaníček et al. , 2004Heck and Seitz 2007;Makhloof and Ilk 2008;Grombein et al. 2013;Lin and Denker 2019;Lin et al. 2020), GGT (Makhloof and Ilk 2008;Lin et al. 2020), and GC (Deng and Shen 2018a, b, 2019). Meanwhile, the two pairs of the three components of the GGT and GC satisfy Laplace's equation as V xx sh + V yy sh + V zz sh = 0 and V xxz sh + V yyz sh + V zzz sh = 0 . Thus, the consistencies among these gravitational quantities of a spherical shell can be confirmed.
Analogously, regarding the other GV ( V x sh and V y sh ), GGT ( V xy sh , V xz sh , and V yz sh ), and These components are the zero terms. The Mathematica codes GV_VxVy. nb, GGT_VxyVxzVyz.nb, and GC_VxxxVxxyVxyzVyyxVyyyVzzxVzzy.nb at the link https:// github. com/ xiaol edeng/ analy tical-tesse roid-spher ical-zonal-band can demonstrate these characteristics.

Setup of the Numerical Experiments
In the following numerical experiments, the basic numerical strategy is to compare the relative errors between the calculated values and the reference values of the gravitational effects. The difference from previous studies is that this paper investigates the relative errors between the calculated values and reference values of the gravitational effects of a single tesseroid, which helps to understand the superposition error (or elimination) effect of a spherical zonal band and spherical discretized into tesseroids. By using tesseroids to discretize a spherical zonal band and a spherical shell, there are no relative errors for the zero terms of the GV, GGT, and GC components ( V x , V y , V xy , V xz , V yz , V xxx , V xxy , V xyz , V yyx , V yyy , V zzx , and V zzy ) in the numerical experiments. Thus, these zero terms are not considered in the following experiments.
The general formula of the relative errors in log 10 scale of the gravitational effects is presented as: where F means the relative errors in log 10 scale of the gravitational effects ( V , V r , V zz , V xx , V yy , V zzz , V xxz , and V yyz ) of the tesseroid, spherical zonal band, and spherical shell. F ref represents the reference values of the gravitational effects of the tesseroid, spherical zonal band, and spherical shell. F cal denotes the calculated values of the gravitational effects not only for a single tesseroid but also for the sum of every individual tesseroids forming the whole spherical zonal band and spherical shell.
Based on the analytical expressions for the GP (V), GV ( V r ), GGT ( V zz , V xx , and V yy ), and GC ( V zzz , V xxz , and V yyz ) of a tesseroid in Eqs. To make the process of the calculated gravitational effects of the tesseroid complete, the expressions for the GP, GV, GGT, and GC of a tesseroid in Cartesian integral kernels are listed in Table 1, which are cited from Appendix 1 of Deng and Shen (2019). In this paper, the calculated values of the gravitational effects of the tesseroid in Table 1 are obtained by using the 3D GLQ approach. The GLQ is a commonly used numerical method to calculate the integrals in the expressions for the gravitational and magnetic effects of a tesseroid in previous studies (Ku 1977;Blakely 1995;Asgharzadeh et al. 2007Asgharzadeh et al. , 2008Asgharzadeh et al. , 2018Wild-Pfeiffer 2008;Hirt et al. 2011;Li et al. 2011;Hinze et al. 2013;Du et al. 2015;Rexer and Hirt 2015;Roussel et al. 2015;Uieda et al. 2016;Baykiev et al. 2016 where f ′ , ′ , r ′ means the Cartesian integral kernels of the GP, GV, GGT, and GC in Table 1. W i , W j , and W r k are the weights and N , N , and N r are the integer orders for the spherical longitude, latitude, and geocentric distance, respectively. f i , j , r k denotes the integral kernels of the GP, GV, GGT, and GC with the point i , j , r k . The detailed formulas of the W i , W j , W r k , i , j , and r k can be referred to in Eqs. (23) -(28) and Table 5 of Deng and Shen (2018b).
The numerical values of the parameters for the computation point, tesseroid, spherical zonal band, and spherical shell are listed in Table 2. It should be noted that the values of the r 2 and r 1 for the tesseroid in Table 2 also apply for the spherical zonal band and spherical shell. Regarding the calculated values of the gravitational effects of the tesseroid, the latitudes (i.e., , ′ , 1 , and 2 ) are applied, whereas the colatitudes (i.e., , ′ , 1 , and 2 ) are adopted to compute the reference values of the gravitational effects of the tesseroid, spherical zonal band, and spherical shell. The height of the computation point is the GOCE-type satellite height 260 km (Grombein et al. 2013;Uieda et al. 2016;Deng and Shen 2018a, b), which can largely reduce the near-zone problem (Deng and Shen 2018a, b;Deng 2022), i.e., large relative errors occur when the computation point approaches the near surface of the mass body. Table 1 Detailed expressions for the GP (V), GV ( V r ), GGT ( V zz , V xx , and V yy ), and GC ( V zzz , V xxz , and V yyz ) of a tesseroid in Cartesian integral kernels with the computation point located on the polar axis, where l =

Influence of GLQ Orders on Gravitational Effects of the Tesseroid, Spherical Zonal Band, and Spherical Shell
The chosen orders of the GLQ are important to the accuracy of the calculated gravitational effects of a tesseroid. In previous studies, the computation errors generally reduced with the increased orders of the GLQ, and the determination of the truncated GLQ orders was based on the computation accuracy of the numerical strategy with a spherical shell discretized into tesseroids, e.g., Uieda et al. (2016) Deng et al. (2020). Different from the previous studies, a spherical zonal band discretized into tesseroids is adopted in this section to investigate the influence of the 3D GLQ orders on the GP, GV, GGT, and GC of the tesseroids. Moreover, the relative errors of the gravitational effects of a single tesseroid between the calculated values and reference values are also studied. Note that the term 'nodes' of the GLQ in Deng and Shen (2018b) and Deng et al. (2020Deng et al. ( , 2022 should be replaced by the 'degrees' or 'orders'. The numerical values of the basic parameters are referred to in Table 2. The computation point P is located on the polar axis with its longitude = 0 • , colatitude = 0 • (or latitude = 90 • ), and geocentric distance r = 6, 638, 137 m, which is 260 km above the surface of the spherical shell. The colatitudes of the spherical zonal band are 1 = 10 • and 2 = 11 • . In other words, the grid size of the discretized tesseroids is 1 • × 1 • not only for the spherical zonal band but also for the spherical shell. To calculate the reference values of the gravitational effects of a single tesseroid, the first tesseroid within the spherical zonal band is chosen, i.e., the colatitude boundaries of this single tesseroid are the same as those of the spherical zonal band as 1 = 10 • and 2 = 11 • and its longitude boundaries are 1 = 0 • and 2 = 1 • . Regarding the calculated values of the gravitational effects of this single tesseroid, its latitude boundaries are 1 = 90 • − 2 = 79 • and 2 = 90 • − 1 = 80 • and longitude boundaries are 1 = 0 • and 2 = 1 • . By using Eq. (34) for the single tesseroid, spherical zonal band, and spherical shell, their relative errors of the gravitational effects in log 10 scale with different 3D GLQ orders from (1, 1, 1) to (7, 7, 7) are illustrated in Fig. 3. In addition, the values of these relative errors are listed in Table 3.
In Fig. 3, the points for all GC components (i.e., Band_ V zzz , Band_ V xxz , and Band_ V yyz ; Shell_ V zzz , Shell_ V xxz , and Shell_ V yyz ) almost overlap together at different 3D GLQ orders for the spherical zonal band and spherical shell, i.e., their relative errors in log 10 scale are equal to each other in Table 3. The same rule can be found for all GGT components (i.e., Band_ V zz , Band_ V xx , and Band_ V yy ; Shell_ V zz , Shell_ V xx , and Shell_ V yy ). In general, regarding the gravitational effects of the single tesseroid and spherical zonal band, as the 3D GLQ orders increase, their relative errors decrease and remain at a certain level of about [ −11 , −9 ]. The relative errors of the gravitational effects of the spherical shell decline with the increased 3D GLQ orders. There is one turning point for the situation of the V with the 3D GLQ orders (3, 3, 3) not only for the single tesseroid in Fig. 3a but also for the spherical zonal band in Fig. 3b. The turning point means that the relative errors of the V at the 3D GLQ orders (3,3,3) are smaller than those at the 3D GLQ orders from (4,4,4) to (7,7,7).

Fig. 3
Illustration of the relative errors of the gravitational effects in log 10 scale for (a) a single tesseroid, (b) a spherical zonal band, and (c) a spherical shell with the influence of the 3D GLQ orders from (1, 1, 1) to (7, 7, 7). The grid size of the discretized tesseroids is 1 • × 1 • Regarding the relative errors of the single tesseroid and spherical zonal band in Table 3, the values of the V , V r , V zz , and V zzz are correspondingly equal to each other at different 3D GLQ orders. These numerical results reveal that the superposition error (or elimination) effect does not exist for these gravitational components at different 3D GLQ orders. Meanwhile, the relative errors in log 10 scale of the gravitational components V xx , V yy , V xxz , and V yyz for a spherical zonal band are smaller than those for the single tesseroid at different 3D GLQ orders. This means that the SEEE exists for these gravitational components (i.e., V xx , V yy , V xxz , and V yyz ) when using the strategy of a spherical zonal band discretized into tesseroids. However, at the 3D GLQ orders (3, 3, 3), there exists the abnormal situation that the relative errors of the V xx and V yy for a spherical zonal band are larger than those for the single tesseroid. This may be due to the chosen place of the single tesseroid, which will be investigated using the 3D GLQ orders (3, 3, 3) in the following experiments.
Analogously, for the comparison of the single tesseroid and spherical shell, the relative errors of the spherical shell are larger than those of the single tesseroid with the 3D GLQ orders 1 ≤ n ≤ 4 for the V and V r , 1 ≤ n ≤ 5 for the V xx and V yy , 1 ≤ n ≤ 6 for the V zz , V xxz , and V yyz , and 1 ≤ n ≤ 7 for the V zzz . When the 3D GLQ orders are 5 ≤ n ≤ 7 for the V and V r , 6 ≤ n ≤ 7 for the V xx and V yy , and n = 7 for the V zz , V xxz , and V yyz , the Table 3 Values of the relative errors in log 10 scale using different 3D GLQ orders (n, n, n) (1 ≤ n ≤ 7 ) in Fig. 3 Fig. 3a −2.9 −6.0 −9.2 −9.7 −9.7 −9.7 −9.7 Tess_ V zz in Fig. 3a −2.9 −6.5 −8.  Fig. 3a −2.7 −5.8 −9.2 −9.9 −9.9 −9.9 −9.9 Tess_ V zzz in Fig. 3a Fig. 3b −2.9 −6.0 −9.2 −9.7 −9.7 −9.7 −9.7 Band_ V zz in Fig. 3b −2.9 −6.5 −8.  1 3 relative errors of the spherical shell are smaller than those of the single tesseroid. These results show that the 3D GLQ orders affect the superposition error (or elimination) effect between the single tesseroid and spherical shell. Regarding the spherical shell in relation to the single tesseroid, the superposition error effect exists for lower 3D GLQ orders, while the SEEE is observable for higher 3D GLQ orders.

Influence of Grid Sizes on Gravitational Effects of the Tesseroid, Spherical Zonal Band, and Spherical Shell
In practical applications, the chosen grid size of the tesseroid has an impact on the accuracy of its calculated gravitational effects. In other words, from the numerical aspect, the grid size of the tesseroids to discretize the spherical zonal band and spherical shell affects the relative errors of the gravitational effects of the single tesseroid, spherical zonal band, and spherical shell. Based on Sect. 3.2, different grid sizes are adopted to investigate the influence of the grid sizes on the gravitational effects of the single tesseroid, spherical zonal band, and spherical shell in this section. Due to the fact that the numerical experiment in this section is an extension of that in Sect. 3.2, the numerical settings are the same as those in Sect. 3.2. In order to control the effect of the 3D GLQ orders, the 3D GLQ orders are set as (3, 3, 3). The grid size of the discretized tesseroids are extended from 1 • × 1 • to 5 • × 5 • , 2 • × 2 • , 30 � × 30 � , and 15 � × 15 � . The longitude boundaries of the corresponding single tesseroid are different as [ 0 The relative errors of the gravitational effects in log 10 scale with different grid sizes for the single tesseroid, spherical zonal band, and spherical shell are shown in Fig. 4, in which their values are listed in Table 4.
In Fig. 4 and Table 4, the relative errors overlap together for the two pairs of V xx , V yy ; V xxz , V yyz for the spherical zonal band, and for the two pairs of V zz , V xx , V yy ; V zzz , V xxz , V yyz for the spherical shell. Specifically, in Fig. 4a and Table 4, with the increased grid size of the single tesseroid, the relative errors decrease first and then increase for the GP ( V ), GV ( V r ), and GGT ( V zz , V xx , and V yy ). Regarding the GC ( V zzz , V xxz , and V yyz ) of the single tesseroid, their relative errors decrease as the grid sizes increase. In Fig. 4b for the spherical zonal band, the relative errors of the V , V r , V zz , V xx , and V yy have the same rules as those for the single tesseroid in Fig. 4a. Regarding the V zzz , V xxz and V yyz for the spherical zonal band, their relative errors decrease with the increased grid sizes. In Fig. 4c for the spherical shell, the relative errors of all gravitational components decrease with the increased grid sizes. Regarding some components (i.e., V , V r , V zz , V xx , and V yy ) of the single tesseroid and the spherical zonal band, the smaller the grid size initially, the smaller the calculation errors, but, once a certain threshold is reached, the smaller the grid size, the larger the calculation errors. One possible explanation is that after a certain threshold is reached, other factors generating the errors are greater than the grid size effects, resulting in the appearance of other faults, e.g., the effects of the GLQ orders.
For the comparison of the single tesseroid and spherical zonal band in Table 4, the relative errors of the V , V r , V zz , and V zzz are correspondingly equal to each other at different grid sizes. These numerical results show that the superposition error (or elimination) effect does not exist for these components at different grid sizes. Regarding the V xx and V yy , their relative errors in log 10 scale for the spherical zonal band are larger than those for the single tesseroid by the grid sizes of 1 3 5 • × 5 • , 2 • × 2 • , 1 • × 1 • and 30 � × 30 � . This indicates the existence of the superposition error effect. Regarding the grid size of 15 � × 15 � for the V xx and V yy , their relative errors in log 10 scale for the spherical zonal band (i.e., −10.1 and −10.10 ) are smaller than those for the single tesseroid (i.e., −9.8 and −10.06 ), which indicates the existence of the SEEE. Similarly, by comparing the relative errors of the V xxz and V yyz between the single tesseroid and spherical zonal band, the SEEE happens for the V xxz and V yyz with the grid sizes 5 • × 5 • , 2 • × 2 • , 1 • × 1 • , and 30 � × 30 � . Regarding the grid size 15 � × 15 � , the superposition error effect exists for the V xxz , while the SEEE occurs for the V yyz .  Fig. 4 Visualization of the relative errors of the gravitational effects in log 10 scale for (a) the single tesseroid, (b) a spherical zonal band, and (c) a spherical shell with the influence of the grid sizes as 5 • × 5 • , 2 • × 2 • , 1 • × 1 • , 30 � × 30 � , and 15 � × 15 � . The applied 3D GLQ orders are (3, 3, 3). Other parameters are the same as in Fig. 3 Regarding the single tesseroid and spherical shell in Fig. 4a, c, and Table 4, the relative errors of the GGT ( V zz , V xx , and V yy ) and GC ( V zzz , V xxz , and V yyz ) for the spherical shell are larger than those for the single tesseroid at different grid sizes (i.e., 5 • × 5 • , 2 • × 2 • , 1 • × 1 • , 30 � × 30 � , and 15 � × 15 � ). Meanwhile, this rule can also be found for the relative errors of the GP ( V ) and GV ( V r ) for the spherical shell with the grid sizes 5 • × 5 • , 2 • × 2 • , 1 • × 1 • , and 30 � × 30 � . These numerical results indicate the existence of the superposition error effect for the numerical strategy of a spherical shell discretized into tesseroids. However, for the grid size 15 � × 15 � of the V and V r , the numerical results of the relative errors between the single tesseroid and spherical shell reveal that the SEEE happens in current situation.

Comparison of Gravitational Effects of Every Individual Tesseroid Forming the Whole Spherical Zonal Band and Spherical Shell
In previous numerical experiments in Sects. 3.2 and 3.3, the location of the single tesseroid was based on a specific tesseroid, i.e., the first tesseroid in the longitude direction at different grid sizes. To comprehensively investigate the superposition error (or elimination) effect of the spherical zonal band and spherical shell discretized into tesseroids, the specific Table 4 Values of the points using different gird sizes in Fig. 4 with the 3D GLQ orders (3, 3, 3) single tesseroid is extended to every individual tesseroid forming the whole spherical zonal band and spherical shell. In other words, the relative errors of the gravitational effects (i.e., GP, GV, GGT, and GC) of every individual tesseroid forming the whole spherical zonal band and spherical shell are studied in this section. The numerical settings in this section are the same as those in Sect. 3.3. The 3D GLQ orders are set as (3, 3, 3), and the grid size of the discretized tesseroids is 1 • × 1 • . This simplifies the numerical experiment in this section, in which the effects of the 3D GLQ orders and grid size are explored numerically in Sects. 3.2 and 3.3, respectively. The relative errors of the gravitational effects in log 10 scale for each individual tesseroid forming the whole spherical zonal band are illustrated in Fig. 5a. Note that these relative errors are calculated based on the reference values and calculated values of each individual tesseroid. After subtracting the relative errors of the gravitational effects in log 10 scale for the spherical zonal band discretized using tesseroid (i.e., these values are listed in Table 4 with the grid size 1 • × 1 • ) based on the corresponding numerical values in Fig. 5a, their differences of the relative errors in log 10 scale of the gravitational effects are shown in Fig. 5b. Meanwhile, the statistical information of the values in Fig. 5 is listed in Table 5.
When calculating the relative errors of the gravitational effects for every individual tesseroid forming the whole spherical shell, the ln(0) occurs for the reference values of the gravitational effects (e.g., V, V r , V zz , V xx , and V yy ) with the colatitude of the integration point 1 = 0 • . Thus, the relative errors of the gravitational effects with the first number of the colatitude for the tesseroids forming the whole spherical shell are dropped. The number of the colatitude for every individual tesseroid is [2,180] with the grid size 1 • × 1 • . Then, the relative errors of the gravitational effects in log 10 scale for each individual tesseroid totally forming the spherical shell along the colatitude direction with the first longitude profile are shown in Fig. 6a, where the first longitude profile is chosen to reveal the variation in the relative errors. In this numerical experiment, the first longitude is from 0 • to 1 • due to the grid size of 1 • × 1 • . Similarly, after subtracting the relative errors of the gravitational effects in log 10 scale for a spherical shell discretized into tesseroids (i.e., these values are listed in Table 4 with the grid size 1 • × 1 • ) Fig. 5 Illustration of (a) the relative errors of the gravitational effects in log 10 scale for each individual tesseroid totally forming the whole spherical zonal band; (b) the difference of the relative errors of the gravitational effects in log 10 scale between each individual tesseroid and a spherical zonal band discretized into tesseroids. The x-axis means the counter of the individual tesseroid along the longitude direction. The grid size of the discretized tesseroids is 1 • × 1 • and the 3D GLQ orders are (3, 3, 3) based on the numerical results in Fig. 6a, their differences of the gravitational effects are shown in Fig. 6b. The statistical information of the values in Fig. 6 is listed in Table 6.
Moreover, the relative errors in log 10 scale of the gravitational effects for every individual tesseroid forming the whole spherical shell are also calculated and their statistical values are listed in Table 7. After subtracting the relative errors of the gravitational effects in log 10 scale for a spherical shell, the differences of the relative errors in log 10 scale of the gravitational effects are obtained, and their statistical values are also listed in Table 7.
In Fig. 5a, the relative errors of the V , V r , V zz , and V zzz appear as distinct straight lines, which indicates that their values are not affected by changes in longitude of each individual tesseroid. Regarding the V xx , V yy , V xxz , and V yyz , their relative errors exhibit periodic changes as longitude changes. The period of change is about 180 • . In Fig. 5b, the differences of the V , V r , V zz , and V zzz are about zero, which reveals that the superposition error (or elimination) effect does not exist for these components by using the strategy of a spherical zonal band discretized into tesseroids. Regarding the V xx , V yy , V xxz , and V yyz , their differences are related to the position of the individual tesseroid along the longitude direction. From Table 5, the mean values of the differences are about 0.3 for the V xx and V yy and −0.1 for the V xxz and V yyz . These numerical values reveal that on the overall average, the SEEE occurs for the V xx and V yy , and the superposition error effect exists for the V xxz and V yyz .
In Fig. 6a, as the colatitude increases, all curves first drop sharply, then oscillates erratically, and finally rise. In Fig. 6b, except for a small number of the values below zero at the beginning, all subsequent values are greater than zero. These results reveal that the SEEE occurs for the spherical shell discretized into tesseroids in the most cases. This effect can be confirmed not only by the mean values of the differences in Table 6 as 1.8 for the V , 3.7 for the V r , 5.3 for the V zz , 4.8 for the V xx , 5.2 for the V yy , 6.9 for the V zzz , V xxz , and 7.0 for the V yyz , but also by the mean values of the differences in Table 7 as 1.8 for the V , 3.7 for the V r , 5.3 for the V zz , 5.1 for the V xx , V yy , and 6.9 for the V zzz , V xxz , V yyz .

Conclusions and Outlook
Previous studies investigated the gravitational effects of a tesseroid by different numerical approaches. Meanwhile, the integrals in the expressions for the gravitational effects of a tesseroid were often considered not to be integrated analytically. In this contribution, the analytical expressions for the gravitational effects (i.e., GP, GV, GGT, and GC) of a tesseroid are derived when the computation point is located on the polar axis. Based on the relation between the tesseroid and the spherical zonal band, the simpler Fig. 6 Visualization of (a) the relative errors of the gravitational effects in log 10 scale for each individual tesseroid totally forming the spherical shell along the colatitude direction with the first longitude profile; (b) the difference of the relative errors of the gravitational effects in log 10 scale between each individual tesseroid and a spherical shell discretized into tesseroids. The x-axis means the colatitude of the discretized tesseroid. Other parameters are the same as in Fig. 5 Table 6 Statistical information of the values in Fig. 6, where other information is the same in Table 5 Quantity in figure Fig. 6a − 11.7 − 5.5 − 9.9 0.8 9.9 V zz in Fig. 6a − 11.8 − 3.6 − 9.5 1.0 9.6 V xx in Fig. 6a − 12.0 − 3.9 − 9.1 1.2 9.1 V yy in Fig. 6a − 13.3 − 5.5 − 9.4 1.0 9.5 V zzz in Fig. 6a − 12.6 − 4.1 − 9.4 1.1 9.4 V xxz in Fig. 6a − 11.8 − 4.1 − 9.3 1.0 9.4 V yyz in Fig. 6a − 11.7 − 5.0 − 9.4 0.9 9.5 V in Fig. 6b − 1.8 3.7 1.8 0.8 2.0 V r in Fig. 6b − 0.7 5.5 3.7 0.8 3.8 V zz in Fig. 6b − 0.7 7.5 5.3 1.0 5.3 V xx in Fig. 6b − 0.3 7.7 4.8 1.2 4.9 V yy in Fig. 6b 1.2 9.0 5.2 1.0 5.3 V zzz in Fig. 6b 1.7 10.1 6.9 1.1 7.0 V xxz in Fig. 6b 1.6 9.3 6.9 1.0 7.0 V yyz in Fig. 6b 2.5 9.3 7.0 0.9 7.1 analytical expressions for the gravitational effects of a spherical zonal band are derived in comparison with previous studies. In addition, the analytical expressions for the gravitational effects of a spherical shell are presented based on the relation between the tesseroid and the spherical shell.
These new analytical formulas of the GGT and GC of the tesseroid and spherical zonal band are confirmed to be correct by using Laplace's equation. Meanwhile, the consistency of the analytical expressions for the gravitational effects of a spherical zonal band between Deng (2022) and this paper is confirmed in Appendix F.
The influence of the 3D GLQ orders from (1, 1, 1) to (7, 7, 7) on the gravitational effects of the single tesseroid, spherical zonal band, and spherical shell is investigated with the grid size 1 • × 1 • . As the 3D GLQ orders increase, the relative errors of the gravitational effects of the single tesseroid and spherical zonal band decrease and remain at a certain level. For the spherical shell, the relative errors decline with the increased 3D GLQ orders. Comparing the single tesseroid and spherical zonal band, the superposition error (or elimination) effect does not exist for these gravitational components (i.e., V , V r , V zz , and V zzz ) at different 3D GLQ orders. Regarding other gravitational components (i.e., V xx , V yy , V xxz , and V yyz ), the SEEE mainly occurs. For the single tesseroid and spherical shell, the 3D GLQ orders have an impact on the superposition error (or elimination) effect. Specifically, the superposition error effect exists for lower 3D GLQ orders, and the SEEE exists for higher 3D GLQ orders.
The influence of the grid sizes on the gravitational effects of the single tesseroid, spherical zonal band, and spherical shell is investigated with the 3D GLQ order (3, 3, 3). As the grid sizes decrease from 5 • × 5 • to 15 � × 15 � , the relative errors decline first and then increase for the V , V r , V zz , V xx , and V yy of the single tesseroid and spherical zonal band. The relative errors decline with the finer grid sizes not only for the V zzz , V xxz , and V yyz of the single tesseroid and spherical zonal band but also for all gravitational effects (i.e., V , V r , V zz , V xx , V yy , V zzz , V xxz , and V yyz ) of the spherical shell. For the comparison between the single tesseroid and spherical zonal band, the superposition error (or Table 7 Statistical information of the values for the relative errors in log 10 scale of the gravitational effects for every individual tesseroid forming the whole spherical shell and the differences of the relative errors in log 10 scale of the gravitational effects between every individual tesseroid and the spherical shell

Parameter
Quantity Min Max Mean STD RMS Relative errors V − 12.0 − 6.6 − 10.2 0.7 10.2 V r − 11.7 − 5.5 − 9.9 0.8 9.9 V zz − 11.8 − 3.6 − 9.5 1.0 9.6 V xx − 14.3 − 3.4 − 9.3 1. 1.7 10.1 6.9 1.1 7.0 V xxz − 0.3 11.9 6.9 1.0 7.0 V yyz − 0.3 12.5 6.9 1.0 7.0 elimination) effect does not occur for these gravitational components (i.e., V , V r , V zz , and V zzz ) with different grid sizes. Regarding the V xx and V yy , the superposition error effect exists with the grid sizes 5 • × 5 • , 2 • × 2 • , 1 • × 1 • , and 30 � × 30 � and the SEEE occurs with the grid size 15 � × 15 � . On the contrary, regarding the V xxz and V yyz , the SEEE happens with the grid sizes 5 • × 5 • , 2 • × 2 • , 1 • × 1 • , and 30 � × 30 � . The superposition error effect occurs for the V xxz and V yyz with the grid size 15 � × 15 � . Regarding the comparison between the single tesseroid and spherical shell, numerical results show that the superposition error effect happens for the gravitational effects of the spherical shell discretized into tesseroids with different grid sizes, except for the V and V r with the grid size 15 � × 15 � . The superposition error (or elimination) effect is investigated for the gravitational effects of every individual tesseroid forming the whole spherical zonal band and spherical shell with the grid size 1 • × 1 • and 3D GLQ orders (3, 3, 3). Numerical results show that the superposition error (or elimination) effect does not exist for the V , V r , V zz , and V zzz by using the strategy of a spherical zonal band discretized into tesseroids. Regarding the V xx , V yy , V xxz , and V yyz of the spherical zonal band, the SEEE occurs for the V xx and V yy , and the superposition error effect exists for the V xxz and V yyz on the overall average. Regarding a spherical shell discretized into tesseroids, the SEEE occurs for this strategy in most cases.
In short, these numerical experiments confirm the existence of the SEEE hypothesis. When using the tesseroids in practical applications of gravity field modeling, the local error assessment by using the analytical formulas of the GP, GV, GGT, and GC of a tesseroid of this paper will be performed necessarily in the future study. Due to the existence of the SEEE and superposition error effect, the benchmark of a spherical zonal band and a spherical shell discretized into small mass bodies should be employed with caution in gravity field modeling in geodesy and geophysics. Their numerical results only pertain to the global situation as a whole spherical zonal band or a whole spherical shell, not the local one.
The new strategy of a reference tesseroid with respect to its calculated tesseroid can be proposed in comparison with the commonly used strategies of a spherical zonal band and a spherical shell discretized into tesseroids. The mathematical derivations using the spherical coordinates , , and r in the paper can also be performed by using the spherical polar coordinates and (Novák et al. 2017(Novák et al. , 2019. Meanwhile, the coordinate system transformation of the higher-order gravitational gradients will be investigated in the near future. In future studies, these new derived analytical formulas of the gravitational effects of a tesseroid can serve as reference values for the evaluation of the calculated gravitational effects of the tesseroid by using not only the above-mentioned numerical methods but also other numerical approaches, e.g., spherical harmonics (Ramillien 2017;Baykiev et al. 2020;Šprlák et al. 2020). Future research will investigate the comparison of these numerical approaches to calculate the volume or surface integrals of the tesseroid. Compared with the spherical zonal band and spherical shell, the variation in the relative errors of the gravitational effects of the tesseroid can be revealed from a more subtle perspective. The SEEE can be extended to an arbitrary mass body that can be discretized into smaller mass bodies not only in the spatial domain but also in the spectral domain. On the basis of this paper, the higher-order gravitational potential gradient (e.g., the fourth-order gravitational potential gradient (Deng and Ran 2021) and firstorder derivatives of the invariants of the GGT ) of a tesseroid can also be derived analytically when the computation point is located on the polar axis. Finally, based on the results of this study in the gravity field, it can be extended to the magnetic field. In other words, the analytical expressions for the magnetic effects of a tesseroid can be derived in the magnetic field.

Appendix A Derivation of the Analytical Expression for the GP of a Tesseroid
When the computation point is located on the polar axis and = 0 , Eq. (3)  where l 2 = √ r 2 + r � 2 − 2rr � cos 2 and l 1 = √ r 2 + r � 2 − 2rr � cos 1 .
In Eq. (40), ∫ r 2 r 1 r ′ l 2 dr ′ is obtained as: and ∫ r 2 r 1 −r � l 1 dr � is obtained as: After substituting Eqs. (41) and (42) into Eq. (40), the analytical formula of the GP of a tesseroid when the computation point is located on the polar axis is obtained in Eq. (4).
Regarding the analytical expressions for the V x and V y of a tesseroid, their derivations are presented in the Mathematica code GV_VxVy.nb at the link https:// github. com/ xiaol edeng/ analy tical-tesse roid-spher ical-zonal-band.