External gravitational field of a homogeneous ellipsoidal shell: a reference for testing gravity modelling software

There are numerous applications in geodesy and other geo-sciences in which the gravitational potential effect or other functions of the potential are computed by forward modelling from a given mass distribution. Different volume discretisations, e.g. prisms, tesseroids or mass layers are used. In order to control the numerical realisation of the forward calculation in the practical application, e.g. in reduction tasks, these evaluation programs should be verified against rigorous analytical solutions. In this contribution, a closed analytical solution for the potential of an ellipsoidal shell as a test body is presented. Furthermore, we derive the respective closed formulae for the gravity vector and the gravity gradient tensor. Program implementations of the tesseroid approach are compared on the basis of this ellipsoidal mass arrangement. For the practical usage, fast-converging expansions in spherical harmonics are provided in addition. The derivation of the formulae is based on a closed solution of the potential of a homogeneous ellipsoid for computation points situated on the rotation axis, which then is extended to the external space.


Introduction
The determination of the gravity field of the Earth and its variation in time is one of the main tasks of Geodesy (Rummel et al. 2005). The most dominant spatial variations on the Earth's surface and in external space close to the Earth are induced by the topographical, isostatic masses and hydrological mass variations. The modelling of these masses is generally performed by volume discretisations, reducing the shape of the volume elements to geometrically simple bodies such as rectangular prisms (Forsberg 1984;Nagy et al. 2000Nagy et al. , 2002, polyhedrons (D'Urso 2014), spherical tesseroids (Anderson 1976;Heck and Seitz 2007 Marotta et al. 2019) or ellipsoidal prisms (Abd-Elmotaal and Kühtreiber 2021; Roussel et al. 2015). The precise computation of the gravitational field of the topographic and isostatic masses is still a numerical challenge, despite of the use of modern parallel computers with distributed kernels, since-due to the irregularity of these mass constituents-a huge amount of volume elements has to be used in the discretisation of the attracting masses. This applies in particular due to the continuing increase in the resolution of volume elements.
Efficient approaches have been developed in the past years for the accurate calculation of the gravitational potential and its functionals and have been implemented by suitable software. In order to control those approaches and the respective computer programs, the numerical results should be compared with and verified against analytical solutions, based on simple mass distributions. The most simple global mass distribution, approximating the Earth's topography, is a spherical shell, which has often been applied in Geodesy for testing evaluation procedures and software (cf. Heck and Seitz 2007;Fukushima 2018). A massive spherical shell bounded by two concentric spheres produces a gravitational field that can be described by a closed analytical solution. However, a spherical shell creates a mathematically sim-ple, isotropic gravitational field. Only the potential and its radial derivatives do not vanish and can therefore be tested in such a simple test embedding. By using Laplace's equation, V xx = V yy = −V zz /2 can be calculated indirectly in the isotropic case.
In the present paper, we propose an ellipsoidal massive shell bounded by two confocal ellipsoids of revolution as a more general and realistic test configuration. Using the ellipsoid instead of the isotropic spherical case leads to a much better approximation of the general shape of the Earth as it depends also on the horizontal position of the computation point. In addition to the derivatives in radial direction also the latitudinal and mixed gradients can be controlled by this ellipsoidal setting. Analytical solutions for the external gravitational potential of a homogeneous ellipsoid of revolution have been derived by various authors, e.g. by Kellogg (1929), Hobson (1931), MacMillan (1958) and Wang (1988). Any of these solutions refers to coordinate systems which are rather unfamiliar in Geodesy. The arctan, sinh, cosh functions used in these approaches are also very time-consuming in terms of numerical efforts.
In this paper, we derive solutions in spherical coordinates and elliptical coordinates which have not been published before, according to our knowledge. Furthermore, we follow a new innovative way of derivation: after having obtained a closed solution for the special case of a computation point on the positive z-axis, this result is extended to the complete external space of the ellipsoid by harmonic continuation. As a by-product, this procedure provides an illustration of the well-known Runge-Krarup theorem related to the validity of the spherical harmonic series. We prove that any homogeneous ellipsoid of revolution with the same total mass M (but different in density and volume) and the same linear eccentricity (but different semi-major axes) produces the same gravitational field at any computation point in external space.
While the solution in terms of spherical coordinates produces an infinite, but very fast-converging series in spherical harmonics, an alternative derivation in elliptical coordinates results in a finite series of ellipsoidal harmonics of degrees zero and two. However, in practice the series solution in spherical harmonics will be preferred due to its higher numerical effectivity. The computing times required in both cases are comparatively small. A significant advantage of the representation in spherical harmonics series expansion is the possibility of obtaining their spectral components and being able to use them for further investigations.
The structure of the paper is as follows: in Sect. 2, the geometrical properties of an ellipsoidal shell are explained, while the central Sect. 3 presents two approaches for the calculation of the gravitational potential of a solid homogeneous ellipsoid of revolution and an ellipsoidal shell, based on spherical and ellipsoidal harmonics. A new methodology for the analytical derivation of the mathematical expressions is presented.
Section 4 contains numerical comparisons between several solution approaches, before a discussion of results and some conclusions finalise the paper. In an appendix, we provide the complete expressions for the gravity vector and the gravity gradient tensor (first-and second-order derivatives of the potential) in terms of spherical and elliptical coordinates, and investigate the convergence behaviour of the spherical harmonic series.

Geometry of an ellipsoidal shell
It is assumed that the ellipsoidal shell is situated between two concentric, confocal ellipsoids ε 1 with semi-axes a 1 , b 1 and ε 2 with semi-axes a 2 , b 2 having the same rotational axis (cf. Fig. 1). The attached Cartesian coordinate system is defined as follows: the z-axis is the rotation axis and the (x, y)−plane is the equatorial plane. For the dimensions of the two considered ellipsoids the relation a 1 > a 2 and b 1 > b 2 holds, this means that ε 2 is completely situated inside ε 1 . External field means that the calculation point P is situated on or outside the ellipsoid ε 1 .
Since both ellipsoids have been assumed as confocal, the linear eccentricity E = √ a 2 − b 2 is the same for ε 1 and ε 2 . Numerical values for the parameters a i , b i and e i = E/a i of various confocal ellipsoids are listed in Table 1, starting from the respective values of the Geodetic Reference System 1980 (Moritz 1980b), where E GRS80 = 521854.0097 m.

Gravitational potential of a homogeneous ellipsoid and ellipsoidal shell
In this chapter, we will present two approaches for the calculation of the potential and its derivatives at points outside   the ellipsoid. The first approach is based on representations in terms of spherical harmonics. The second approach is related to ellipsoidal harmonics. However, both representations describe the same gravitational field and lead to the same numerical values (besides rounding effects) for the potential and its derivatives outside the ellipsoid ε 1 . The derivation of the formulae relies on the following procedure: first, we derive a closed analytical solution for the special case of a computation point situated on the positive z-axis, the rotation axis of the ellipsoid. In a second step, this result is extended to computation points in the complete external space, using harmonic continuation by representations in terms of spherical and ellipsoidal harmonics. As a by-product, an illustration of the well-known Runge-Krarup theorem related to the validity of the spherical harmonic series in the spatial domain below the Brillouin sphere is provided.
The gravitational potential of a homogeneous ellipsoidal shell results from the difference of the external potential fields of the two confocal ellipsoids ε 1 and ε 2 , due to the linear superposition principle of potentials.

Approach I: analytical expansion in spherical harmonics
The gravitational potential of an ellipsoidal shell with constant mass density ρ can be expressed by Newton's integral (Heiskanen and Moritz, 1967, (1-11)) at the computation point P outside the domain containing the total mass M. P ∈ is the running integration point. Newton's gravitational constant is denoted by G. The Euclidean distance between P and P can be expressed in Cartesian coordinates (x, y, z) or in geocentric spherical coordinates (r , ϑ, λ) or in cylindrical coordinates (η, ζ , λ) The respective differential volume elements d in the mentioned coordinate systems, which are used within this paper, are d = dx dy dz , The spherical distance ψ between the position vectors of P and P which occurs in Eq.
3.1.1 Special case: gravitational potential of a homogeneous ellipsoid at P ∈ z-axis First, let us consider a computation point P on the z-axis, with z ≥ b and x P = y P = 0. It is tempting to use spherical polar coordinates (r , ϑ, λ) for the evaluation of the Newton integral (1), but this procedure leads to insurmountable difficulties in the integration process. As an alternative we prefer cylindrical coordinates (η , ζ , λ ) for the position of the integration point P , where λ is the geographical longitude, η the distance from the rotation axis of the ellipsoid and ζ the distance from the equatorial plane (see Fig. 2). The distance from the computation point to the integration point can be expressed in cylindrical coordinates by setting ζ P ≥ max(ζ ), η = 0 in Eq. (4) Expressing Newton's integral (1) by replacing the volume element d in cylindrical coordinates (see Eq. (5)) results in the explicit form of Newton's integral The upper limit of η is the orthogonal distance of the ellipsoidal surface points from the z-axis, depending on ζ : p(ζ ), which can be derived from the equation of the ellipsoid. From (Heck 2003) results Integration with respect to λ gives simply while the integral with respect to η yields where is the square of the second numerical eccentricity. In the following, we also refer to the first numerical eccentricity e and the linear eccentricity E with related by the well-known formulae (1 − e 2 )(1 + e 2 ) = 1 , where X = αx 2 + βx + γ and α < 0, yields The difference between two arcsin terms, which occur in (18), can be expressed by Together with yields after some tedious, but elementary manipulations Defining x := ae z and using the relationship Finally, let us consider the total mass M of the homogeneous ellipsoid resulting in the closed expression For an evaluation point P outside the ellipsoid which is located on the (positive) z-axis, z > b holds. It is straightforward that E/z = be /z < e = 0.082 094 438 (Moritz 1980b). This leads to (E/z) 2n+1 < 2 · 10 −8 even for n = 3 which indicates the very fast numerical convergence.

General case: gravitational potential at arbitrary points outside the ellipsoid
Expression (23) holds for any computation point on the positive z-axis (or on the negative z-axis if z is changed to |z|) for z ≥ b. This equation can easily be extended to the whole space outside the ellipsoid, considering that V in the external space is a harmonic function, with rotational symmetry with respect to the z-axis and equatorial symmetry with respect to the (x, y)−plane: with unknown coefficients c 2n . P 2 and P 4 denote Legendre polynomials of degrees 2 and 4, depending on the polar distance ϑ (complement to geocentric latitude). Specifying this expression for points on the positive z-axis, i.e. r = z, ϑ = 0 yields which can be compared with the series expansion (23). Termwise comparison results in the coefficients c 2n Thus, we obtain the final result for the gravitational potential of the homogeneous ellipsoid at any computation point outside the ellipsoid The terms of degree zero and two (n = 0 and n = 1) are consistent with the respective expressions by Moritz (1990, p. 36), approximately derived in a different context.  (27) that the potential of a homogeneous ellipsoid only depends on its total mass M and linear eccentricity E, not on any other ellipsoidal parameter. Thus, any homogeneous ellipsoid of revolution with the same mass M and the same linear eccentricity E (confocal ellipsoids) produces the same gravitational field at any point P(r , ϑ, λ) in the external (mass-free) space, i.e. outside of ε 1 (see Fig. 1).
Usually, the spherical harmonic expansion of the gravitational potential is scaled to a fixed radius R > E (e.g. R = a, Brillouin sphere, or R = b, Bjerhammar sphere). For a mass distribution with rotational and equatorial symmetry, the standard expression is provided by the formula Since Eqs. (27) and (28) represent the same function, the Stokes' coefficients V n in Eq. (28) can be analytically expressed as or in fully normalised coefficients Numerical values for the fully normalised coefficients derived in Eq. (30) are listed in Table 2. They are based on the linear eccentricity E G RS80 = 521854.0097 m and semi-major axis a = R = 6378137 m of the Geodetic Reference System 1980 (GRS80) (Moritz 1980b).
Remark 2 The construction of Eqs. (27) and (28) by the extension of the closed expression Eq. (22) from the z-axis to the external space can be considered as an illustration of the validity of the Runge property of harmonic functions, which in the present context is formulated by the Runge-Krarup theorem (Moritz, 1980a, p. 67 ff). Loosely spoken the Runge-Krarup theorem, applied to our case, guarantees that any harmonic function, regular outside the ellipsoid, can be approximated arbitrarily well by harmonic functions regular outside a sphere situated completely within the ellipsoid. Due to our constructive approach based on Eq. (24), we have extended the domain of harmonicity from the region outside a Brillouin sphere, enclosing the ellipsoid completely, to the larger region outside a Bjerhammar sphere, situated completely within the ellipsoid; as Eq. (27) shows, the radius of this internal sphere may be as small as R = E, the linear eccentricity of the ellipsoid. As a matter of fact, although the spherical harmonic series (27) and (28) converge in the domain outside the focal sphere R = E, they do not represent the potential inside the ellipsoid, which is not a harmonic function but fulfils the Poisson equation inside the masses. In contrast, the spherical harmonic series (27) and (28) represent the analytical (harmonic) continuation in the internal domain.
Based on Eq. (27), it is straightforward to develop the first and second derivatives of the potential at a point in the space outside the ellipsoid.
The first and second radial derivatives of the potential are A general representation of the radial derivative of any degree k > 0 is provided in Appendix A.3. Applying the relation between the partial derivatives in spherical coordinates and Cartesian coordinates in a local topocentric system which are given in Tscherning (1976) or Grombein et al. (2013), we elaborated the components of the gravity vector and the gravity gradient tensor g i and M i j of the potential of a homogeneous ellipsoid of revolution. The expressions are compiled in Appendix A.1.
The gravitational potential and its derivatives at computation points on the surface of the ellipsoid can be calculated by inserting r E (ϑ) for r in Eqs. (27), (31) and (32) (cf. Heck (2003))

Gravitational field of a homogeneous ellipsoidal shell
In the following, we will derive the gravity field of a homogeneous shell between two confocal ellipsoids of revolution. Confocal ellipsoids ε 1 and ε 2 have the same linear eccentricity (see Fig. 1 and compare the definition in Eq. (15) From Eqs. (27), (31) and (32), it is obvious that just the difference of the masses M 1 and M 2 has to be considered: With this premise, the final expression for the potential of a homogeneous ellipsoidal shell results in It should be noted that the radial distance between two confocal ellipsoids is not constant. Figure 3 presents the variation of the thickness of the ellipsoidal shell from the pole to the equator, where a 1 corresponds to the semi-major axis of the GRS80 ellipsoid and a 2 = a 1 − 10 km. The maximum difference occurs at the pole with d(ϑ = 0 • ) = 34 m; see also Table 1.
The lower bound r LB (ϑ) of the shell is the ellipsoidal radius of the confocal ellipsoid ε 2 and the upper bound r E (ϑ) is defined in Eq. (33). Throughout this chapter, we fix the dimension of ε 2 by definition of the thickness of the ellipsoidal shell at the equator with d = a 1 − a 2 > 0. The latitude-dependent thickness of the ellipsoidal shell is sufficiently accurately approximated by Elliptical (u, θ, λ) and spherical (r , ϑ, λ) coordinates, the dimensions of confocal ellipsoids as well as the thickness of the ellipsoidal shell are depicted in Fig. 4.
The semi-major and semi-minor axis of the confocal ellipsoid ε 1 is denoted by a and b, and the semi-minor axis of any smaller confocal ellipsoid ε 2 is denoted by the elliptic coordinate u. The thickness d(ϑ) of the ellipsoidal layer can be computed from the difference of Eqs. (36) and (33). Its slight increase from the equator to the pole is visible in Fig. 3. At the equator d(ϑ = π/2) = d holds. Figure 5 shows the contribution of each degree n to the potential, assuming a 1 −a 2 = 10 km. The rapid convergence of the series (35) is obvious. From the extremal numerical example (thickness of the shell is approx. 10 km, P on ε 1 ), Fig. 5 Contributors to the series expansion for the gravitational potential of an ellipsoidal homogeneous shell, calculated on the surface of ε 1 along a meridian. a 1 − a 2 = 10 km, ρ = 2670 kg m −3 the truncation error is less than 10 −5 m 2 s −2 if the series is expanded only up to degree n = 6.
The radial derivative of the potential of a homogeneous shell between two confocal ellipsoids of revolution results from Eqs. (31) and (34) by the findings that the potential of the shell is only a function of the mass difference between the two confocal ellipsoids: In Fig. 6, the rapid convergence of the series (38) is obvious. From the extremal numerical example (thickness of the shell is approx. 10 km, P on ε 1 ), the truncation error is less than 10 −6 mGal if the series is expanded only up to degree n = 6. Even if the series is truncated after n = 4, the approximation error is still below 10 −3 mGal.

Approach II: expansion in ellipsoidal harmonics
As an alternative to spherical harmonics the external gravitational potential of a mass distribution can be represented in ellipsoidal harmonics, which will be based on elliptical coordinates (u, θ, λ). According to Heiskanen and Moritz (1967, p. 43), specified for a body with rotational and equatorial symmetry, which is independent on longitude λ, the potential can be expressed by the infinite series Fig. 6 Contributors to the series expansion for the negative radial derivative of the gravitational potential of an ellipsoidal homogeneous shell, calculated on the surface of ε 1 along a meridian. a 1 −a 2 = 10 km, ρ = 2670 kg m −3 whereι = √ −1 is the imaginary unit and E is the linear eccentricity as defined in (15); a 2n are complex-valued coefficients. For n = 0 and n = 1, the Legendre functions of second kind Q 2n can be taken from Heiskanen and Moritz (1967, p. 66) For the derivation of the ellipsoidal harmonic coefficients a 2n , we use a procedure analogous to chapter 3.1; i.e. we first consider a computation point P situated on the z-axis and then generalise the result to the complete external space by harmonic extension. For points on the rotational (= z) axis of the homogeneous ellipsoid of revolution θ = 0 • holds. From this, it follows that for all Legendre polynomials P 2n (cos θ) = 1 ∀ n holds. Furthermore, it is obvious that for the points P located on the rotational axis the coordinate z in the polar system of the point P is identical to u; thus, u = z.
Based on this new approach (ansatz), Eqs. (39) and (40) Now, let us compare this expression with the closed form (22). It can easily be recognised that which can be inserted in (39) in order to obtain the development of the external potential in ellipsoidal harmonics: Obviously, the ellipsoidal harmonic series is finite and is composed of degree 0 and 2 terms only. As a result, we get the gravitational potential of the homogeneous ellipsoid of revolution in the space outside the ellipsoid with semi-minor axis b and linear eccentricity E in closed form In Appendix A.2, we present the partial derivatives of first and second order with respect to the elliptical coordinates u and θ and compose them into the components of the gravity vector and gravity gradient tensor with respect to the local triad (−θ, λ, u), corresponding to a North/East/Up coordinate system. It is obvious that-although the formulae in elliptical coordinates are closed analytical expressions-in practical applications the fast-converging, truncated spherical harmonic series are more convenient and will be preferred due to their simpler programming feasibility. In addition, the representation in spherical harmonics provides information about the spectral content, which can be of great importance in practice.
In detail, the components of the gravity vector are Since the linear eccentricity E is the same for all confocal ellipsoids, the gravitational potential of a homogeneous ellipsoidal shell results from Eq. (43) by replacing the mass M by the mass M of the shell, provided by Eq. (34).
On the surface of the homogeneous ellipsoid (u = b), the expression (43) reduces to where e = E b denotes the second numerical eccentricity of the ellipsoid (see Eq. (16)). Similarly, the gravitational potential of an ellipsoidal shell results from Eq. (47) by replacing M by M, when b and e are the parameters of the upper boundary surface containing the computation point.
Throughout this paper, geocentric spherical coordinates are denoted by r , ϑ, λ and elliptical coordinates by u, θ , λ. In elliptical coordinates, a point P in space is fixed by the intersection of a confocal ellipsoid u = const., a hyperboloid θ = const. and the meridian plane λ = const. passing through P. Furthermore, geodetic (geographical) coordinates are denoted by h, ϕ, λ. The longitude λ is identical in any coordinate system used. The "vertical" coordinates are accordingly r , u and the ellipsoidal height h in the case of geodetic coordinates. The geodetic latitude ϕ of a point P on the surface of an ellipsoid of revolution is the angle between the ellipsoidal normal, which is running through P, and the projection of the ellipsoidal normal onto the equatorial plane (see Fig. 4).
Equations (48) and (49) result from the manipulation of (Heiskanen and Moritz, 1967, (1-103)) and (Heck, 2003, (3.35), (3.37)). The geodetic latitude ϕ, the spherical colatitude ϑ and the elliptical co-latitude θ are related to each other by the following equations Remark 3 Equations (48) or (49) can easily be inserted in Eq. (47). Obviously, the gravitational potential on the surface of a homogeneous ellipsoid can be represented in closed analytical form not only in elliptical coordinates but also when geodetic or spherical latitudes are used.

Numerical comparisons
The formulae derived in Sect. 3 provide helpful analytical tools for testing algorithms and software in gravity field modelling. In Geodesy, a spherical shell has often been applied to testing software for calculating topographic gravity field functionals (cf. Heck and Seitz, 2007;Fukushima, 2018). However, a spherical shell creates a mathematically simple gravitational field that is independent from the horizontal position; in contrast, an ellipsoidal shell provides a more geo-realistic test setting which is also dependent on the geographical latitude because it takes the flattening of the Earth into account. As a very concrete practical consequence of modelling the test masses in the form of an ellipsoidal shell, not only the potential and its radial derivatives but also the latitude-dependent (and mixed ones with r ) derivatives are different from zero, so that the software is also controllable for these functionals. The massive ellipsoid or ellipsoidal shell can be discretised in various ways, as performed by different groups, resulting in numerical values for the potential and its derivatives (see, for example, Šprlák et al, 2020). By the derived analytical formulae, we have a means to check the performance of those approaches and software, in space as well as in frequency domain. Using a global discretization of an ellipsoidal topography, we demonstrate in Chapter 4.2 such a benchmark check related to the tesseroid approach according to Heck and Seitz (2007) and Grombein et al. (2013) which has been implemented in software used by different groups. Before, in Chapter 4.1, we validate the spherical harmonic expansion described in Chapter 3.1 with respect to the rigorous closed formula (43).
All numerical investigations are carried out for computation points P i running along a meridian at λ i = ϑ/2 = constant with 0 • < ϑ i < π/2, ϑ i = ϑ/2 + (i − 1) ϑ. A step size of ϑ = 5 is chosen. P i (ϑ i , λ i ) is located on the surface of the ellipsoid ε 1 . In all numerical inves- Fig. 7 Difference in potential between the solution in ellipsoidal and spherical harmonics of an ellipsoidal homogeneous shell, calculated on the surface of ε 1 along a meridian. a 1 − a 2 = 10 km, ρ = 2670 kg m −3 tigations carried out, the Newtonian gravitation constant G = 6.6743 · 10 −11 m 3 kg −1 s −2 is used. The mass density is set to ρ = 2670 kg m −3 without restriction of generality.
Since the potential of an ellipsoidal shell and also its radial derivative are rotationally symmetrical, and there is equatorial symmetry too, this test setting is obvious and sufficiently general.

Spherical versus ellipsoidal harmonic expansion
In this chapter, the analytical representations developed in Sect. 3 for the potential of a homogeneous ellipsoidal shell in terms of spherical and ellipsoidal harmonics are compared numerically.
We compare the results given in Eqs. (30) and (35) with the exact closed solution in two terms based on ellipsoidal harmonics (47). To ensure that the evaluation of the potential is taken at the same points P i on the ellipsoid, the geocentric polar distance ϑ i is converted to the elliptical co-latitude θ i by the use of Eq. (48). Both latitude parameters are visualised in Fig. 4.
In Fig. 7, the difference between the potential of the ellipsoidal shell given in ellipsoidal harmonics and the series solution in spherical harmonics is provided for N ∈ {0, 2, 4, 6, 8, 10} in semi-logarithmic scale.
The approximation error in spherical harmonics compared to the closed (exact) solution in ellipsoidal harmonics is as follows: for the potential, it is less than 10 −5 m 2 s −2 if the series is expanded only up to degree N = 6, as can be seen in Fig. 7. For N ≥ 10, the results are numerically identical indicated by the difference smaller than 10 −9 m 2 s −2 which primarily is due to rounding errors and corresponds to a relative accuracy of 10 −14 .
The difference in gravitational attraction between the solution in ellipsoidal (Eq. (38)) and spherical (Eq. (46)) harmonics of an ellipsoidal homogeneous shell of thickness d = 10 km and homogeneous density ρ = 2670 kg m −3 is less than 10 −3 mGal if the series is expanded only up to degree N = 6. This can be seen in Fig. 8 indicated by the dashed lines. If δg s (N = 20) is used as reference, the approximation error becomes much smaller (cf. Fig. 8 solid lines). From these results, it can be concluded that the closed solution faces some numerical limitations, rather than the solution in spherical harmonics.
Statistical values for the potential and gravitational attraction in ellipsoidal and spherical harmonics are provided in Table 3.
The equivalent representation of the potential and its attraction in spherical and elliptical coordinates (Eqs. (35), (43) and (38), (46), respectively) is confirmed by the identical numerical results (besides negligible numerical effects) presented in Figs. 7 and 8.

Tesseroid discretisation versus spherical harmonic expansion
Usually, a computational implementation of a spatial mass discretisation is validated against an analytical solution.
Here, we apply extended spherical volume elements, called tesseroids. A spherical tesseroid is bounded by two meridians in East-West direction with geocentric longitudes λ 1 and λ 2 . In North-South direction, the tesseroid is bounded by two concentric conical surfaces with opening angles ϑ 1 = π/2 − ϕ 1 and ϑ 2 = π/2 − ϕ 2 . In the spherical case, the top and bottom bounding surfaces are spheres of radii r 1 = R and r 2 = R + d. The spherical tesseroid is depicted in Fig. 9. An ellipsoidal tesseroid is bounded by the confocal ellipsoids ε 1 and ε 2 , respectively, where the distance d between Fig. 9 Geometry of a spherical tesseroid (Heck and Seitz 2007) the bounding surfaces is latitude-dependent (Roussel et al. 2015).
Routines for the computation of the gravitational field of a spherical tesseroid have been developed and tested in many studies (Marotta and Barzaghi 2017;Hirt et al. 2014;Grombein et al. 2013;Uieda et al. 2016). In Heck and Seitz (2007) and Fukushima (2018), the potential and the gravitational attraction (negative radial derivative) was validated against a spherical shell, using the closed formulae (Vaníček et al. (2001(Vaníček et al. ( , 2004) In the spherical case, the thickness of the shell is H = d = r 2 − r 1 . R = r 1 = 6368137 m is the constant radius of the sphere on which the shell is elongated. Already in Heck and Seitz (2007), it is pointed out that the spherical tesseroid formulae can also be used in an ellipsoidal arrangement. An ellipsoidal tesseroid can be approximated by a spherical one by using r 1 = r E (ε 2 ) and r 2 = r E (ε 1 ), respectively, referring Eq. (33) to the mean latitude (ϕ 1 + ϕ 2 )/2.
In order to check the tesseroid approach in the ellipsoidal context, we have discretised an ellipsoidal shell between the confocal ellipsoids ε 1 and ε 2 , using a grid size ϕ = λ = 30 and a 1 −a 2 = 10 km, identical with the tesseroid dimensions. The differences between the results of the tesseroid approach and the rigorous formulae for the potential and its radial derivative are presented in Figs. 10 and 11, showing the approximation error along a meridian on the upper boundary surface ε 1 .  The absolute difference between the potential values calculated from Eq. (35) and the tesseroid method (Heck and Seitz 2007;Grombein et al. 2013) indicates that generally the tesseroid approach agrees better than 10 −5 m 2 s −2 .
The maximal numerical discrepancy occurs at the pole and reaches an absolute value of 0.027 m 2 s −2 . Keep in mind, that at the poles, the tesseroids degenerate into triangular columns.
The same conclusions can be drawn for the radial derivative of the potential (gravity attraction). For the chosen numerical example (a 1 − a 2 = d = 10 km, ϑ = λ = 30 , h P = 0 m), the numerical tesseroid solution agrees with Eq. (38) to an order of magnitude better than 10 −6 mGal, besides in the direct vicinity of the poles. Both the analytical solution for the ellipsoidal shell and the gravity forward modelling by tesseroids-as well as the respective software implementations-are thus confirmed. This mutual confirmation is of course dependent on the used geometrical fineness of discretization of the ellipsoidal shell. In order to touch this important topic, which is not in the focus of this paper, the approximation errors obtained on a coarser discretisation with ϕ × λ = {3 × 3 , 2 × 2 , 1 × 1 } are presented in Figs. 10 and 11. The effect of the different discretisations can be clearly seen here. While the coarser discretisation leads to systematic deviations with a very smooth progression, the result with a discretisation of 30 × 30 merely shows random behaviour. This is caused by rounding errors and reflects only numerical effects. In addition, this confirms that the used fine discretisation of 30 × 30 is fine enough for our numerical investigations. Fig. 11 Approximation error for the radial derivative of the potential of an ellipsoidal homogeneous shell, calculated on the surface of ε 1 along a meridian when using tesseroids. a 1 − a 2 = 10 km, ρ = 2670 kg m −3 , ϕ = λ = {3 (red), 2 (blue), 1 (green), 30 (black) }. The exact solution which the tesseroid evaluation is compared to, is calculated using Eq. (38) in spherical harmonics with a maximum degree of expansion N = 20

Conclusion
We have derived the external gravitational potential of a homogeneous ellipsoidal shell, bounded by two confocal ellipsoids of revolution, by solving the respective Newton integral in terms of spherical and elliptical coordinates, resulting in spherical and ellipsoidal harmonics. While the spherical harmonic representation produces a numerically fast-converging series for a global homogeneous ellipsoidal shell of, for example, 10 km thickness, which is used as an extreme case in the presented numerical investigations, the development in ellipsoidal harmonics consists of a finite series of degrees zero and two. As a matter of fact, both representations are theoretically equivalent, providing the same numerical results (up to rounding effects) if the spherical harmonic series is truncated at a suitable maximal degree. The series development in spherical harmonics enables the representation of the spectral components of the potential and the gradients, which may be considered as an advantage compared with the representation in ellipsoidal harmonics.
In Sects. 3.1 and 3.2, the derivation procedure starts from the solution of the Newton integral for computation points situated on the rotation axis, for which a closed analytical solution exists in spherical as well as in elliptical coordinates; in a second step, this solution is extended to the external space by harmonic continuation. The comparison between the rigorous finite ellipsoidal harmonic solution and the spherical harmonic expansion proves that the series expansion up to degree n = 6 is sufficient for guaranteeing an approximation error less than 10 −5 m 2 s −2 . This truncated series is easy to handle and contains only four terms. Thus, in practical applications the truncated spherical harmonic series will be preferred due to its higher numerical efficiency in comparison with the closed, but complicated analytical expressions in elliptical coordinates. The latter shows in the presented investigations some numerical problems, especially for the gravity gradient. The closed solution in elliptical coordinates consists only of two terms, but time-consuming arctan(x) function calls appear. Another disadvantage is that the spectral components cannot be used explicitly.
The derived formulae can be applied for testing and verification of computer programs and software implementations used for forward and inverse modelling in gravity field studies. In contrast to the often used isotropic test scenario based on a homogeneous spherical shell, the gravitational field of an ellipsoidal shell is also dependent on the horizontal position of the computation point, thus providing a more general control. As a matter of fact, with the presented ellipsoidal approach it is now possible to check also the gradients in latitudinal direction.
In Sect. 4.2, we have checked the well-established tesseroid approach (Heck and Seitz 2007;Grombein et al. 2013) and respective software by using a homogeneous ellipsoidal shell as test body representing the Earth's topography. Discretising the ellipsoidal shell by tesseroidal volume elements of different sizes (3 , 2 , 1 , 30 ) and summing up their contributions to the gravitational potential provide a numerical solution that has been compared with the rigorous analytical solution. In this way, the validity of the approach and the performance of the evaluation software could be proved. In the case of a 10 km-thick ellipsoidal shell and a discretisation of 30 × 30 , the agreement at a computation point on the ellipsoidal surface is generally better than 10 −5 m 2 s −2 ; the maximal discrepancy occurs at the poles, amounting to 0.027 m 2 s −2 , due to the degeneration of the tesseroids in the near zone to (quasi-)triangular elements. The results obtained with the coarser resolution clearly show the obvious dependency on the geometrical discretisation used. The latter is also an important outcome of the presented investigations. The derived solutions for potential and gravitational attraction for an ellipsoidal shell can be used to get a reliable statement on the results derived from DTMs of different resolutions.
Finally, also software for spherical harmonic analysis and synthesis can be controlled by the aid of the formulae derived in our paper. comments. The effort and cooperation of the Associate Editor is kindly acknowledged.
Funding Open Access funding enabled and organized by Projekt DEAL.
Data availability All data generated or analysed during this study are included in this published article.
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/.

A.3 Higher-order radial derivatives of the potential in spherical coordinates
Based on Eq. (27), the radial derivative of any degree k > 0 at a point in the space outside the ellipsoid can be expressed in the general formula: The first and second radial derivatives are explicitly given in Eqs. (31) and (32).

A.4 Convergence of the infinite spherical harmonic series
The fast convergence of the potential coefficients can be derived analytically as shown in the following. Starting with Eq. (30), the analytical expression of the potential coefficients in the basis of spherical harmonics is: 2n + 1 = 3(−1) n/2 (n + 1)(n + 3) √ 2n + 1 E R n , 0 ≤ n even .
The ratio between V n+2 and V n is given straightforward by V n+2 V n = 3(−1) n/2 3(−1) n/2 (−1) (n + 1)(n + 3) (n + 3)(n + 5) It can be seen that the magnitude of the n-dependent prefactor in Eq. (101) is always smaller than one. If we assume that ε 1 is the surface of the Geodetic Reference System 1980 (Moritz 1980b), we set R = a, without loss of generality, and note that E = ae, the final expression of the convergence of the potential coefficients is given by: The numerical behaviour of the coefficients expressing the potential of an ellipsoidal shell in spherical harmonics is clearly shown in Fig. 5. The rapid decay is obvious from the provided numbers in Table 2.