Evaluation of gravitational curvatures for a tesseroid and spherical shell with arbitrary-order polynomial density

In recent years, there are research trends from constant to variable density and low-order to high-order gravitational potential gradients in gravity field modeling. Under the research circumstances, this paper focuses on the variable density model for gravitational curvatures (or gravity curvatures, third-order derivatives of gravitational potential) of a tesseroid and spherical shell in the spatial domain of gravity field modeling. In this contribution, the general formula of the gravitational curvatures of a tesseroid with arbitrary-order polynomial density is derived. The general expressions for gravitational effects up to the gravitational curvatures of a spherical shell with arbitrary-order polynomial density are derived when the computation point is located above, inside, and below the spherical shell. When the computation point is located above the spherical shell, the general expressions for the mass of a spherical shell and the relation between the radial gravitational effects up to arbitrary-order and the mass of a spherical shell with arbitrary-order polynomial density are derived. The influence of the computation point’s height and latitude on gravitational curvatures with the polynomial density up to fourth-order is numerically investigated using tesseroids to discretize a spherical shell. Numerical results reveal that the near-zone problem exists for the fourth-order polynomial density of the gravitational curvatures, i.e., relative errors in log10\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\log _{10}$$\end{document} scale of gravitational curvatures are large than 0 below the height of about 50 km by a grid size of 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'\times 15'$$\end{document}. The polar-singularity problem does not occur for the gravitational curvatures with polynomial density up to fourth-order because of the Cartesian integral kernels of the tesseroid. The density variation can be revealed in the absolute errors as the superposition effects of Laplace parameters of gravitational curvatures other than the relative errors. The derived expressions are examples of the high-order gravitational potential gradients of the mass body with variable density in the spatial domain, which will provide the theoretical basis for future applications of gravity field modeling in geodesy and geophysics.


