# Geopotential spherical harmonic model

In geophysics and physical geodesy, a geopotential model is the theoretical analysis of measuring and calculating the effects of Earth's gravitational field (the geopotential). The Earth is not exactly spherical, mainly because of its rotation around the polar axis that makes its shape slightly oblate. However, a spherical harmonics series expansion captures the actual field with increasing fidelity.
If Earth's shape were perfectly known together with the exact mass density ρ = ρ(x, y, z), it could be integrated numerically (when combined with a reciprocal distance kernel) to find an accurate model for Earth's gravitational field. However, the situation is in fact the opposite: by observing the orbits of spacecraft and the Moon, Earth's gravitational field can be determined quite accurately. The best estimate of Earth's mass is obtained by dividing the product GM as determined from the analysis of spacecraft orbit with a value for the gravitational constant G, determined to a lower relative accuracy using other physical methods.

## Background

From the defining equations (1) and (2) it is clear (taking the partial derivatives of the integrand) that outside the body in empty space the following differential equations are valid for the field caused by the body:
| ${\frac {\partial F_{x}}{\partial x}}+{\frac {\partial F_{y}}{\partial y}}+{\frac {\partial F_{z}}{\partial z}}=0$ | (5) |
| ${\frac {\partial ^{2}u}{\partial x^{2}}}+{\frac {\partial ^{2}u}{\partial y^{2}}}+{\frac {\partial ^{2}u}{\partial z^{2}}}=0$ | (6) |
Functions of the form $\phi =R(r)\,\Theta (\theta )\,\Phi (\varphi )$ where (r, θ, φ) are the spherical coordinates which satisfy the partial differential equation (6) (the Laplace equation) are called spherical harmonic functions.
They take the forms:
| ${\begin{aligned}g(x,y,z)&={\frac {1}{r^{n+1}}}P_{n}^{m}(\sin \theta )\cos m\varphi \,,&0&\leq m\leq n\,,&n&=0,1,2,\dots \\h(x,y,z)&={\frac {1}{r^{n+1}}}P_{n}^{m}(\sin \theta )\sin m\varphi \,,&1&\leq m\leq n\,,&n&=1,2,\dots \end{aligned}}$ | (7) |
where spherical coordinates (r, θ, φ) are used, given here in terms of cartesian (x, y, z) for reference:
| ${\begin{aligned}x&=r\cos \theta \cos \varphi \\y&=r\cos \theta \sin \varphi \\z&=r\sin \theta \,,\end{aligned}}$ | (8) |
also P$^0$$_n$ are the Legendre polynomials and P$^m$$_n$ for 1 ≤ m ≤ n are the associated Legendre functions.
The first spherical harmonics with n = 0, 1, 2, 3 are presented in the table below. [Note that the sign convention differs from the one in the page about the associated Legendre polynomials, here $P_{2}^{1}(x)=3x{\sqrt {1-x^{2}}}$ whereas there $P_{2}^{1}(x)=-3x{\sqrt {1-x^{2}}}$.]
    nSpherical harmonics 0 ${\frac {1}{r}}$ 1 ${\frac {1}{r^{2}}}P_{1}^{0}(\sin \theta )={\frac {1}{r^{2}}}\sin \theta$ ${\frac {1}{r^{2}}}P_{1}^{1}(\sin \theta )\cos \varphi ={\frac {1}{r^{2}}}\cos \theta \cos \varphi$ ${\frac {1}{r^{2}}}P_{1}^{1}(\sin \theta )\sin \varphi ={\frac {1}{r^{2}}}\cos \theta \sin \varphi$ 2 ${\frac {1}{r^{3}}}P_{2}^{0}(\sin \theta )={\frac {1}{r^{3}}}{\frac {1}{2}}(3\sin ^{2}\theta -1)$ ${\frac {1}{r^{3}}}P_{2}^{1}(\sin \theta )\cos \varphi ={\frac {1}{r^{3}}}3\sin \theta \cos \theta \ \cos \varphi$ ${\frac {1}{r^{3}}}P_{2}^{1}(\sin \theta )\sin \varphi ={\frac {1}{r^{3}}}3\sin \theta \cos \theta \sin \varphi$ ${\frac {1}{r^{3}}}P_{2}^{2}(\sin \theta )\cos 2\varphi ={\frac {1}{r^{3}}}3\cos ^{2}\theta \ \cos 2\varphi$ ${\frac {1}{r^{3}}}P_{2}^{2}(\sin \theta )\sin 2\varphi ={\frac {1}{r^{3}}}3\cos ^{2}\theta \sin 2\varphi$ 3 ${\frac {1}{r^{4}}}P_{3}^{0}(\sin \theta )={\frac {1}{r^{4}}}{\frac {1}{2}}\sin \theta \ (5\sin ^{2}\theta -3)$ ${\frac {1}{r^{4}}}P_{3}^{1}(\sin \theta )\cos \varphi ={\frac {1}{r^{4}}}{\frac {3}{2}}\ (5\ \sin ^{2}\theta -1)\cos \theta \cos \varphi$ ${\frac {1}{r^{4}}}P_{3}^{1}(\sin \theta )\sin \varphi ={\frac {1}{r^{4}}}{\frac {3}{2}}\ (5\ \sin ^{2}\theta -1)\cos \theta \sin \varphi$ ${\frac {1}{r^{4}}}P_{3}^{2}(\sin \theta )\cos 2\varphi ={\frac {1}{r^{4}}}15\sin \theta \cos ^{2}\theta \cos 2\varphi$ ${\frac {1}{r^{4}}}P_{3}^{2}(\sin \theta )\sin 2\varphi ={\frac {1}{r^{4}}}15\sin \theta \cos ^{2}\theta \sin 2\varphi$ ${\frac {1}{r^{4}}}P_{3}^{3}(\sin \theta )\cos 3\varphi ={\frac {1}{r^{4}}}15\cos ^{3}\theta \cos 3\varphi$ ${\frac {1}{r^{4}}}P_{3}^{3}(\sin \theta )\sin 3\varphi ={\frac {1}{r^{4}}}15\cos ^{3}\theta \sin 3\varphi$

## Formulation

The model for Earth's gravitational potential is a sum
| $u=-{\frac {\mu }{r}}+\sum _{n=2}^{N_{z}}{\frac {J_{n}P_{n}^{0}(\sin \theta )}{r^{n+1}}}+\sum _{n=2}^{N_{t}}\sum _{m=1}^{n}{\frac {P_{n}^{m}(\sin \theta )\left(C_{n}^{m}\cos m\varphi +S_{n}^{m}\sin m\varphi \right)}{r^{n+1}}}$ | (9) |
where $\mu =GM$ and the coordinates (8) are relative to the standard geodetic reference system extended into space with origin in the center of the reference ellipsoid and with z-axis in the direction of the polar axis.
The zonal terms refer to terms of the form:
    ${\frac {P_{n}^{0}(\sin \theta )}{r^{n+1}}}\quad n=0,1,2,\dots$
and the tesseral terms terms refer to terms of the form:
    ${\frac {P_{n}^{m}(\sin \theta )\cos m\varphi }{r^{n+1}}}\,,\quad 1\leq m\leq n\quad n=1,2,\dots$
    ${\frac {P_{n}^{m}(\sin \theta )\sin m\varphi }{r^{n+1}}}$
The zonal and tesseral terms for n = 1 are left out in (9). The coefficients for the n=1 with both m=0 and m=1 term correspond to an arbitrarily oriented dipole term in the multi-pole expansion. Gravity does not physically exhibit any dipole character and so the integral characterizing n = 1 must be zero.
The different coefficients J$_n$, C$_n$$^m$, S$_n$$^m$, are then given the values for which the best possible agreement between the computed and the observed spacecraft orbits is obtained.
As P$^0$$_n$(x) = −P$^0$$_n$(−x) non-zero coefficients J$_n$ for odd n correspond to a lack of symmetry "north–south" relative the equatorial plane for the mass distribution of Earth. Non-zero coefficients C$_n$$^m$, S$_n$$^m$ correspond to a lack of rotational symmetry around the polar axis for the mass distribution of Earth, i.e. to a "tri-axiality" of Earth.
For large values of n the coefficients above (that are divided by r$^{(n+1)}$ in (9)) take very large values when for example kilometers and seconds are used as units. In the literature it is common to introduce some arbitrary "reference radius" R close to Earth's radius and to work with the dimensionless coefficients
    ${\begin{aligned}{\tilde {J_{n}}}&=-{\frac {J_{n}}{\mu \ R^{n}}},&{\tilde {C_{n}^{m}}}&=-{\frac {C_{n}^{m}}{\mu \ R^{n}}},&{\tilde {S_{n}^{m}}}&=-{\frac {S_{n}^{m}}{\mu \ R^{n}}}\end{aligned}}$
and to write the potential as
| $u=-{\frac {\mu }{r}}\left(1+\sum _{n=2}^{N_{z}}{\frac {{\tilde {J_{n}}}P_{n}^{0}(\sin \theta )}{{({\frac {r}{R}})}^{n}}}+\sum _{n=2}^{N_{t}}\sum _{m=1}^{n}{\frac {P_{n}^{m}(\sin \theta )({\tilde {C_{n}^{m}}}\cos m\varphi +{\tilde {S_{n}^{m}}}\sin m\varphi )}{{({\frac {r}{R}})}^{n}}}\right)$ | (10) |

### Derivation

A compact derivation of the spherical harmonics used to model Earth's gravitational field. :
| $\phi =R(r)\ \Theta (\theta )\ \Phi (\varphi )$ | (16) |
| ${\frac {\partial ^{2}f}{\partial x^{2}}}\ +\ {\frac {\partial ^{2}f}{\partial y^{2}}}\ +\ {\frac {\partial ^{2}f}{\partial z^{2}}}\ =\ {1 \over r^{2}}{\partial \over \partial r}\left(r^{2}{\partial f \over \partial r}\right)+{1 \over r^{2}\cos \theta }{\partial \over \partial \theta }\left(\cos \theta {\partial f \over \partial \theta }\right)+{1 \over r^{2}\cos ^{2}\theta }{\partial ^{2}f \over \partial \varphi ^{2}}$ | (17) |
| ${\frac {r^{2}}{\phi }}\left({\frac {\partial ^{2}\phi }{\partial x^{2}}}+{\frac {\partial ^{2}\phi }{\partial y^{2}}}+{\frac {\partial ^{2}\phi }{\partial z^{2}}}\right)\ ={\frac {1}{R}}{\frac {d}{dr}}\left(r^{2}{\frac {dR}{dr}}\right)+{\frac {1}{\Theta \cos \theta }}{\frac {d}{d\theta }}\left(\cos \theta {\frac {d\Theta }{d\theta }}\right)+{\frac {1}{\Phi \cos ^{2}\theta }}{\frac {d^{2}\Phi }{d\varphi ^{2}}}$ | (18) |
| ${\frac {1}{R}}{\frac {d}{dr}}\left(r^{2}{\frac {dR}{dr}}\right)=\lambda$ | (19) |
| ${\frac {1}{\Theta \cos \theta }}{\frac {d}{d\theta }}\left(\cos \theta {\frac {d\Theta }{d\theta }}\right)+{\frac {1}{\Phi \cos ^{2}\theta }}{\frac {d^{2}\Phi }{d\varphi ^{2}}}=-\lambda$ | (20) |
| ${\frac {1}{\Phi }}{\frac {d^{2}\Phi }{d\varphi ^{2}}}=-m^{2}$ | (21) |
| ${\frac {1}{\Theta }}\ \cos \theta \ {\frac {d}{d\theta }}\left(\cos \theta {\frac {d\Theta }{d\theta }}\right)\ +\lambda \ \cos ^{2}\theta =m^{2}$ | (22) |
| $\Phi (\varphi )\ =a\ \cos m\varphi \ +\ b\ \sin m\varphi$ | (23) |
| ${\frac {d}{dx}}\left(\left(1-x^{2}\right){\frac {d\Theta }{dx}}\right)+\left(\lambda -{\frac {m^{2}}{1-x^{2}}}\right)\Theta =0$ | (24) |
| ${\frac {d}{dx}}\left(\left(1-x^{2}\right)\ {\frac {dP_{n}}{dx}}\right)\ +\ n(n+1)\ P_{n}=0$ | (25) |
| ${\frac {d}{dx}}\left(\left(1-x^{2}\right)\ {\frac {dP_{n}^{m}}{dx}}\right)\ +\ \left(n(n+1)-{\frac {m^{2}}{1-x^{2}}}\right)\ P_{n}^{m}\ =\ 0$ | (26) |
| $\phi ={\frac {1}{r^{n+1}}}\ P_{n}^{m}(\sin \theta )\ (a\ \cos m\varphi \ +\ b\ \sin m\varphi )$ | (27) |
| ${\begin{aligned}P_{0}(x)&=1\\P_{n}(x)&={\frac {1}{2^{n}n!}}\ {\frac {d^{n}\left(x^{2}-1\right)^{n}}{dx^{n}}}\quad n\geq 1\\\end{aligned}}$ | (28) |
| ${\begin{aligned}P_{0}(x)&=1\\P_{1}(x)&=x\\P_{2}(x)&={\frac {1}{2}}\left(3x^{2}-1\right)\\P_{3}(x)&={\frac {1}{2}}\left(5x^{3}-3x\right)\\P_{4}(x)&={\frac {1}{8}}\left(35x^{4}-30x^{2}+3\right)\\P_{5}(x)&={\frac {1}{8}}\left(63x^{5}-70x^{3}+15x\right)\\\end{aligned}}$ | (29) |
| $P_{n}^{m}(x)=\left(1-x^{2}\right)^{\frac {m}{2}}\ {\frac {d^{m}P_{n}}{dx^{m}}}\quad 1\leq m\leq n$ | (30) |

## Largest terms

The dominating term (after the term −μ/r) in (9) is the J$_2$ coefficient, the second dynamic form factor representing the oblateness of Earth:
    $u={\frac {J_{2}\ P_{2}^{0}(\sin \theta )}{r^{3}}}=J_{2}{\frac {1}{r^{3}}}{\frac {1}{2}}(3\sin ^{2}\theta -1)=J_{2}{\frac {1}{r^{5}}}{\frac {1}{2}}(3z^{2}-r^{2})$
Relative the coordinate system
| ${\begin{aligned}{\hat {\varphi }}&=-\sin \varphi {\hat {x}}+\cos \varphi {\hat {y}}\\{\hat {\theta }}&=-\sin \theta \,(\cos \varphi {\hat {x}}+\sin \varphi {\hat {y}})+\cos \theta {\hat {z}}\\{\hat {r}}&=\cos \theta \,(\cos \varphi {\hat {x}}+\sin \varphi {\hat {y}})+\sin \theta {\hat {z}}\end{aligned}}$ | (11) |
illustrated in figure 1 the components of the force caused by the "J$_2$ term" are
| ${\begin{aligned}F_{\theta }&=-{\frac {1}{r}}\ {\frac {\partial u}{\partial \theta }}=-J_{2}\ {\frac {1}{r^{4}}}3\ \cos \theta \ \sin \theta \\F_{r}&=-{\frac {\partial u}{\partial r}}=J_{2}\ {\frac {1}{r^{4}}}{\frac {3}{2}}\ \left(3\sin ^{2}\theta \ -\ 1\right)\end{aligned}}$ | (12) |
In the rectangular coordinate system (x, y, z) with unit vectors (x̂ ŷ ẑ) the force components are:
| ${\begin{aligned}F_{x}&=-{\frac {\partial u}{\partial x}}=J_{2}\ {\frac {x}{r^{7}}}\left(6z^{2}-{\frac {3}{2}}(x^{2}+y^{2})\right)\\F_{y}&=-{\frac {\partial u}{\partial y}}=J_{2}\ {\frac {y}{r^{7}}}\left(6z^{2}-{\frac {3}{2}}(x^{2}+y^{2})\right)\\F_{z}&=-{\frac {\partial u}{\partial z}}=J_{2}\ {\frac {z}{r^{7}}}\left(3z^{2}-{\frac {9}{2}}(x^{2}+y^{2})\right)\end{aligned}}$ | (13) |
The components of the force corresponding to the "J$_3$ term"
    $u={\frac {J_{3}P_{3}^{0}(\sin \theta )}{r^{4}}}=J_{3}{\frac {1}{r^{4}}}{\frac {1}{2}}\sin \theta \left(5\sin ^{2}\theta -3\right)=J_{3}{\frac {1}{r^{7}}}{\frac {1}{2}}z\left(5z^{2}-3r^{2}\right)$
are
| ${\begin{aligned}F_{\theta }&=-{\frac {1}{r}}{\frac {\partial u}{\partial \theta }}=-J_{3}{\frac {1}{r^{5}}}{\frac {3}{2}}\cos \theta \left(5\sin ^{2}\theta -1\right)\\F_{r}&=-{\frac {\partial u}{\partial r}}=J_{3}{\frac {1}{r^{5}}}2\sin \theta \left(5\sin ^{2}\theta -3\right)\end{aligned}}$ | (14) |
and
| ${\begin{aligned}F_{x}&=-{\frac {\partial u}{\partial x}}=J_{3}{\frac {xz}{r^{9}}}\left(10z^{2}-{\frac {15}{2}}(x^{2}+y^{2})\right)\\F_{y}&=-{\frac {\partial u}{\partial y}}=J_{3}{\frac {yz}{r^{9}}}\left(10z^{2}-{\frac {15}{2}}(x^{2}+y^{2})\right)\\F_{z}&=-{\frac {\partial u}{\partial z}}=J_{3}{\frac {1}{r^{9}}}\left(4z^{2}\ \left(z^{2}-3(x^{2}+y^{2})\right)+{\frac {3}{2}}(x^{2}+y^{2})^{2}\right)\end{aligned}}$ | (15) |
The exact numerical values for the coefficients deviate (somewhat) between different Earth models but for the lowest coefficients they all agree almost exactly.
For the JGM-3 model (see below) the values are:
    μ = 398600.440 km$^3$⋅s$^{−2}$
    J$_2$ = 1.75553 × 10$^{10}$ km$^5$⋅s$^{−2}$
    J$_3$ = −2.61913 × 10$^{11}$ km$^6$⋅s$^{−2}$
For example, at a radius of 6600 km (about 200 km above Earth's surface) J$_3$/(J$_2$r) is about 0.002; i.e., the correction to the "J$_2$ force" from the "J$_3$ term" is in the order of 2 permille. The negative value of J$_3$ implies that for a point mass in Earth's equatorial plane the gravitational force is tilted slightly towards the south due to the lack of symmetry for the mass distribution of Earth's "north–south".

## Recursive algorithms used for the numerical propagation of spacecraft orbits

Spacecraft orbits are computed by the numerical integration of the equation of motion. For this the gravitational force, i.e. the gradient of the potential, must be computed. Efficient recursive algorithms have been designed to compute the gravitational force for any $N_{z}$ and $N_{t}$ (the max degree of zonal and tesseral terms) and such algorithms are used in standard orbit propagation software.

## Available models

The earliest Earth models in general use by NASA and ESRO/ESA were the "Goddard Earth Models" developed by Goddard Space Flight Center (GSFC) denoted "GEM-1", "GEM-2", "GEM-3", and so on. Later the "Joint Earth Gravity Models" denoted "JGM-1", "JGM-2", "JGM-3" developed by GSFC in cooperation with universities and private companies became available. The newer models generally provided higher order terms than their precursors. The EGM96 uses N$_z$ = N$_t$ = 360 resulting in 130317 coefficients. An EGM2008 model is available as well.
For a normal Earth satellite requiring an orbit determination/prediction accuracy of a few meters the "JGM-3" truncated to N$_z$ = N$_t$ = 36 (1365 coefficients) is usually sufficient. Inaccuracies from the modeling of the air-drag and to a lesser extent the solar radiation pressure will exceed the inaccuracies caused by the gravitation modeling errors.
The dimensionless coefficients ${\tilde {J_{n}}}=-{\frac {J_{n}}{\mu \ R^{n}}}$, ${\tilde {C_{n}^{m}}}=-{\frac {C_{n}^{m}}{\mu \ R^{n}}}$, ${\tilde {S_{n}^{m}}}=-{\frac {S_{n}^{m}}{\mu \ R^{n}}}$ for the first zonal and tesseral terms (using $R$ = 6378.1363 km and $\mu$ = 398600.4415 km$^3$/s$^2$) of the JGM-3 model are
Zonal coefficients :
| n | |
| 2 | −0.1082635854×10$^{−2}$ |
| 3 | 0.2532435346×10$^{−5}$ |
| 4 | 0.1619331205×10$^{−5}$ |
| 5 | 0.2277161016×10$^{−6}$ |
| 6 | −0.5396484906×10$^{−6}$ |
| 7 | 0.3513684422×10$^{−6}$ |
| 8 | 0.2025187152×10$^{−6}$ |
Tesseral coefficients :
| n | m | C | S |
| 2 | 1 | −0.3504890360×10$^{−9}$ | 0.1635406077×10$^{−8}$ |
| 2 | 2 | 0.1574536043×10$^{−5}$ | −0.9038680729×10$^{−6}$ |
| 3 | 1 | 0.2192798802×10$^{−5}$ | 0.2680118938×10$^{−6}$ |
| 3 | 2 | 0.3090160446×10$^{−6}$ | −0.2114023978×10$^{−6}$ |
| 3 | 3 | 0.1005588574×10$^{−6}$ | 0.1972013239×10$^{−6}$ |
| 4 | 1 | −0.5087253036×10$^{−6}$ | −0.4494599352×10$^{−6}$ |
| 4 | 2 | 0.7841223074×10$^{−7}$ | 0.1481554569×10$^{−6}$ |
| 4 | 3 | 0.5921574319×10$^{−7}$ | −0.1201129183×10$^{−7}$ |
| 4 | 4 | −0.3982395740×10$^{−8}$ | 0.6525605810×10$^{−8}$ |
According to JGM-3 one therefore has that J$_2$ = 0.1082635854×10$^{−2}$ × 6378.1363$^2$ × 398600.4415 km$^5$/s$^2$ = 1.75553×10$^{10}$ km$^5$/s$^2$ and J$_3$ = −0.2532435346×10$^{−5}$ × 6378.1363$^3$ × 398600.4415 km$^6$/s$^2$ = −2.61913×10$^{11}$ km$^6$/s$^2$.