Introduction
Generally, different structures of the Earth (e.g., atmosphere, ice, water, topography, sediment, crust, mantle, and core) have different densities. For example, the density from the surface to the inner core of the Earth varied from 1020.0 to 13088.5 kg m −3 in the preliminary reference Earth model (Dziewonski and Anderson 1981) and from 2651.0 to 13012.2 kg m −3 in the radial ek137 model (Kennett 2020). The density hypothesis and digital models (e.g., CRUST1.0 (Laske et al. 2013) and UNB_TopoDens (Sheng et al. 2019)) X.-L. Deng arbitrary-order polynomial density (Zhang and Jiang 2017), GV and GGT in the Fourier domain with the depth-dependent polynomial density (Wu and Chen 2016;Wu 2018), GGT with the depth-dependent density (Jiang et al. 2018), and GP, GV, and GGT with the depth-dependent nth-order polynomial density (Karcol 2018;Fukushima 2018b). Similarly, the gravitational effects of a spherical shell were derived in the form of the polynomial density, i.e., the GP and GV with the radial fifth-order polynomial density (Karcol 2011) and the GP, GV, and GGT with the linear, quadratic, and cubic order polynomial density (Lin et al. 2020). Regarding a tesseroid (Anderson 1976;Heck and Seitz 2007), its gravitational effects up to the GGT were derived in the form of the linear density in Fukushima (2018a), Lin and Denker (2019), and Soler et al. (2019). Lin et al. (2020) presented the general expressions for the GP, GV, and GGT of a tesseroid in the form of the polynomial density.
With the first measurement of radial gravitational curvatures (or gravity curvatures, GC, third-order derivatives of the gravitational potential) using the atom interferometry sensors in the laboratory condition (Rosi et al. 2015), the subsequent study of the GC was focused on their theory and simulation applications in geodesy and geophysics (Šprlák and Novák 2015Hamáčková et al. 2016;Sharifi et al. 2017;Novák et al. 2017;Pitoňák et al. 2017Pitoňák et al. , 2018Pitoňák et al. , 2020Shen 2018a, b, 2019;Novák et al. 2019;Du and Qiu 2019;Romeshkani et al. 2020Romeshkani et al. , 2021. The GC were more sensitive to field sources than the loworder gravitational effects (i.e., GP, GV, and GGT) (Heck 1984;Du and Qiu 2019). Numerical results showed that they were more sensitive to local characteristics of the gravity field (Romeshkani et al. 2020), can better describe areas with large terrain undulations and the junction of the land and sea (Deng and Shen 2019), and can more accurately identify the near-underground shallow density structure in the spatial domain (e.g., caves, salt domes, sediment base morphology, continental margins, and hidden fault systems) (Novák et al. 2019) than the low-order gravitational effects. These advantages of the GC are important for studying the distribution of materials and geological structures within the Earth.
Currently, the expressions for the GC of a tesseroid and a spherical shell with the constant density were derived in Deng and Shen (2018a). Meanwhile, Deng (2022) recently derived the analytical solutions for the gravitational curvatures of a spherical cap and spherical zonal band with the constant density. With the increased depth inside the Earth, the actual density shows a complex rule of variable density. The constant density is not enough to represent the actual density variation of the earth for the GC. In other words, the existing constant density assumption cannot meet the demand for a real complex environment with the variable density.
Thus, it is essential to investigate their variable density forms, especially the polynomial density. This study focuses on the derivation of arbitrary-order (i.e., N ≥ 0) polynomial density for the GC of a tesseroid and spherical shell.
This study goes beyond the previous studies in that the general formulae of the GC of a tesseroid with N th-order polynomial density and the general analytical expressions for the GC of a spherical shell with N th-order polynomial density are derived. In addition, we derive the general analytical expressions for gravitational effects (e.g., GP, GV, and GGT) of a spherical shell with N th-order polynomial density when the computation point is located above, inside, and below the spherical shell. The derived polynomial density expressions of the GC would more accurately describe the complex density characteristics of the Earth's interior than the constant density expressions, which will lay a theoretical foundation for the applications of gravity field modeling in geodesy and geophysics.
This paper is organized as follows. Section 2 presents the theoretical aspects, in which Sect. 2.1 provides the general formula of the GC of a tesseroid with N th-order polynomial density. In Sect. 2.2, the analytical expressions for gravitational effects up to the GC of a spherical shell with N th-order polynomial density are derived when the computation point is located above, inside, and below the spherical shell. Section 3 performs the numerical experiments, where the setup of the numerical experiments is provided in Sect. 3.1 and the effects of density values on the gravitational effects of a spherical shell are studied in Sect. 3.2. The influences of the computation point's height and latitude on the GC with the polynomial density up to fourth-order are investigated in Sects. 3.3 and 3.4, respectively. Finally, the conclusions of this paper and outlooks on future research work are summarized in Sect. 4.

Theoretical aspects
In this section, the theoretical contents are presented in detail. Based on the previous works of Shen (2018b, 2019) and Lin et al. (2020), the formula of the GC of a tesseroid with arbitrary-order polynomial density is derived in Sect. 2.1. Section 2.2 derives the expressions for the GC of a spherical shell with arbitrary-order polynomial density based on the work of Lin et al. (2020). The general formula of the mass for a spherical shell with arbitrary-order polynomial density is derived in Sect. 2.2. Moreover, Sect. 2.2 derives the general relation between the radial gravitational effects up to arbitrary-order and the mass of a spherical shell with arbitrary-order polynomial density when the computation point P is outside the spherical shell.  Heck and Seitz (2007). P and Q are the computation (or field, observation) and running integration points with their longitudes (λ and λ ), latitudes (ϕ and ϕ ), and geocentric distances (r and r ), respectively. Δλ = λ 2 − λ 1 , Δϕ = ϕ 2 − ϕ 1 , and Δr = r 2 − r 1 are the longitude, latitude, and radial dimensions of a tesseroid, respectively. P(x, y, z) is in the local north-oriented frame system with x, y, and z pointing to the north, east, and radial up directions, respectively

The GC of a tesseroid with arbitrary-order polynomial density
Regarding the variable density model for a tesseroid (see Fig. 1), Lin et al. (2020) derived the GP, GV, and GGT of a tesseroid with a polynomial density model up to N th-order in the vertical direction and evaluated these gravitational effects with a polynomial density model up to cubic order (i.e., N = 3) in the numerical experiments. Due to the potential benefits of the GC in geodesy and geophysics revealed in the introduction, it is essential to investigate the variable polynomial density model for the GC of a tesseroid. In this section, based on the previous works of Shen (2018b, 2019) and Lin et al. (2020), the general expression for the GC of a tesseroid with N th-order polynomial density is derived, and the detailed components of the GC are provided as well.
Replacing the constant density ρ in Eq.
, Δ x , Δ y , and Δ z in Eqs. (3-6) are from Eqs. (4) and (15) of Grombein et al. (2013). (λ, ϕ, r ) and (λ , ϕ , r ) are the spherical longitudes, latitudes, and geocentric distances of the computation point P and integration point Q, respectively. The computation point P is in the local north-oriented frame system with its coordinate origin at the computation point and x, y, and z pointing to the north, east, and radial up directions, respectively.
The formula of the polynomial density ρ(r ) of a tesseroid is presented as: where n and N are the integers and n = 0, 1, 2, 3, ..., N with N ≥ 0. ρ n is the polynomial density coefficient with the unit kg m −(n+3) . For example, when n = 0, ρ 0 is the constant density with the unit kg m −3 . Another method to obtain the expressions for the ρ n can be referred in Table 1 of Lin et al. (2020). Substituting Eq. (8) into (1), the optimized formula of the GC (V αβγ ) of a tesseroid in Cartesian integral kernels with N th-order polynomial density is presented as: Including the low-order gravitational potential gradients (i.e., V , V α , and V αβ with α, β = x, y, z) presented in Eq. (3) of Lin et al. (2020), the general expressions for the gravitational effects up to the GC of a tesseroid in Cartesian integral kernels with N th-order polynomial density are presented as: where I (λ , ϕ , r ), I α (λ , ϕ , r ), and I αβ (λ , ϕ , r ) are the Cartesian integral kernels of the GP, GV, and GGT of a tesseroid, which can be referred to in Eq. (21) of Grombein et al. (2013) and Eqs. (A.1-A.10) of Deng and Shen (2019). The (V tess ) n , (V tess α ) n , (V tess αβ ) n , and (V tess αβγ ) n are the nthorder polynomial density parts of the GP, GV, GGT, and GC of a tesseroid in Cartesian integral kernels, respectively. In order to show the completeness of the expressions for gravitational effects up to the GC, the detailed expressions for a tesseroid in Cartesian integral kernels having a polynomial density model up to N th-order are presented in Table 1.
The expressions in the second column of Table 1 for the GP, GV, and GGT are the same as those in Eq. (21) of Grombein et al. (2013) and the GC are the same as those in Eqs. (A.11-A.20) of Deng and Shen (2019). It should be noted that low-order gravitational potential gradients (i.e., GP (V ), GV (V α ), and GGT (V αβ ) with α, β = x, y, z) presented in Table 1 are consistent with those in Table 2 and Eq.

Gravitational effects up to the GC of a spherical shell with arbitrary-order polynomial density
The detailed derivation process of analytical expressions for the GP of a spherical shell having a polynomial density model from constant up to cubic order was presented in the electronic supplementary material of Lin et al. (2020). Furthermore, the general expressions for gravitational quantities Table 1 Detailed expressions for the GP (V ), GV (V α ), GGT (V αβ ), and GC (V αβγ ) (α, β, γ = x, y, z) of a tesseroid in Cartesian integral kernels having a polynomial density model up to N th-order, where = r 2 + r 2 − 2r r (sin ϕ sin ϕ + cos ϕ cos ϕ cos(λ − λ)), Δ x = r cos ϕ sin ϕ − sin ϕ cos ϕ cos(λ − λ) , Δ y = r cos ϕ sin(λ − λ), Δ z = r sin ϕ sin ϕ + cos ϕ cos ϕ cos(λ − λ) − r , and κ n = r n+2 cos ϕ (i.e., GP, radial GV, radial-radial GGT, and radial-radialradial GC) and mass of a spherical shell with N th-order polynomial density are derived in this section. A spherical zonal band (or a spherical layer in Fig. 15 of MacMillan (1930), Fig. 1.1 of Tsoulis (1999), and a circular ring element in Fig. E1 of Lin et al. (2020)) is shown in Fig. 2. When the spherical distance θ of a spherical zonal band is integrated from θ = 0 to π , a spherical shell can be obtained (Vaníček et al. , 2004. The formula of the GP of a spherical shell at the computation point P is given as (MacMillan 1930;Tsoulis 1999; where ρ(r ) is the polynomial density up to N th-order of a spherical shell, which is assumed as the same as that of a tesseroid in Eq. (8). r 1 and r 2 are the inner and outer radii of the spherical shell. r ∈ [r 1 , r 2 ] is the geocentric distance, α ∈ [0, 2π ] is the azimuth, and θ ∈ [0, π] is the spherical distance of the integration point in the local polar coordinate of the spherical coordinate systems. is the Euclidean distance between the computation point P and integration point. Substituting Eqs. (8) and (14) into Eq. (12), the expression for the GP (V (r )) of a spherical shell is simplified as: where the computation point P is located above (i.e., r > r 2 ), inside (i.e., r 1 ≤ r ≤ r 2 ), and below (i.e., r < r 1 ) the spherical shell. The detailed derivation process of analytical expressions for the GP of a spherical shell is presented in Appendix 1.
Combining Eqs. (47), (48), and (49) in Appendix 1 together, the analytical expressions for the GP (V (r )) of a spherical shell with N th-order polynomial density can be presented as: where the analytical expressions for the GP of a spherical shell with different nth-order polynomial density parts beginning from constant (i.e., (V sh ) n with n = 0, 1, 2, 3, ..., N ) are listed in Table 2 to reveal the detailed variation of the polynomial density parts. By performing the differentiation with respect to the geocentric distance r of the computation point at one, two, and three times on the analytical expressions for the GP (V (r )) of a spherical shell in Eq. (16), the analytical expressions for radial GV (V z (r )), radial-radial GGT (V zz (r )), and radialradial-radial GC (V zzz (r )) of a spherical shell with N th-order polynomial density are derived as: where the analytical expressions for radial GV, radial-radial GGT, and radial-radial-radial GC of a spherical shell with different nth-order polynomial density parts (i.e., (V sh z ) n , (V sh zz ) n , and (V sh zzz ) n with n = 0, 1, 2, 3, ..., N ) are similarly listed in Tables 3, 4, and 5, respectively.
The nonzero components of the GGT (i.e., V x x (r ), V yy (r ), and V zz (r )) and GC (i.e., V x xz (r ), V yyz (r ), and V zzz (r )) satisfy Laplace's equation with the computation point locating above and below the spherical shell (i.e., r > r 2 and r < r 1 ) in Eqs. (20-21) Poisson's equation, and Laplace's equation with the computation point locating inside the spherical shell (i.e., r 1 < r < r 2 ) in Eqs. (22-23) as: where V x x (r ) = V yy (r ) and V x xz (r ) = V yyz (r ). Substituting Eqs. (18)(19) into Eqs. (20-23), the analytical expressions for other nonzero components of the GGT and GC of a spherical shell with N th-order polynomial density are presented as: (25) where the analytical expressions for other nonzero GGT and GC of a spherical shell with different nth-order polynomial density parts (i.e., (V sh Table 3 Analytical expressions for radial GV of a spherical shell with different nth-order polynomial density parts (i.e., (V sh z ) n ) Quantity Above (r > r 2 ) I n s i d e ( r 1 ≤ r ≤ r 2 ) B e l o w ( r < r 1 ) Other parameters are the same as in Table 2 Table 4 Analytical expressions for radial-radial GGT of a spherical shell with different nth-order polynomial density parts (i.e., (V sh zz ) n ), where the computation point P is located above (r > r 2 ), inside (r 1 < r < r 2 ), and below (r < r 1 ) the spherical shell Table 5 Analytical expressions for radial-radial-radial GC of a spherical shell with different nth-order polynomial density parts (i.e., (V sh zzz ) n ) Quantity Above (r > r 2 ) I n s i d e ( r 1 < r < r 2 ) B e l o w ( r < r 1 ) Other parameters are the same as in Table 4 with n = 0, 1, 2, 3, ..., N ) are listed in Tables 6 and 7, respectively.
When the computation is located on the top (i.e., r = r 2 ) or bottom (i.e., r = r 1 ) of the spherical shell, there are discontinuities as the density jumps. Regarding the GP component V in Eq. (16) or Table 2 and GV component V z in Eq. (17) or Table 3, r = r 2 and r = r 1 are applied for the situation where the computation point is inside the spherical shell (i.e., r 1 ≤ r ≤ r 2 ) (Fukushima 2018a, Eqs. (46) and (47)). For the GGT components (V zz in Eq. (18) and Table 4; V x x and V yy in Eq. (24) and Table 6) and GC components (V zzz in Eq. (19) and Table 5; V x xz and V yyz in Eq. (25) and Table 7), the computation point that is located on the top or bottom of the spherical shell is slightly moved inside the spherical shell. This strategy is the same as that in Fukushima (2018a). In detail, the geocentric distance r of the computation point when it is on the top (r top ) or bottom (r bottom ) of the spherical shell can be presented as: where R is the radius of the reference sphere. is the machine epsilon as = 2 −53 ≈ 1.11 × 10 −16 for double precision and = 2 −113 ≈ 9.63 × 10 −35 for quadruple precision (Fukushima 2012(Fukushima , 2018a. Another strategy to treat the special condition of the computation point located on the top of the spherical shell is to X.-L. Deng Table 6 Analytical expressions for other nonzero gravity gradient tensor of a spherical shell with different nth-order polynomial density parts (i.e., Other parameters are the same as in Table 4 Table 7 Analytical expressions for other nonzero gravitational curvatures of a spherical shell with different nth-order polynomial density parts Other parameters are the same as in Table 4 move the computation point 1 m above the spherical shell for the evaluation of the GGT components Lin et al. 2020). This strategy can also be applied to the GC components of this special condition. Furthermore, the general formula of the mass for a spherical shell with N th-order polynomial density is provided as: where the formulae of the mass for a spherical shell with different nth-order polynomial density parts (i.e., (M sh ) n ) are listed in Table 8. The expressions for the (M sh ) 0 , (M sh ) 1 , (M sh ) 2 , and (M sh ) 3 in Table 8 are the same as those in Eqs. (E6b), (E7b), (E8b), and (E9b) of Lin et al. (2020), respectively.
When the computation point P is outside the spherical shell (i.e., r > r 2 ), the relation between Eq. (28) and Eqs. (16)(17)(18)(19) can be found as: Based on Eqs. (30-32), the general relation between the radial gravitational effects (V t (r )) up to arbitrary-order and the mass (M) of a spherical shell with N th-order polynomial density where the computation point P is above the spherical shell (i.e., r > r 2 ) is derived as: where t is the nonzero positive integer (i.e., t = 1, 2, 3, ...), which means the number of derivation in the radial direction of the computation point.
In this section, the work goes beyond the previously cited studies (MacMillan 1930;Tsoulis 1999;Deng and Shen 2018a;Lin et al. 2020) in that we derive the general analytical expressions for the GP (V (r ) in Eq. (16) and Table 2), radial GV (V z (r ) in Eq. (17) and Table 3 (24) and Table 6, and V zz (r ) in Eq. (18) and Table 4), and nonzero GC (V x xz (r ), V yyz (r ) in Eq. (25) and Table 7, and V zzz (r ) in Eq. (19) and Table 5) of a spherical shell with N th-order polynomial density when the computation point P is located above, inside, and below the spherical shell. Moreover, the consistent expressions for the GP, GV, GGT, and GC of a spherical shell are listed in Table 9.
The general formula of the mass for a spherical shell having a polynomial density model up to N th-order is derived in Eq. (28) and Table 8 as well. Moreover, we derive the general relation between the radial gravitational effects up to arbitrary-order and the mass of a spherical shell with N thorder polynomial density in Eq. (33) when the computation point P is outside the spherical shell.

Numerical investigations
In this section, the numerical experiments are performed to investigate the effects of the GC of the tesseroids and spherical shell with the polynomial density up to fourth-order. The setup of the numerical experiments, including numerical strategy, layout, method, and chosen values, is described in Sect. 3.1. The effects of selected numerical density values on the gravitational effects of a spherical shell are investigated in Sect. 3.2. Sections 3.3 and 3.4 examine the effects of the computation point's height and latitude on the GC with the polynomial density up to the fourth-order.

Setup of the numerical experiments
In the following experiments, the commonly used numerical strategy of using tesseroids to discretize a whole spherical shell is adopted due to the simple analytical solutions of the gravitational effects of a spherical shell. As the continuous work of Deng and Shen (2018a), the numerical layouts of the experiments are the same as each other. The main difference lies in that the focus of the experiments in this paper is the effect of variable density models on the gravitational effects with the influence of height and latitude. This effect is revealed by the relative errors between the calculated gravitational effects of the discretized tesseroids and the analytical solutions of the gravitational effects of a spherical shell.
Specifically, the relative errors ((δ F) N ) of the tesseroids discretizing the whole spherical shell with N th-order polynomial density are given in log 10 scale, where F represents the function of the nonzero gravitational effects from the GP to GC, i.e., The reference values are the analytical nonzero gravitational effects of a spherical shell, which can be calculated from Tables 2, 3, 4, 5, 6 and 7. The calculated values are the sum of the nonzero gravitational effects of the discretized tesseroids forming the total spherical shell.
Laplace's equation of the GGT and GC is applied to reveal the internal numerical quality as: where (δΔL 1 ) N and (δΔL 2 ) N are the Laplace parameters of the GGT and GC with N th-order polynomial density. They are theoretically equal to zero and denoted as absolute errors in the numerical experiments. In the following experiments, the numerical calculation with respect to relative and absolute errors is performed in quadruple precision. The 3D Gauss-Legendre quadrature method is applied to obtain the numerical values of gravitational effects up to  Table 2 Equation (15a) in Lin et al. (2020) GP with linear, quadratic, and cubic density (V sh ) 1 , (V sh ) 2 , (V sh ) 3 in Table 2 Equations (16a), (17a), (18a) in Lin et al. (2020) GV with constant density (V sh z ) 0 in Table 3 Equation (15b) in Lin et al. (2020) GV with linear, quadratic, and cubic density ( Table 3 Equations (16b), (17b), (18b) in Lin et al. (2020) GGT with constant density (V sh zz ) 0 in Table 4 Equation (15c) in Lin et al. (2020) GGT with linear, quadratic, and cubic density ( Table 4 Equations (16c), (17c), (18c) in Lin et al. (2020) GC with constant density (V sh xxz ) 0 ,(V sh yyz ) 0 in Table 7 and (V sh zzz ) 0 in Table 5 Equations (8) and (9) in Deng and Shen (2018a) where F means the function of detailed expressions for gravitational effects up to the GC in Table 1. I (λ i , ϕ j , r k ) denotes the function of detailed integral kernels of gravitational effects up to the GC in Table 1.
is the integration domain; N λ , N ϕ , and N r are the integer degrees; W λ i , W ϕ j , and W r k are the weights; x i , x j , and x k are the nodes with respect to the longitude, latitude, and radial directions, respectively. With larger degrees, better computational accuracy can be obtained, whereas the computation time will be longer. Based on Sect. 4.2 of Deng and Shen (2018b), the degrees of N λ , N ϕ , and N r are chosen as 2, 2, and 2 to present the acceptable numerical results at the height of 260 km, and related weights and nodes can be referred in Table 4 of Wild-Pfeiffer (2008) and Table 5 of Deng and Shen (2018b).
The numerical values of the parameters of the spherical shell and tesseroid are listed in Table 10. It should be noted that the constant density of the spherical shell and tesseroid is the same as each other as ρ 0 = 2670 kg m −3 , which is often assumed as the mean density for the topography of the earth (Heiskanen and Moritz 1967;Hinze 2003;Tenzer et al. 2021).

Effects of density values on gravitational effects of a spherical shell
The chosen numerical values of the density coefficient ρ n in the different expressions for gravitational effects (e.g., V , V z , V zz , and V zzz ) of a spherical shell and tesseroid are important for the evaluation of the gravitational effects of a spherical shell and tesseroid. Due to the analytical solutions of the gravitational effects of a spherical shell as the reference values for discretized tesseroids, the effects of density values on the gravitational effects of a spherical shell are investigated in this section.
In the numerical experiments, the values of gravitational effects of a spherical shell with the chosen density coefficient ρ n (n ≥ 1) are assumed to be at the same level as these of a Fig. 3 Visualization of the reference values in log 10 scale for the GP ((V sh ) 0 blue curve), radial GV ((V sh z ) 0 red curve), radial-radial GGT ((V sh zz ) 0 green curve), and radial-radial-radial GC ((V sh zzz ) 0 purple curve) of a spherical shell using the constant density ρ 0 = 2670 kg m −3 with the influence of the height h ∈ [0, 300] km with an interval of 1 km Table 11 Statistic information of the reference values in log 10 scale of the GP (V sh ) 0 , radial GV (V sh z ) 0 , radial-radial GGT (V sh zz ) 0 , and radialradial-radial GC (V sh zzz ) 0 of a spherical shell with constant density ρ 0 in Fig Table 10. The computation point P is located above the spherical shell (i.e., r > r 2 ) with the height h ∈ [0, 300] km using an interval of 1 km. The relation between geocentric distance r and height h of the computation point P is r = R + h, where R is the radius of the reference sphere and equal to the outer radius of the spherical shell r 2 in Table 10. In other words, the spherical case is adopted. Substituting the above numerical values into the first expressions with the constant density ρ 0 when the computation point is located above the spherical shell in Tables 2, 3 4 and 5, the values of the GP (V sh ) 0 , radial GV (V sh z ) 0 , radialradial GGT (V sh zz ) 0 , and radial-radial-radial GC (V sh zzz ) 0 of a spherical shell using the constant density ρ 0 with the influence of the height are obtained and illustrated in log 10 scale in Fig. 3. In Fig. 3, the x-axis means the height h of the com- Table 12 Mean values in log 10 scale of gravitational potential (V sh ) 1 , radial gravity vector (V sh z ) 1 , radial-radial gravity gradient tensor (V sh zz ) 1 , and radial-radial-radial gravitational curvatures (V sh zzz ) 1 of a spherical shell with linear density coefficient ρ 1 = ρ 0 × 10 −m (m ∈ [2, 4, 6, 8]), where the unit of ρ 1 is kg m −4 and the values are truncated to one decimal place Quantity ρ 0 × 10 −2 ρ 0 × 10 −4 ρ 0 × 10 −6 ρ 0 × 10 −8 (V sh ) 1 8.9 6 .9 4 .9 2 .9 (V sh z ) 1 2.  Table 13 Mean values in log 10 scale of gravitational potential (V sh ) 2 , radial gravity vector (V sh z ) 2 , radial-radial gravity gradient tensor (V sh zz ) 2 , and radial-radial-radial gravitational curvatures (V sh zzz ) 2 of a spherical shell with quadratic density coefficient ρ 2 = ρ 0 × 10 −m (m ∈ [10, 12, 14, 16]), where the unit of ρ 2 is kg m −5 and the values are truncated to one decimal place (V sh zz ) 2 −5.6 −7.6 −9.6 −11.6 (V sh zzz ) 2 −11.9 −13.9 −15.9 −17.9 putation point above the reference sphere. The y-axis means the reference values in log 10 scale of the GP ((V sh ) 0 ), radial GV ((V sh z ) 0 ), radial-radial GGT ((V sh zz ) 0 ), and radial-radial-X.-L. Deng Table 14 Mean values in log 10 scale of gravitational potential (V sh ) 3 , radial gravity vector (V sh z ) 3 , radial-radial gravity gradient tensor (V sh zz ) 3 , and radial-radial-radial gravitational curvatures (V sh zzz ) 3 of a spherical shell with cubic density coefficient ρ 3 = ρ 0 × 10 −m (m ∈ [18, 20, 22, 24]), where the unit of ρ 3 is kg m −6 and the values are truncated to one decimal place

Table 15
Mean values in log 10 scale of gravitational potential (V sh ) 4 , radial gravity vector (V sh z ) 4 , radial-radial gravity gradient tensor (V sh zz ) 4 , and radial-radial-radial gravitational curvatures (V sh zzz ) 4 of a spherical shell with quartic density coefficient ρ 4 = ρ 0 × 10 −m (m ∈ [26, 28, 30, 32]), where the unit of ρ 4 is kg m −7 and the values are truncated to one decimal place radial GC ((V sh zzz ) 0 ) of a spherical shell using the constant density ρ 0 = 2670 kg m −3 . Meanwhile, the related statistic information (i.e., minimum, maximum, mean, and variance) of gravitational effects is shown in Table 11, where the units of (V sh ) 0 , (V sh z ) 0 , (V sh zz ) 0 , and (V sh zzz ) 0 are m 2 s −2 , m s −2 , s −2 , and m −1 s −2 , respectively. From Fig. 3, with the increased heights, the reference values of the (V sh ) 0 , (V sh z ) 0 , (V sh zz ) 0 , and (V sh zzz ) 0 of a spherical shell with the constant density are gradually getting smaller. The statistic information in Table 11 reveals that the ranges of the values of the (V sh ) 0 , (V sh z ) 0 , (V sh zz ) 0 , and (V sh zzz ) 0 in log 10 scale are small especially with the small variances. Consequently, the magnitude of the numerical results for the spherical shell in Fig. 3 is stable, and they can be used as the reference values for the calculated values of the tesseroids in the following numerical experiments.
Furthermore, the mean values in Table 11 using the constant density ρ 0 are set as the reference values. In other words, by changing the density coefficients ρ 1 , ρ 2 , ρ 3 , and ρ 4 at the level of ρ 0 × 10 −m (m ≥ 1), the values of different gravitational effects of a spherical shell can be obtained to reach the same level as the mean values using the constant density ρ 0 in Table 11.
Herein, the relation between ρ n (n ≥ 0) and ρ 0 is theoretically derived. The initial assumption is that the nonzero gravitational effects (i.e., and (V sh zzz ) n ) of a spherical shell with different nth-order polynomial density parts in Tables 2,  3, 4, 5, 6 and 7 are equal when the computation point is located above the spherical shell. The physical meaning of this assumption is the uniform equality of the polynomial density parts. It helps to understand the magnitude of each term of the polynomial density. The formulae of this assumption can be expressed as: Numerical relation between ρ n (n = 1, 2, 3, 4) and ρ 0 with the values of r 1 and r 2 in Table 10, where the nonzero gravitational effects of a spherical shell with different nth-order polynomial density parts are equal when the computation point is located above the spherical shell and the values are truncated to one decimal place Quantity ρ 0 ×Coefficient ρ 1 1.6 × 10 −7 ρ 2 2.5 × 10 −14 ρ 3 3.9 × 10 −21 ρ 4 6.1 × 10 −28 After simplifying Eqs. (40)(41)(42)(43)(44)(45), the relation between ρ n (n ≥ 0) and ρ 0 is obtained as: where this same relation exists for the mass of a spherical shell with different nth-order polynomial density parts if M sh 0 = M sh n in Table 8. Regarding Eq. (46), once the values of inner and outer radii (i.e., r 1 and r 2 ) of the spherical shell are known, the numerical relation between ρ n and ρ 0 can be obtained. Substituting the values of r 1 and r 2 in Table 10 into Eq. (46), the relation between ρ n (n = 1, 2, 3, 4) and ρ 0 is listed in Table 16. The coefficients in Table 16 are in the ranges of [ρ 0 × 10 −7 , ρ 0 × 10 −6 ], [ρ 0 × 10 −14 , ρ 0 × 10 −13 ], [ρ 0 × 10 −21 , ρ 0 × 10 −20 ], and [ρ 0 × 10 −28 , ρ 0 × 10 −27 ] for density coefficients ρ 1 , ρ 2 , ρ 3 , and ρ 4 presented in Tables 12, 13, 14 and 15. Regarding the following experiments in Sects. 3.3 and 3.4, the numerical density coefficients are set according to the formula in Eq. (46) in quadruple precision to show the extreme situation that the values of gravitational effects are correspondingly equal due to the different polynomial density parts.

Influence of height on the GC with polynomial density up to fourth-order
The near-zone problem (or very-near-area problem) of the computation point for the tesseroid mass body (i.e., there are large approximation errors when the computation point approaches the surface of the tesseroid mass body) has been investigated by using tesseroids to discretize the spherical shell for the GP, GV, GGT (Uieda et al. 2016;Shen and Deng 2016), GC Shen 2018a, b, 2019), and invariants of the gravity gradient tensor and their first-order derivatives (Deng et al. 2021) with the constant density. Based on the same numerical situation to investigate the influence of the geocentric distance on the GC (Deng and Shen 2018a, Sect.3.3), we extend the constant density to the polynomial density up to fourth-order to reveal the density variation on the effects of the GC using the relative and absolute errors. The GP, GV, and GGT are included for comparison as well.
For the position of the computation point, the height h is chosen as h ∈ [0, 300] km (Lin et al. 2020, Sect. 3.3) with an interval of 1 km and the spherical longitude and latitude are (0 • , 0 • ) (Deng and Shen 2018a, Sect. 3.3). Other numerical values of the spherical shell and tesseroid can be referred to in Table 10.
The reason for the same performance of relative errors in log 10 scale from zero-to fourth-order polynomial density for different gravitational effects using tesseroids to discretize the spherical shell lies in that the assumption of the density coefficient in Eq. (46) of Sect. 3.2. In other words, the same coefficients are cancelled both in the numerator and denominator in the relative errors for different order polynomial density, whereas the absolute errors can reflect the error superposition effect caused by different order polynomial density in Eqs. (34)(35), which are revealed under the log 10 operation in Fig. 4b.
Furthermore, the existence of the rapid drop zone, especially below the height of 10 km for the GGT and GC in Fig. 4a, shows that the near-zone problem still occurs for different order polynomial density, i.e., the choice of polynomial density model does not affect the existence of the near-zone problem. Solving the near-zone problem will be focused on improving the numerical method for triple (or double) integrals and refined geometrical shape (or grid size) of the tesseroid mass body.

Influence of latitude on the GC with polynomial density up to fourth order
The influence of the computation point's latitude on the GC of the tesseroid mass body was investigated with the constant density using the spherical integral kernels (Deng and Shen 2018a) and Cartesian integral kernels (Deng and Shen 2018b). Herein, the constant density is extended up to fourth-order polynomial density to study the density variation on the effects of the GC with the influence of latitude using the relative and absolute errors. Meanwhile, the low-order gravitational potential gradients (i.e., GP, GV, and GGT) are also investigated. The latitude variation of the computation point indirectly reveals the effects of the geometrical shape of the tesseroid mass body on the gravitational effects.
The computation point is at the satellite height of h = 260 km above the surface of the spherical shell, where the near-zone effects of different gravitational effects are largely reduced as revealed in Fig. 4. The longitude and latitude of the computation point are λ = 0 • and ϕ ∈ [0 • , 90 • ] with an interval of 1 • because of the spherical symmetry of the spherical shell. Table 10 lists other numerical values of the spherical shell and tesseroid.
From Fig. 5a, all the relative errors show the similar rule that as the latitude increases, the relative errors increase as well, especially in the polar region (i.e., ∼ 80 Although the curves between Fig. 5a and Fig. 8 of Deng and Shen (2018a) are similar, it should be noted that there are some differences between them. Firstly, the coordinate system of gravitational effects in Fig. 5a is the local northoriented frame system with x, y, z pointing to the north, east, and radial up and in Fig. 8 of Deng and Shen (2018a) is the local East-North-Up topocentric reference system with x, y, z pointing to the east, north, and radial up. In other words, (δV x x ) 4 , (δV yy ) 4 , (δV x xz ) 4 , and (δV yyz ) 4 in Fig. 5a are correspondingly with respect to GGT_δV yy , GGT_δV x x , GC_δV yyz , and GC_δV x xz in Fig. 8 of Deng and Shen (2018a). Secondly, there are no gaps for (δV yy ) 4 and (δV yyz ) 4 at latitude ϕ = 90 • in Fig. 5a, whereas the curves of the GGT _δV x x and GC_δV x xz have gaps at latitude ϕ = 90 • in Fig. 8 of Deng and Shen (2018a). The reason is that the tesseroid formulae of the GGT and GC in this paper are in Cartesian integral kernels, whereas they were in spherical integral kernels in Deng and Shen (2018a). Using Cartesian integral kernels can avoid the gap at the polar point (i.e., the polar-singularity problem). Thirdly, the curve of the (δV yyz ) 4 is smooth at the low latitude region (i.e., ϕ < 44 • ) in Fig. 5a, whereas the GC_δV x xz has a slightly rough curve at the same low latitude region. Fourthly, the numerical methods are Gauss-Legendre quadrature with 3D degrees (2, 2, 2) for Fig. 5 in this section and Taylor series expansion with a second-order for Fig. 8 of Deng and Shen (2018a). Finally, Fig. 5 in this section adopts the fourth-order polynomial density model, whereas Fig. 8 of Deng and Shen (2018a) applied the constant density.
X.-L. Deng   Fig. 5 a Visualization of the relative errors in log 10 scale of the GP ((δV ) 4 blue curve), GV ((δV z ) 4 red curve), GGT ((δV xx ) 4 dark-blue curve, (δV yy ) 4 green curve, and (δV zz ) 4 yellow curve), and GC ((δV xxz ) 4 thistle curve, (δV yyz ) 4 deep-sky-blue curve, and (δV zzz ) 4 magenta curve) with fourth-order polynomial density by a grid size of 15 × 15 with the influence of the latitude ϕ from 0 to 90 • with an interval of 1 • at the satellite height of h = 260 km; b the absolute errors in log 10 scale of the Laplace parameters of the GGT ((δΔL 1 ) N ) and GC ((δΔL 2 ) N ) using different order polynomial density (i.e., N = 0, 1, 2, 3, 4), where other parameters are the same as in Fig. 4b Regarding Fig. 5b, the ten curves of the absolute errors in log 10 scale do not change particularly as the latitude rises. Similar superposition effects of absolute errors both for Laplace parameters of the GGT ((δΔL 1 ) N ) and GC ((δΔL 2 ) N ) are obvious as the increased polynomial density from zero-to fourth-order. The absolute errors in log 10 scale are in ranges of about [−42, −39] for (δΔL 1 ) N and [−46, −43] for (δΔL 2 ) N , which are all below the machine epsilon of quadruple precision in log 10 scale as −34. It can be confirmed that under the current numerical condition, the components of the GGT and GC with zero-up to fourthorder polynomial density in Eqs. (34-35) satisfy Laplace's equation at the satellite height of h = 260 km.

Conclusions and outlooks
In recent years, the research trends in gravity field modeling have been from constant to variable density and from low-order to higher-order gravitational potential gradients (e.g., the GC). In this contribution, the formulae of the GC of a tesseroid and spherical shell are extended from constant density to arbitrary-order polynomial density. In detail, the general expression for the GC of a tesseroid with N thorder polynomial density is derived in Cartesian integral kernel, and the detailed components of the GC are presented in the spherical coordinate system. The general analytical expressions for gravitational effects up to the GC of a spherical shell with N th-order polynomial density are derived when the computation point is located above, inside, and below the spherical shell. Meanwhile, the general formula of the mass for a spherical shell having a polynomial density model up to N th-order is derived. Furthermore, we derive the general relation between radial gravitational effects up to arbitrary-order and the mass of a spherical shell with N thorder polynomial density with the computation point located above the spherical shell.
Under the assumption that the nonzero gravitational effects of a spherical shell with different nth-order polynomial density parts are equal when the computation point is located above the spherical shell, the effects of density values on gravitational effects up to the GC of a spherical shell are investigated. The relation between the constant density ρ 0 and density coefficient ρ n (n ≥ 0) is derived. The difference of the magnitude between the density coefficients ρ n and ρ n+1 is about at the level of −7.
The influence of the computation point's height on gravitational effects up to the GC with polynomial density up to fourth-order is investigated. Numerical experiments show that the near-zone problem occurs for the GC with different order polynomial density. In other words, the change in density does not affect the existence of the near-zone problem. The relative errors in log 10 scale of the GGT and GC are large than 0 below the height of about 24 km and 50 km, respectively. The key to solving the near-zone problem of the GC of the tesseroid lies in the improvement of the numerical algorithm to calculate the triple or double integrals and the selection of the geometrical shape of the tesseroid mass body, e.g., the rotation method (Marotta and Barzaghi 2017;Marotta et al. 2019), splitting line method using the double exponential quadrature (Fukushima 2018a), different types of the regular, adaptive and combined subdivision (Li et al. 2011;Grombein et al. 2013;Shen and Deng 2016;Uieda et al. 2016;Deng and Shen 2019;Lin and Denker 2019;Soler et al. 2019;Zhong et al. 2019;Zhao et al. 2019;Lin et al. 2020;Chen 2020, 2021).
In addition, we study the influence of the latitude on gravitational effects up to the GC with the polynomial density up to fourth-order at the satellite height of h = 260 km. The polar-singularity problem does not occur for the GC and GGT with different order polynomial density because of the applied Cartesian integral kernels of the tesseroid. The relative errors in log 10 scale of gravitational effects up to the GC increase with the increased latitude. In other words, the geometrical shape of the tesseroid in the high latitude region, especially in the polar region, leads to an increase in the relative errors, whereas the relative errors in log 10 scale are still within an acceptable range below −3. Under the assumption in Eq. (46), the density variation can be revealed in the superposition effects of the absolute errors of Laplace parameters of the GGT and GC in log 10 scale with the influence of the height and latitude.
Regarding the potential applications of the formulae of the GC of a tesseroid in Cartesian integral kernels with arbitraryorder polynomial density in Eqs. (9-10), the atmospheric, topographic, and crustal effects of the Earth and other celestial bodies can be studied for the GC components in the gravity field modeling, which are similar to the GGT components (Eshagh 2009a(Eshagh , b, c, 2010(Eshagh , 2021. For example, the contribution of the binomial expansion up to degree four in topographic and atmospheric effects for the GC at the satellite altitude will be investigated in the future based on the similar application of the contribution of the second and third terms of the topographic and atmospheric effects for the GGT (Eshagh 2009a). In the future, the chosen exponents of the polynomial density coefficient of the GC will be examined for their topographic and atmospheric effects.
The analytical expressions for the gravitational effects up to the GC of a spherical shell with arbitrary-order polynomial density in Eqs. (19) and (25) can be regarded as the reference values for the tesseroid and other spherical mass bodies (e.g., spherical triangular tessellation (Zhang et al. 2018a)) in the numerical experiments. Other variable density models (e.g., exponential and sinusoidal density functions (Soler et al. 2019)) for the GC of the tesseroid and spherical shell will be focused on compared to the arbitrary order poly-nomial density model. Moreover, the residual topographic effects of the Earth can be investigated for the GC components using the spherical shell in the gravity field modeling, which are analogous to the GGT components . Note that the spherical harmonic expansion of the potential for the spherical shell even with irregular boundaries and laterally variable density was determined in the spectral domain (Eshagh 2010). The spherical harmonic coefficients of the potential of the irregular and heterogeneous layers (e.g., topographic and atmospheric effects in Eshagh (2009b); crustal effects in Eshagh (2021)) can simply be computed and inserted into the spherical harmonic expansions of the GC. Further investigations on the consistencies of the GC of the spherical shell between this paper in spatial domain and Eshagh (2009bEshagh ( , 2010Eshagh ( , 2021 in spectral domain will be performed in the future, c.f. Kuhn and Seitz (2005), Hirt and Kuhn (2014), and Hirt et al. (2016).
When modeling the GC with an arbitrary-order polynomial density of the atmosphere, water, ice, topography, crust, mantle, and core of the Earth or other celestial bodies, the choice of density coefficient in the polynomial density in the practical situation should be made carefully according to the available density models (e.g., CRUST1.0 (Laske et al. 2013), GEMMA (Reguzzoni and Sampietro 2015), and UNB_TopoDens (Sheng et al. 2019)), seismic waves, and geological survey results. The mentioned density models include the density jumps at the boundaries of the mass bodies. The technique to treat the density jumps can be referred to Eqs. (26)(27). In other words, the computation point can be slightly moved inside the mass bodies to avoid the density jumps based on the chosen double or quadruple precision. When adopting these density models in the application of the topographic effects, the density gradients can be taken, and the grid sizes of the tesseroids are 1 • × 1 • for the CRUST1.0, 0.5 • × 0.5 • for the GEMMA, and 30 × 30 for the UNB_TopoDens.
The modeling of the GC of the tesseroid and spherical shell with different layers of the Earth or other celestial bodies belongs to the forward problem in geophysics. The low-order gravitational quantities (e.g., GV and GGT) of the tesseroid have been widely applied for the inversion in studies of the internal structures of the Earth or the Moon (Liang et al. 2014;Uieda and Barbosa 2017;Zhang et al. 2018b;Zhao et al. 2019Zhao et al. , 2021. In the future, the inverse problem of the GC using different mass bodies (e.g., the tesseroid and spherical shell) that are likely to appear in applications will be investigated.
Author Contributions XLD derived the formulae, performed the numerical experiments, wrote the manuscript, and approved the final submitted version of the manuscript.
Funding Open Access funding enabled and organized by Projekt DEAL.

Data availability
The data that support the findings of this study are available from the corresponding author on reasonable request.

Conflict of interest
The author declares no competing financial interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

Appendix 1: Derivation of analytical expressions for the GP of a spherical shell with arbitrary-order polynomial density
The analytical expressions for the GP (V (r )) of a spherical shell with arbitrary-order polynomial density are derived with the computation point P located above (i.e., r > r 2 ) inside (i.e., r 1 ≤ r ≤ r 2 ), and below (i.e., r < r 1 ) the spherical shell.
If the computation point P is located above the spherical shell (i.e., r > r 2 ), 1 = r − r and 2 = r + r . Then, Eq. (15) can be presented as: When the computation point P is located below the spherical shell (i.e., r < r 1 ), 1 = r − r and 2 = r + r . Under this condition, Eq. (15) can be derived as: If the computation point P is located inside the spherical shell (i.e., r 1 ≤ r ≤ r 2 ), the original spherical shell can be divided as two spherical shell (Lin et al. 2020), i.e., the first spherical shell with inner and outer radii as [r 1 , r ] where the computation point is above this spherical shell (i.e., 1 = r − r and 2 = r + r ) and the second spherical shell with inner and outer radii as [r , r 2 ] where the computation point is below the spherical shell (i.e., 1 = r −r and 2 = r +r ). In this case, Eq. (15) can be derived as: