Frontiers reaches 6.4 on Journal Impact Factors

Original Research ARTICLE

Front. Astron. Space Sci., 25 February 2016 |

Reference Ellipsoid and Geoid in Chronometric Geodesy

  • 1Department of Physics and Astronomy, University of Missouri, Columbia, MO, USA
  • 2Department of Physical Geodesy and Remote Sensing, Siberian State University of Geosystems and Technologies, Novosibirsk, Russia

Chronometric geodesy applies general relativity to study the problem of the shape of celestial bodies including the earth, and their gravitational field. The present paper discusses the relativistic problem of construction of a background geometric manifold that is used for describing a reference ellipsoid, geoid, the normal gravity field of the earth and for calculating geoid's undulation (height). We choose the perfect fluid with an ellipsoidal mass distribution uniformly rotating around a fixed axis as a source of matter generating the geometry of the background manifold through the Einstein equations. We formulate the post-Newtonian hydrodynamic equations of the rotating fluid to find out the set of algebraic equations defining the equipotential surface of the gravity field. In order to solve these equations we explicitly perform all integrals characterizing the interior gravitational potentials in terms of elementary functions depending on the parameters defining the shape of the body and the mass distribution. We employ the coordinate freedom of the equations to choose these parameters to make the shape of the rotating fluid configuration to be an ellipsoid of rotation. We derive expressions of the post-Newtonian mass and angular momentum of the rotating fluid as functions of the rotational velocity and the parameters of the ellipsoid including its bare density, eccentricity and semi-major axes. We formulate the post-Newtonian Pizzetti and Clairaut theorems that are used in geodesy to connect the parameters of the reference ellipsoid to the polar and equatorial values of force of gravity. We expand the post-Newtonian geodetic equations characterizing the reference ellipsoid into the Taylor series with respect to the eccentricity of the ellipsoid, and discuss the small-eccentricity approximation. Finally, we introduce the concept of relativistic geoid and its undulation with respect to the reference ellipsoid, and discuss how to calculate it in chronometric geodesy by making use of the anomalous gravity potential.

1. Introduction

Accurate definition, determination and realization of terrestrial reference frame is essential for fundamental astronomy, celestial mechanics, geophysics as well as for precise satellite and aircraft navigation, positioning, and mapping. The International Terrestrial Reference Frame (ITRF) is materialized by coordinates of a number of geodetic points (stations) located on the Earth's surface, and spread out across the globe. Definition of the coordinates is derived from the International Terrestrial Reference System (ITRS) which is a theoretical concept. ITRF is used to measure plate tectonics, regional subsidence or loading and to describe the Earth's rotation in space by measuring its rotational parameters (Petit and Luzum, 2010).

Nowadays, four main geodetic techniques are used to compute the accurate terrestrial coordinates and velocities of stations—GPS, VLBI, SLR, and DORIS, for the realization of ITRF. The observations are so accurate that geodesists have to model and include to the data processing the secular Earth's crust changes to reach self-consistency between the various ITRF realizations referred to different epoch. It is recognized that in order to maintain the up-to-date ITRF realization as accurate as possible the development of the most precise theoretical model and relationships between parameters of the model is of a paramount importance.

This model is not of a pure kinematic origin because the gravity field plays an essential role in geodetic network computations (Heiskanen and Moritz, 1967; Hofmann-Wellenhof and Moritz, 2006) making a particular theory of gravity a part of definition of the ITRS. Until recently, such a theory was the Newtonian theory of gravity. However, current SLR/GPS/VLBI techniques allow us to determine the transformation parameters between the coordinates and velocities of the collocation points of the ITRF realization with the precision approaching to ~1 mm and ~1 mm/year, respectively (Petit and Luzum, 2010, Table 4.1). At the same time, general relativity predicts that the post-Newtonian effects in measuring the ITRF coordinates, as compared with the Newtonian theory of gravity, is of the order of the Earth's gravitational radius that is about 1 cm (Kopejkin, 1991; Müller et al., 2008). This general-relativistic effect distorts the coordinate grid of ITRF at the noticeable level which can be detected with the currently available geodetic techniques and, hence, should be taken into account.

ITRS is specified by the Cartesian equatorial coordinates rotating rigidly with the Earth. For the purposes of geodesy and gravimetry the Cartesian coordinates are often converted to a geographical system of elliptical coordinates (h, θ, λ) referred to an international reference ellipsoid which is a non-trivial solution of the Newtonian gravity field equations found by Maclaurin (Chandrasekhar, 1967). It describes the figure of a fluid body with a homogeneous mass density rotating around a fixed z-axis with a constant angular velocity ω. The post-Newtonian effects deforms the shape of Maclaurin's ellipsoid (Chandrasekhar, 1965; Bardeen, 1971) and modify the basic equations of classic geodesy (Müller et al., 2008; Kopeikin et al., 2011, 2016). These relativistic effects must be properly calculated to ensure the adequacy of the geodetic coordinate transformations at the millimeter level of accuracy. A pioneering study of relativistic effects in geodesy have been carried out by Bjerhammar (1985).

Post-Newtonian equilibrium configurations of uniformly rotating fluids have been discussed in literature by researchers from USA (Chandrasekhar, 1965, 1967a,b,c, 1971a,b; Bardeen, 1971; Chandrasekhar and Elbert, 1974, 1978; Chandrasekhar and Miller, 1974), Lithuania (Bondarenko and Pyragas, 1974; Pyragas et al., 1974, 1975), USSR (Tsirulev and Tsvetkov, 1982a,b; Tsvetkov and Tsirulev, 1983; Galtsov et al., 1984) and, the most recently, scientists from the Fridrich Schiller University of Jena in Germany (Petroff, 2003; Ansorg et al., 2004; Meinel et al., 2008; Gürlebeck and Petroff, 2010, 2013). These papers focused primarily on studying the astrophysical aspects of the problem like stability of the rotating stars, the points of bifurcations, exact axially-symmetric spacetimes, emission of gravitational waves, etc. The present paper focuses on the post-Newtonian effects in physical geodesy. More specifically, we extend the research on the figures of equilibrium into the realm of relativistic geodesy and pay attention to the possible geodetic applications for the adequate processing of the high-precise data obtained by various geodetic techniques that include but not limited to SLR, LLR, VLBI, DORIS, and GNSS (Drewes, 2009; Plag and Pearlman, 2009).

Important stimulating factor for pursuing more advanced theoretical research on relativistic effects in geodesy and the Earth figure of equilibrium is related to the recent breakthrough in manufacturing ring laser gyroscopes (Schreiber, 2013; Beverini et al., 2014; Hurst et al., 2015). Earth rotation and orientation provide the link between the terrestrial (ITRF) and celestial reference frames (ICRF). Traditionally, the Earth orientation parameters (EOPs) are observed by the International Earth Rotation Service (IERS) that operates the global VLBI network. VLBI technique determines EOPs off-line after accumulating the data over sufficiently long period of time. Moreover, VLBI depends on suitable models of Earth's troposphere and geophysical factors which are difficult to predict. Ring laser gyroscopes provide direct access to the Earth rotation axis and a high resolution in the short-term. Currently, the modern quantum optics has matured to a point where the ring laser gyroscopes can resolve rotation rates of 10−12 rad/s after 1 h of integration and demonstrate an impressive stability over several month (Schreiber, 2013). The combination of VLBI and ring laser measurements offers an improved sensitivity for the EOPs in the short-term and the direct access to the Earth rotation axis. Another application of the ring laser gyroscopes is to measure the Lense-Thirring effect predicted by Einstein's general relativity (Ciufolini and Wheeler, 1995), in a terrestrial laboratory environment (Beverini et al., 2014).

Another, rapidly emerging branch of relativistic geodesy is called chronometric geodesy (Petit et al., 2014). This new line of development is pushed forward by the fascinating progress in making up quantum clocks (Poli et al., 2013), ultra-precise time-scale dissemination over the globe (Kómár et al., 2014), and geophysical applications of the clocks (Bondarescu et al., 2012, 2015). Clocks at rest in a gravitational potential tick slower than clocks outside of it. On Earth, this translates to a relative frequency change of 10−16 per meter of height difference (Falke et al., 2014). Comparing the frequency of a probe clock with a reference clock provides a direct measure of the gravity potential difference between the two clocks. This novel technique has been dubbed chronometric leveling. It is envisioned as one of the most promising application of the relativistic geodesy in a near future (Mai and Müller, 2013; Mai, 2014; Petit et al., 2014). Optical frequency standards have recently reached stability of 2.2 × 10−16 at 1 s, and demonstrated an overall fractional frequency uncertainty of 2.1 × 10−18 (Nicholson et al., 2015) which enables their use for chronometric geodesy at the absolute level of one centimeter.

The chronometric leveling directly measures the equipotential surface of gravity field (geoid) without conducting a complicated gravimetric survey and solving the differential equations for the anomalous gravity potentials (Stoke's or Molodensky's problem Hofmann-Wellenhof and Moritz, 2006; Torge and Müller, 2012). Combining the data of the chronometric leveling with those of the conventional geodetic techniques will allow us to determine the normal heights of reference points with an unprecedented accuracy (Bondarescu et al., 2012). An adequate physical interpretation of this type of measurements requires more precise theoretical definition of geoid (Kopeikin et al., 2015; Oltean et al., 2015) and is inconceivable without an accompanying development of the corresponding mathematical algorithms accounting for the major relativistic effects in the description of the equipotential level surface and its, so-called, normal) gravity field (Hofmann-Wellenhof and Moritz, 2006; Torge and Müller, 2012).

The objective of the present paper is to suggest a solution of Einstein's gravity field equations that gives rise to the post-Newtonian equipotential surface described by the second-order polynomial which exactly coincides with the shape of the Maclaurin ellipsoid in the Newtonian gravity. Our previous paper (Kopeikin et al., 2016) has explored the post-Newtonian case of a homogeneous density distribution of a uniformly rotating perfect fluid. We have shown that the post-Newtonian effects inevitably distort the Maclaurin ellipsoid of rotation to a surface described by the fourth-order polynomial which cannot be reduced to the ellipsoidal surface of the second order by doing coordinate transformations. In the present paper we assume that the mass distribution is not homogeneous but an ellipsoidal one with a free parameter, and prove that under this assumption the equipotential level surface can be made as an exact ellipsoid of the second order which may be found more convenient for doing the reduction of high-precise geodetic measurements with taking into account relativistic corrections.

The paper consists of several sections, bibliography and appendix. Section 2 explains briefly the principles of the post-Newtonian approximations and describes the post-Newtonian metric tensor. In Section 3 we give definition of the post-Newtonian geoid. Section 4 discusses the post-Newtonian ellipsoid which generalizes the Maclaurin ellipsoid and is, in general, the surface of the fourth order. It also introduces the reader to the concept of the post-Newtonian coordinate freedom and shows how this freedom can be used to simplify the mathematical description of the PN-ellipsoid. Sections 5–7 calculate respectively the Newtonian and post-Newtonian gravitational potentials inside the rotating PN-ellipsoid. Sections 8, 9 give the post-Newtonian definitions of the conserved mass and angular momentum. Section 10 derives the post-Newtonian equations of the equipotential level surface. Sections 11, 12 provide the reader with the relativistic generalization of the Pizzetti and Clairaut theorem of classical geodesy (Pizzetti, 1913). Section 13 approximates the relativistic formulas which can be used in practical applications of relativistic geodesy and shows how to build the exact reference-ellipsoid in the post-Newtonian approximation. Finally, Section 14 introduces the concept of the post-Newtonian geoid's undulation and its relation to the current problem of the measurement of the global average sea level. Appendix contains technical details of calculations of the mathematical integrals which appear in the main part of the paper.

The following notations are used throughout the paper:

• the Greek indices α, β, … run from 0 to 3,

• the Roman indices i, j, … run from 1 to 3,

• repeated Greek indices mean Einstein's summation from 0 to 3,

• repeated Roman indices mean Einstein's summation from 1 to 3,

• the unit matrix (also known as the Kroneker symbol) is denoted by δij=δij,

• the fully antisymmetric symbol Levi-Civita is denoted as εijk=εijk with ε123 = + 1,

• the bold letters a = (a1, a2, a3) ≡ (ai), b = (b1, b2, b3) ≡ (bi), and so on, denote spatial 3-dimensional vectors,

• a dot between two spatial vectors, for example a·b=a1b1+a2b2+a3b3=δijaibj, means the Euclidean dot product,

• the cross between two vectors, for example (a × b)i ≡ εijkajbk, means the Euclidean cross product,

• we use a shorthand notation for partial derivatives α=xα,

• covariant derivative with respect to a coordinate xα is denoted as ∇α;

• the Minkowski (flat) space-time metric ηαβ = diag(−1, +1, +1, +1),

gαβ is the physical spacetime metric,

• ḡαβ is the background spacetime metric describing the gravitational field of the reference ellipsoid,

• the Greek indices are raised and lowered with the metric ḡαβ,

G is the universal gravitational constant,

c is the speed of light in vacuum,

• ω is a constant rotational velocity,

• ρ is a bare constant density of the fluid,

a is a semi-major axis of the ellipsoid of revolution,

b is a semi-minor axis of the ellipsoid of revolution,

• ρ0 is a constant bare density of matter,

κπGρ0a2c2 is a dimensional parameter characterizing the strength of gravitational field,

• a bar above any quantity indicates that it belongs to the background spacetime manifold with the metric ḡαβ.

Other notations are explained in the text as they appear.

2. Post-Newtonian Spacetime and Metric Tensor

Discussion of the chronometric geodesy starts from the construction of the spacetime manifold for the case of a rotating Earth. We shall employ Einstein's general relativity to build such a manifold though some other alternative theories of gravity discussed, for example in Will (1993), can be used as well. Einstein's field equations represent a system of ten non-linear differential equations in partial derivatives for the ten components of the metric tensor, gαβ, which represent gravitational potentials. Because the equations are difficult to solve exactly due to their non-linearity, we apply the post-Newtonian approximations (PNA) for their solution (Chandrasekhar, 1965).

The PNA are applied in case of slowly-moving matter having a weak gravitational field. This is exactly the situation in the solar system which makes PNA highly appropriate for constructing relativistic theory of reference frames (Soffel et al., 2003) and relativistic celestial mechanics in the solar system (Soffel, 1989; Brumberg, 1991; Kopeikin et al., 2011). The PNA are based on the assumption that a Taylor expansion of the metric tensor can be done in the inverse powers of the speed of gravity c that is equal to the speed of light in vacuum in general relativity. Exact mathematical formulation of a set of basic axioms required for doing the post-Newtonian expansion was given by Rendall (1990). Practically, it requires to have several small parameters characterizing the source of gravity. They are: εi ~ vic, εe ~ vec, and ηi~Uic2, ηe~Uec2, where vi is a characteristic velocity of motion of matter inside a body, ve is a characteristic velocity of the relative motion of the bodies with respect to each other, Ui is the internal gravitational potential of each body, and Ue is the external gravitational potential between the bodies. If one denotes a characteristic radius of a body as L and a characteristic distance between the bodies as R, the internal and external gravitational potentials will be UiGML and UeGMR, where M is a characteristic mass of the body. Due to the virial theorem of the Newtonian gravity (Chandrasekhar, 1965) the small parameters are not independent. Specifically, one has εi2~ηi and εe2~ηe. Hence, parameters εi and εe are sufficient in doing post-Newtonian approximations. Because within the solar system these parameters do not significantly differ from each other, we shall not distinguished them. Quite often we shall use notation κ ≡ πGρa2c2 to mark the powers of the fundamental speed c in the post-Newtonian terms.

We assume that physical spacetime has the metric tensor denoted as gαβ. This spacetime is well-approximated by a background manifold with the metric tensor denoted ḡαβ. The perturbation of gravitational field is denoted with ϰαβ and is called the anomalous gravity potential. We, first, build the metric tensor of the background manifold by solving Einstein's equations in harmonic coordinates xα = (x0, xi), where x0 = ct, and t is the coordinate time. The class of the harmonic coordinates is defined by imposing the de Donder gauge condition on the metric tensor (Fock, 1964; Weinberg, 1972),

α(-αβ)=0.    (1)

Einstein equations for the metric tensor are a complicated non-linear system of differential equations in partial derivatives. Because gravitational field of the solar system is weak and motion of matter is slow, we focus on the first post-Newtonian approximation of general relativity. We also assume that the source of the background metric is a perfect fluid. Under these assumptions the post-Newtonian background metric has the following form (Kopeikin et al., 2011)

00=-1+2V̄c2+2c4(Φ̄-V̄2)+O(c-6),    (2)
0i=-4V̄ic3+O(c-5),    (3)
ij=δij(1+2V̄c2)+O(c-4),    (4)

where the gravitational potentials entering the metric, satisfy the Poisson equations,

ΔV̄=-4πGρ̄,    (5)
ΔV̄i=-4πGρ̄vi,    (6)
ΔΦ̄=-4πGρ̄(2v2+2V̄+Π̄+3p̄ρ̄),    (7)

with p¯ and vi being pressure and velocity of matter respectively, and Π¯ is the internal energy of matter per unit mass. We emphasize that ρ¯ is the local mass density of baryons per a unit of invariant (3-dimensional) volume element dV=-u0d3x, where u0 is the time component of the 4-velocity of matter. The local mass density, ρ¯, relates in the post-Newtonian approximation to the invariant mass density ρ¯*=-u0ρ¯, namely (Kopeikin et al., 2011)

ρ̄*=ρ̄+ρ̄c2(12v2+3V̄).    (8)

The internal energy, Π¯, is related to pressure, p¯, and the local density, ρ¯, through the thermodynamic equation

dΠ̄+p̄d(1ρ̄)=0,    (9)

and the equation of state, p¯=p¯(ρ¯).

We shall further assume that the background matter rotates rigidly around fixed z axis with a constant angular velocity ω. This makes the background spacetime stationary with the background metric being independent on time. In the stationary spacetime, the mass density ρ* obeys the exact equation of continuity

i(ρ̄*vi)=0.    (10)

Velocity of rigidly rotating fluid is

vi=εijkωjxk,    (11)

where ωi = (0, 0, ω) is the constant angular velocity. Replacing velocity vi in Equation (10) with Equation (11), and differentiating, reveals that

viiρ̄=0,    (12)

which means that velocity of the fluid is tangent to the surfaces of a constant density ρ¯.

3. Relativistic Geoid

The Newtonian concept of Earth's geoid was extended to the post-Newtonian approximation of general relativity in Soffel (1989) and Kopejkin (1991). Exact concept of the relativistic geoid in general relativity that is not limited to the first post-Newtonian approximation has been discussed in Kopeikin et al. (2015) and Oltean et al. (2015).

We consider the Earth as made up of the background matter with small perturbations which reflect the actual mass distribution, stresses, and velocity flow of Earth's matter with respect to the rigidly rotating frame of reference. The angular velocity, ω, of Earth's rotation also slowly changes because of precession, nutation, polar motion and variations in the length-of-day. The physical metric, gαβ, depends on time and spatial coordinates. It is decomposed into an algebraic sum of the background metric ḡαβ, and its perturbation, ϰαβ,

gαβ=αβ+ϰαβ.    (13)

In the present paper, we shall neglect dependence of the perturbation ϰαβ on time because it produces very tiny relativistic effects that are currently unobservable.

Terrestrial reference frame is formed by the world lines of observers being fixed with respect to the crust of the Earth. Each observer moves in spacetime with four-velocity uα = c−1dxα where xα = {x0, x1, x2, x3} are coordinates of the observer, and τ is its proper time defined in terms of the metric tensor Equation (13) as follows,

c2dτ2=-gαβdxαdxβ.    (14)

Physical space of observers at each instant of time is represented by a three-dimensional hypersurface of constant proper time that is orthogonal everywhere to the world lines of the observers. The metric tensor on this hypersurface is denoted as hαβ and is given by Zelmanov et al. (2006) and Landau and Lifshitz (1975)

hαβgαβ+uαuβ.    (15)

It is used to measure the spatial distances. Rising and lowering of Greek indices of geometric objects residing on the perturbed manifold are done with the help of the full metric gαβ.

Similarly to classic geodesy, general relativity offers two definitions of relativistic geoid (Soffel, 1989; Kopejkin, 1991)

Definition 1: The relativistic u-geoid represents a two-dimensional surface at any point of which the rate of the proper time, τ, of an ideal clock carried out by a set of fiducial terrestrial observers is constant.

The u-geoid is determined by equation W = const., where the physical gravity potential

W(x)=c2(1-dτdt).    (16)

Because dτdt=-gαβαβ, where the dot means the total time derivative with respect to the coordinate time t, Equation (16) becomes,

W(x)c2[1-(-g00-2cg0ivi-1c2gijvivj)12].    (17)

This matches the Newtonian definition of the geoid after decomposition of the metric tensor in the post-Newtonian series and discarding all post-Newtonian terms.

Definition 2: The relativistic a-geoid represents a two-dimensional surface at any point of which the direction of a plumb line measured by the terrestrial observer, is orthogonal to the tangent plane of geoid's surface.

In order to derive equation of a-geoid, we notice that the direction of the plumb line is given by a four-vector of the physical acceleration of gravity, gα-c2aα where aα=uββuα is a four-acceleration of the terrestrial observer being orthogonal to uα. In case of a stationary rotating configuration of matter, we get (Lightman et al., 1975, Problems 10.14 and 16.17)

gα=-c2αln(1-Wc2),    (18)

where W is calculated in Equation (17). We consider an arbitrary displacement, dxαhαβdxβ, on the spatial hypersurface being orthogonal to uα everywhere, and make a scalar product of dxα with the direction of the plumb line. It gives,

dxαgα=dxαgα=-c2dln(1-Wc2).    (19)

From definition of a-geoid the left-hand side of Equation (19) must vanish due to the condition of orthogonality of the two vectors, dxα and gα. Therefore, it makes dln(1−Wc2) = 0, which means the constancy of the gravity potential W on the three-dimensional surface of the a-geoid. Thus, the surface of a-geoid coincides with that of u-geoid.

Calculation of the geoid is a difficult mathematical problem of geodesy. It is solved by introducing a reference surface of reference ellipsoid. Geoid's surface is defined in terms of the geoid undulation (height) with respect to the reference ellipsoid.

4. Post-Newtonian Reference Ellipsoid

In classical geodesy the reference figure for calculation of geoid's undulation is the Maclaurin ellipsoid which is a surface of the second order formed by a rigidly rotating fluid of constant density ρ. Maclaurin's ellipsoid is a surface of the second order (Chandrasekhar, 1969)

x2+y2a2+z2b2=1,    (20)

where a and b are semi-major and semi-minor axes of the ellipsoid. We also assume a > b, and define the eccentricity of the Maclaurin ellipsoid as

ea2-b2a2.    (21)

Physically, the ellipsoidal shape of rotating, self-gravitating fluid is formed because the Newtonian gravity potential is a scalar function represented by a polynomial of the second order with respect to the Cartesian spatial coordinates, and the differential Euler equation defining the equilibrium of the gravity and pressure is of the first order partial different equation which leads to the quadratic (w.r.t. the coordinates) equation of the level surface.

We shall demonstrate in the following sections that in the post-Newtonian approximation the gravity potential, W, of the rotating fluid is a polynomial of the fourth order as was first noticed by Chandrasekhar (1965). Hence, the post-Newtonian equation of the level surface of a rigidly-rotating fluid is expected to be a surface of the fourth order. We shall assume that the surface remains axisymmetric in the post-Newtonian approximation and dubbed the body with such a surface as PN-ellipsoid.

We shall denote all quantities taken on the surface of the PN-ellipsoid with a bar to distinguish them from the coordinates outside of the surface. Let the barred coordinates x¯i={x¯,ȳ,z¯} denote a point on the surface of the PN-ellipsoid with the axis of symmetry directed along the rotational axis and with the origin located at its post-Newtonian center of mass 1. Let the rotational axis coincide with the direction of z axis. Then, the most general equation of the PN-ellipsoid is

σ̄2a2+z̄2b2=1+κF(x̄),    (22)

where κ ≡ πGρa2c2, σ¯2x¯2+ȳ2, is the post-Newtonian parameter which is convenient in the calculations that follows,

F(x̄)K1σ̄2a2+K2z̄2b2+B1σ̄4a4+B2z̄4b4+B3σ̄2z̄2a2b2,    (23)

and K1, K2, B1, B2, B3 are arbitrary numerical coefficients.

Each cross-section of the PN-ellipsoid being orthogonal to the rotational axis, represents a circle. The equatorial cross-section has an equatorial radius, σ¯=re, being determined from Equation (22) by the condition z¯=0. It yields

re=a(1+12κB1).    (24)

The meridional cross-section of the PN-ellipsoid is no longer an ellipse (as it was in case of the Maclaurin ellipsoid) but a curve of the fourth order. Nonetheless, we can define the polar radius, z¯=rp, of the PN-ellipsoid by the condition, σ¯=0. Equation (22) yields

rp=b(1+12κB2).    (25)

The equatorial and polar radii of the PN-ellipsoid should be used in the post-Newtonian approximation instead of the equatorial and polar radii of the Maclaurin reference-ellipsoid for calculation of observable physical effects like the force of force, etc. We characterize the “oblateness” of the PN-ellipsoid by the post-Newtonian eccentricity

ϵre2-rp2re.    (26)

It differs from the eccentricity Equation (21) of the Maclaurin ellipsoid by relativistic correction

ϵ=e-κ1-e22e(B2-B1).    (27)

The PN-ellipsoid vs. Maclaurin's ellipsoid is visualized in Figure 1.


Figure 1. Meridional cross-section of the PN-ellipsoid (a red curve in on-line version) vs. the Maclaurin ellipsoid (a blue curve in on-line version). The left panel represents the most general case with arbitrary values of the shape parameters B1, B2, B3 when the equatorial, re, and polar, rp, radii of the PN-ellipsoid differ from the semi-major, a, and semi-minor, b, axes of the Maclaurin ellipsoid, rea, rpb. The right panel shows the most important physical case of B1 = B2 = 0 when the equatorial and polar radii of the PN-ellipsoid and the Maclaurin ellipsoid are equal. The angle φ is the geographic latitude (−90° ≤ φ ≤ 90°), and the angle θ is a complementary angle (co-latitude) used for calculation of integrals in appendix of the present paper (0 ≤ θ ≤ π). In general, when B1B2 ≠ 0, the maximal radial difference (the “height” difference) between the surfaces of the PN-ellipsoid and the Maclaurin ellipsoid depends on the choice of the post-Newtonian coordinates, and may amount to a few cm (see Section 13). Carefully operating with the residual gauge freedom of the post-Newtonian theory allows us to make the difference between the two surfaces much less than that.

Theoretical formalism for calculation of the post-Newtonian level surface can be worked out in arbitrary coordinates. For mathematical and historic reasons the most convenient are harmonic coordinates which are also used by the IAU (Soffel et al., 2003) and IERS (Petit and Luzum, 2010). The class of the harmonic coordinates is selected by the gauge condition Equation (1). Different harmonic coordinates are interrelated by coordinate transformations which are not violating the gauge condition 1. This freedom is known as a residual gauge (or coordinate) freedom. The field Equations (5)–(7) and their solutions are form-invariant with respect to the residual gauge transformations.

The residual gauge freedom is described by a post-Newtonian coordinate transformation,

xα=xα+κξα(x),    (28)

where functions, ξα, obey the Laplace equation,

Δξα=0.    (29)

Solution of the Laplace equation which is convergent at the origin of the coordinate system, is given in terms of the harmonic polynomials which are selected by the conditions imposed by the statement of the problem. In our case, the problem is to determine the shape of the PN-ellipsoid which has the surface described by the polynomial of the fourth order Equation (32) with yet unknown coefficients B1, B2, B3. The form of the Equation 32 does not change (in the post-Newtonian approximation) if the functions ξα in Equation (28) are polynomials of the third order. It is straightforward to show that the admissible harmonic polynomials of the third order have the following form

ξ1=hx+pxa2(σ2-4z2),    (30a)
ξ2=hy+pya2(σ2-4z2),    (30b)
ξ3=kz+qzb2(3σ2-2z2),    (30c)

where h, k, p, and q are arbitrary constant parameters. Polynomials Equations (30a)–(30c) represent solutions of the Laplace Equation (29). We choose ξ0 = 0 because we consider stationary spacetime so all functions are time-independent.

Coordinate transformation Equation (28) with ξi taken from Equations (30a)–(30c) does not violate the harmonic gauge condition Equation (1) but it changes the numerical post-Newtonian coefficients in the mathematical form of Equations (22) and Equation (23)

K1K1+2h,    (31a)
K2K2+2k,    (31b)
B1B1+2p,    (31c)
B2B2-4q,    (31d)
B3B3-8pb2a2+6qa2b2,    (31e)

Thus, it makes evident that only one out of the five coefficients K1, K2, B1, B2, B3 is algebraically independent while the four others can be chosen arbitrary. To simplify our calculations and eliminate the gauge-dependent terms from mathematical equations we choose K1 = K2 = B1 = B2 = 0. The constant B3B is left free. It will be fixed later on. With the accepted choice of the coordinates the polar radius rp = b, the equatorial radius re = a, the eccentricity of the PN-ellipsoid ϵ = e, and function

F(x̄)Bσ̄2z̄2a2b2.    (32)

Let xi = {x, y, z} is any point inside the PN-ellipsoid. We introduce a quadratic polynomial

C(x)σ2a2+z2b2-1.    (33)

This polynomial vanishes on the boundary surface of the PN-ellipsoid in the Newtonian approximation. However, in the post-Newtonian approximation we have on the boundary the following condition

C(x̄)=κF(x̄),    (34)

as a consequence of Equation (22). In terms of the polynomial C(x¯) function F(x¯) in Equation (32) can be formally recast to

F(x̄)=Bz̄2b2[1+C(x̄)-z̄2b2].    (35)

The term being proportional to C(x¯) can be discarded on the boundary of the PN-ellipsoid. Now we shall calculate the gravitational potentials of the rotating PN ellipsoid which enter the metric tensor Equations (2)–(4).

5. Newtonian Potential v¯

Newtonian gravitational potential V¯ satisfies the inhomogeneous Poisson's equation

ΔV̄(x)=-4πGρ̄(x),    (36)

inside the mass. Its particular solution is given by

V̄(x)=GVρ̄(x)d3x|x-x|,    (37)

where V is the volume occupied by the matter distribution. We shall assume that the density has an ellipsoidal distribution inside the body which differs from the homogeneous density, ρ0, only in the post-Newtonian approximation. Moreover, we shall consider the most simple, linearized case of the ellipsoidal distribution

ρ̄(x)=ρ0[1+κA(σ2a2+z2b2)],    (38)

where ρ0 is a constant bare density, and A is as yet undetermined numerical constant. This assumption allows us to perform integration in Equation (37) explicitly. Inside the mass the integral Equation (37) can be calculated by making use of spherical coordinates. The procedure is as follows (Chandrasekhar, 1969).

Let us consider a point xi = {x, y, z} inside the PN-ellipsoid Equation (22). It is connected to the point x¯i on the surface of the ellipsoid by vector Ri=x¯i-xi where Ri = Rℓi, R=δijRiRj, and a unit vector, ℓi = {sinθcosλ, sinθsinλ, cosθ}. We have

x̄i=xi+iR,    (39)

Substituting Equations (39)–(22) yields a quadratic equation

AR2+2BR+C=κF(x+R),    (40)


Asin2θa2+cos2θb2,  Bsinθ(xcosλ+ysinλ)a2+zcosθb2,    Cσ2a2+z2b2-1.    (41)

We solve Equation (40) by iterations by expanding R=R^+c-2ΔR, where R^ is either R^+ or R^- corresponding to two solutions of the quadratic Equation (40) with the right hand side being nil. The bare solution R^ is used, then, to calculate the right side of Equation (40) for performing the second iteration. We obtain

R±=-BA±B2-AC+κAF±A,    (42)


F±B(z+cosθR^±b)2-B(z+cosθR^±b)4,    (43)

where the term C(x¯)=AR2+2BR+C vanishes because of Equation (40).

We make replacement of variable x′ in Equation (37) to r = xx′, and use spherical coordinates to perform integration with respect to the radial coordinate r = |xx′| from the point r = 0 to the point lying on the surface of the PN ellipsoid, r = R(θ, λ) or r = R+(θ, λ). For the internal point, the angular integration in the remaining integral over the surface of the PN ellipsoid is performed over the solid angle 4π, and the integration with the point R(θ, λ) equals to that with the radial direction R+(θ, λ). This observation makes the angular integration easier because the integral Equation (37) can be written down in the following form (Chandrasekhar, 1969)

V̄=14Gρ0S2(R+2+R-2)dΩ+κAĪ1,    (44)

where R+ and R are defined in Equation (42), and the integral

Ī1=Gρ0Vd3x|x-x|(σ2a2+z2b2)    (45)

takes into account the ellipsoidal distribution of the density. After making use of Equation (42) and expanding the integrand in Equation (44) with respect to the post-Newtonian parameter κ, the Newtonian potential takes on the following form

V¯=12Gρ0S2{2B2ACA2+κ2A[F++F  BB2AC(F+F)]}dΩ+κAI¯1,    (46)

where all post-Newtonian terms of the higher order in κ have been discarded, and integration is performed over a unit sphere S2 with respect to the angles λ and θ, ≡ sinθdθdλ is the element of the solid angle.

Now, we expand F± in polynomial w.r.t. R±,

F±=α0 + α1cosθR± + α2(cosθR±)2 + α3(cosθR±)3    + α4(cosθR±)4,    (47)

where the coefficients

α0=Bz2b2(1-z2b2),    (48)
α1=2Bzb2(1-2z2b2)    (49)
α2=Bb2(1-6z2b2)    (50)
α3=-4Bzb4,    (51)
α4=-Bb4,    (52)

are the polynomials of z only. We also notice that on the surface of the PN-ellipsoid, F(x¯)=α0, as follows from Equation (35) where we can use, C(x¯)=0, in the post-Newtonian terms.

Replacing Equations (47) in (46) transforms it to

V̄=V̄N+κV̄pN+κAĪ1,    (53)


V̄N=Gρ0S2(B2A2-C2A)dΩ    (54)
V¯pN=Gρ0S2[α02Aα1cosθBA2+α2cos2θ(2B2A3C2A2)2α3cos3θ(2B3A4BCA3)+α4cos4θ(8B4A56B2CA4+C22A3)]dΩ,    (55)

Equations (52)–(55) describe the Newtonian potential exactly both on the surface of the PN-ellipsoid and inside it. These equations cannot be extended to the external space which requires a separate integration which will be discussed somewhere else.

The integrals entering V in Equations (54), (55) are discussed in Appendix A in Supplementary Material. The integrals I1 appears besides Equation (53) also in the post-Newtonian potentials of the gravitational field in Equation (80), and its discussion is postponed to section 7. After evaluating the integrals in the internal point and reducing similar terms, potentials VN and VpN take on the following form:

V̄N=πGρ0a2[(1-z2b2)0-(1-3z2b2)1-C(x)1],           (56)
V̄pN=πGρ0a2[F1(z)+b2F2(z)C(x)+b4F3(z)C2(x)],    (57)


F1(z)=α002α1z1+2α2b2[(1z2b2)1(13z2b2)2]4α3b2z[3(1z2b2)2(35z2b2)3]+6α4b4[(1z2b2)222(16z2b2+5z4b4)3+(110z2b2+353z4b4)4],    (58)
F2(z)=α2(122)4α3z(2233)+6α4b2[(1z2b2)23(13z2b2)3+2(15z2b2)4],    (59)
F3(z)=α4(263+64),    (60)

and the polynomial coefficients α0, α1, α2, α3, α4 are given in Equations (48)–(52). It is worth noticing that VpN obeys the Laplace equation

ΔV̄pN=0,    (61)

and Equation (57) represents the harmonic polynomial of the fourth order.

6. Post-newtonian Vector Potential v¯i

Vector potential Vi obeys the Poisson equation

ΔV̄i=-4πGρ̄(x)vi(x),    (62)

which has a particular solution

V̄i=GVρ̄(x)vi(x)|x-x|d3x.    (63)

For a rigidly rotating configuration, vi(x) = εijkωjxk so that

V̄i=εijkωjD̄k,    (64)


D¯i=Gρ¯(x)xid3x|x-x|.    (65)

It can be recast to the following form

Di=xiV¯N+Gρ¯(x)(xi-xi)|x-x|d3x,    (66)

where V¯N is the Newtonian potential given in Equation (56). Because the potential V¯i appears only in the post-Newtonian terms we can consider the density ρ¯ as constant, ρ¯=ρ0. In this case the second term in the right hand side of Equation (66) can be integrated over the radial coordinate, yielding

Vρ¯(x)(xi-xi)|x-x|d3x=ρ06S2(R+3+R-3)lidΩ.    (67)

After making use of Equation (42) to replace R+ and R, we obtain

Vρ¯(xi-xi)|x-x|d3x=ρ0S2(-43B3A3+BCA2)lidΩ,    (68)

where we have omitted the post-Newtonian terms (with κ) since the vector potential V¯i itself appears only in the post-Newtonian terms. Integrals entering Equation (68) are given in Appendix A in Supplementary Material. Calculation reveals

Vρ¯(x)(xx)|xx|d3x=πρ0a2x[(1z2b2)0       2(13z2b2)1+(15z2b2)2]       +xπρ0a2C(x)(12),    (69)
Vρ¯(x)(yy)|xx|d3x=πρ0a2y[(1z2b2)02(13z2b2)1+(15z2b2)2]+yπρ0a2C(x)(12),    (70)
Vρ¯(x)(zz)|xx|d3x=4πρ0a2z[(1z2b2)1(15z23b2)2]2zπρ0a2C(x)(122).    (71)

Substituting this result to Equation (66) and making use of Equation (56) yields

D¯i(D¯x,D¯y,D¯z)=(xD¯1,yD¯1,zD¯2),    (72)


D¯1πGρ0a2[(1-3z2b2)1-(1-5z2b2)2-C(x)2],    (73)
D¯2πGρ0a2[(1z2b2)0(57z2b2)1+4(15z23b2)2+C(x)(4231)],    (74)

7. Post-Newtonian Scalar Potential Φ¯

Potential Φ¯ is defined by equation

ΔΦ¯=-4πGρ¯(x)ϕ¯(x),    (75)

where function

ϕ¯(x)2ω2σ2+3p¯ρ¯+2V¯N.    (76)

Potential Φ¯ enters only the post-Newtonian equations. Therefore, everywhere in calculation of Φ¯ we can assume the density ρ¯ being approximated as constant. The pressure p¯ inside the massive body can be calculated from the hydrostatic equilibrium equation with a constant value of the bare density ρ¯=ρ0. The solution is Chandrasekhar (1969)

p¯ρ0=-πGρ0a2C(x)(0-21).    (77)

Making use of Equation (77) and Equation (56) we can write down function ϕ¯(x) as

ϕ¯(x)=a2[2ω2πGρ0(3041)](σ2a2+z2b2)2a2[ω2πGρ0(031)]z2b2+πGρ0a2(5061).    (78)

Particular solution of Equation (75) can be written as

Ф¯=2ω2a2(I¯1I¯2)πGρ0a2[(3I¯1+2I¯25V¯N)0+2(2I¯1+3I¯23V¯N)1],    (79)

where we have introduced two integrals

Ī1=Gρ0Vd3x|x-x|(σ2a2+z2b2),    (80)
Ī2=Gρ0b2Vz2|x-x|d3x.    (81)

Notice that the integral Ī1 has already appeared in the calculation of the Newtonian potential in Equation (45) while the integral Ī2 is a new one.

The integrals can be split in several algebraic pieces,

I¯1=(σ2a2+z2b2)(2D¯1V¯N)+2z2b2(D¯2D¯1)+Gρ08S2(R+4+R4)(sin2θa2+cos2θb2)dΩ,    (82)
I¯2=z2b2(2D¯2V¯N)+Gρ08b2S2(R+4+R4)cos2θdΩ,    (83)

where the integrals

18S2(R+4+R4)(sin2θa2+cos2θb2)dΩ=S2(2B4A32B2CA2+C24A)dΩ,    (84)
18S2(R+4+R4)cos2θdΩ=S2(2B4A42B2CA3+C24A2)cos2θdΩ,    (85)

We use the results of Appendix A in Supplementary Material to calculate these integrals, and obtain

S2(2B4A32B2CA2+C24A)dΩ=32πa2[(1z2b2)202(16z2b2+5z4b4)1+(110z2b2+353z4b4)2]+πa2[(1z2b2)04(13z2b2)1+3(15z2b2)2]C(x)πa2[1322]C2(x)    (86)
S2(2B4A42B2CA3+C24A2)cos2θb2dΩ=32πa2[(1z2b2)212(16z2b2+5z4b4)2+(110z2b2+353z4b4)3]+πa2(1z2b2)14(13z2b2)2+3(15z2b2)3]C(x)+πa2[2323]C2(x)    (87)

Substituting these results in Equations (82) and (83) yields

I¯1=12πGρ0a2[(1z4b4)06z2b2(153z2b2)1(110z2b2+353z4b4)2]πGρ0a2[3z2b21+(15z2b2)2+122C(x)]C(x),    (88)
I¯2=πGρ0a2[(1z2b2)z2b20+(3212z2b2+252z4b4)1(326z2b2+853z4b4)2+(3215z2b2+352z4b4)3]+πGρa2[(16z2b2)14(15z2b2)2+3(15z2b2)3+(2323)C(x)]C(x).    (89)

Replacing these expressions to Equation (79) results in

Φ¯=Φ0+Φ1C(x)+Φ2C2(x),    (90)


Φ0=12π2G2ρ02a4[70230(8152+23)+21(151202+93)]π2G2ρ02a4[702+0(601+672303)+21(511882+453)]z2b2+16π2G2ρ02a4[2102+0(2881+44522103)+101(5711162+633)]z4b4+πGρ0a4ω2[031+52332(091+212153)z2b2+(05(3192+73))z4b4],    (91)
Φ1=π2G2ρ02a4[0(71+11263)+21(61142+93)]+π2G2ρ02a4[0(211552+303)21(241702+453)]z2b22πGρ0a4ω2[(132+33)3(152+53)z2b2]    (92)
Φ2=12πGρ0a4[(02+812+6031813)πGρ06(23)ω2],    (93)

This finalizes the calculation of the post-Newtonian potentials inside the mass.

8. Post-Newtonian Mass

The post-Newtonian conservation laws have been discussed by a number of researchers the most notably in textbooks (Fock, 1964; Will, 1993; Kopeikin et al., 2011). General relativity predicts that the integrals of energy, linear momentum, angular momentum and the center of mass of an isolated system are conserved in the post-Newtonian approximation. In the present paper we deal with a single isolated body so that the integrals of the center of mass and the linear momentum are trivial and we can always choose the origin of the coordinate system at the center of mass of the body with the linear momentum being nil. The integrals of energy and angular momentum are less trivial and requires detailed calculations which are given below.

The law of conservation of energy yields the post-Newtonian mass of a rotating fluid ball that is defined as follows (Fock, 1964; Will, 1993; Kopeikin et al., 2011)

M=MN+1c2MpN,    (94)

where V is the volume of PN-ellipsoid, and

MN=Vρ¯(x)d3x,    (95)

is the rest mass of baryons comprising the body, and

MpN=Vρ¯(x)(v2+Π+52VN)d3x,    (96)

is the post-Newtonian correction.

In order to calculate the rest mass, MN, we introduce normalized spherical coordinates r, θ, ϕ related to the Cartesian (harmonic) coordinates x, y, z as follows,

x=arsinθcosϕ,  y=arsinθsinϕ,  z=brcosθ.    (97)

In these coordinates the integral Equation (95) is given by

MN=a2bρ00r(θ)0π02π(1+κAr2)r2sinθdrdθdϕ,    (98)

where r(θ) describes the surface of the PN-ellipsoid defined above in Equations (22), (32)

r2(θ)=1+κBsin2cos2θ.    (99)

Integration in Equation (98) results in

MN=4π3ρ0a2b[1+κ5(A+B)].    (100)

The post-Newtonian contribution, MpN, to the rest mass reads

MpN=ρ0a4b010π02π{ω2r2sin2θ+52πGρ0[0(1r2cos2θ)1r2(13cos2θ)]}r2sinθdrdθdϕ=8π15ρ0a4b(ω2+5πGρ00).    (101)

The total mass Equation (94) becomes

M=4π3ρ0a2b[1+κ5(2ω2πGρ0+A+B+100)].    (102)

The inverse relation is used to convert the density ρ0 to the total mass,

ρ0=3M4πa2b[1-κ5(A+B+100+2ω2πGρ0)].    (103)

The total post-Newtonian mass, M, depends on the sum of the parameters, A+B.

9. Post-Newtonian Angular Momentum

Vector of the post-Newtonian angular momentum, Si = (Sx, Sy, Sz), is defined by Fock (1964), Will (1993), and Kopeikin et al. (2011)

Si=SNi+1c2SpNi,    (104)

where SNi and SpNi are the Newtonian and post-Newtonian contributions respectively,

SNi=ρ¯(x)(x×v)id3x,    (105)
SpNi=ρ¯(x)(v2+6V¯+p¯ρ¯)(x×v)id3x4ρ¯(x)(x×V¯)id3x,    (106)

and vector-potential V¯V¯i has been given in Equations (64) and (72).

It can be checked up by inspection that in case of an axisymmetric mass distribution with the ellipsoidal density distribution Equation (38), the only non-vanishing component of the angular momentum is S3 = SzS. Indeed, v = (ω × x), and (x × v)i = (x × (ω × x))i = ωi(x2 + y2 + z2) − xi ω z. It results in

SNx=ωρ¯xzd3x,SNy=ωρ¯yzd3x,SNz=ωρ¯(x2+y2)d3x.    (107)

By making use of the coordinates Equation (97) the reader can easily confirm that SNx=SNy=0, and SNzSN reads

SN=a4bω0r(θ)0π02πρ¯(x)r4sin3θdrdθdϕ,    (108)

where r(θ) is defined in Equation (99). After integration we obtain,

SN=8π15a4bρ0ω[1+κ7(5A+2B)].    (109)

Replacing the density ρ0 by the total mass M with the help of Equation (102), makes the Newtonian part of the angular momentum take on the following form,

SN=25Ma2ω[1-κ35(2A-3B+700+14ω2πGρ0)].    (110)

It is straightforward to prove that x and y components of SpNi also vanish due to the axial symmetry, and only its z component, SpNzSpN, remains. We notice that

ρ¯(x)(x×V¯)zd3x=ρ¯(x)D¯1(x×v)zd3x,    (111)

where D1 is taken from Equation (73). Therefore,

SpN=ρ¯(x)(v2+6V¯N+p¯ρ¯-4D¯1)(x×v)zd3x.    (112)

Making transformation to the coordinates Equation (97) and approximating ρ¯=ρ0, yields

SpN=ωa6bρ0010π02π{ω2r2sin2θ+7061(0+4142)r22(30151+102)r2cos2θ}r4sin3θdrdθdϕ=435Ma2ω[2ω2+πGρ0(190161)].    (113)

The total angular momentum becomes

S=25Ma2ω{1-κ35[2A-3B-40(30-41)-6ω2πGρ0]}.    (114)

10. Post-Newtonian Equation of the Level Surface

The figure of the rotating fluid body is defined by the boundary condition of vanishing pressure, p = 0. The surface p = 0 is called the level surface. Relativistic Euler equation derived for the rigidly rotating fluid body tells us that the level surface coincides with the equipotential surface of the post-Newtonian gravitational potential W¯ which is expressed in terms of the centrifugal and gravitational potentials by the following expression (Kopejkin, 1991; Kopeikin et al., 2011, 2015)

W¯=12ω2σ2+V¯N+κ(V¯pN+AI¯1)+1c2(18ω4σ4+32ω2σ2V¯N4ω2σ2D¯112V¯N2+Φ¯),    (115)

where κπGρ0a2c2, and the potentials V¯N,V¯pN,Ī1,D¯1,Φ¯ have been explained in Sections 5–7. After substituting these potentials to equation Equation (115) it can be presented as a quadratic polynomial with respect to function C(x),

W¯(x)=W¯0+W¯1C(x)+W¯2C2(x),    (116)

where the coefficients of the expansion are the polynomials of the z coordinate only. In particular, the coefficient W¯0 is a polynomial of the fourth order,

W¯0=K0+K1z2b2+K2z4b4,    (117)


K0=12ω2a2+πGρ0a2(01)+18c2ω4a4+12κω2a2(50171+18263)+12κπGρ0a2[6020(221152+63)+1(291402+183)]+κπGρ0a2[2(142+6334)B+12(02)A],    (118)
K1=12ω2a2πGρ0a2(031)14c2ω4a4κω2a2(50401+662303)κπGρ0a2(6025601+9912+6702176123003+9013)+κπGρ0a2[(0181+7821203+604)B(3152)A],    (119)
K2=18c2ω4a4+12κω2a2(50631+1302703)+16κπGρ0a2[180250(541892+423)+1(543111602+6303)]κπGρ0a2[(0201+9021403+704)B+12(0101+3532)A].    (120)

Coefficient W¯1 in Equation (116) is a polynomial of the second order,

W¯1=P+P1z2b2,    (121)


P=12ω2a2[ 1+12c2ω2a2+κ(30181+282123) ]πGρ0a21+κπGρ0a2[(182+183124)B2A]κπGρ0a2(60111121102+2812+6031813),    (122)
P1=14c2ω2a432κω2a2(0161+362203)+κπGρ0a2[200145125502+14012+30039013]κπGρ0a2[2(31252+513304)B+(3152)A].    (123)

Coefficient W¯2 in Equation (116) is constant,

W¯2=18c2ω4a412κω2a2(312263)κπGρ0a2 [ (263+64)B+122A+12(12+02812603+1813) ].    (124)

On the level surface of the PN-ellipsoid2 we have all three coordinates interconnected by Equation (34) of the PN-ellipsoid, C(x¯)=κα0(z¯), so that Equation (116) becomes

W¯W¯0+κW¯1α0(z¯),    (125)

and the term with W¯2~O(κ2), is discarded as negligibly small. After reducing similar terms, the potential W¯ on the level surface is simplified to the polynomial of the fourth order,

W¯=K¯0+K¯1z¯2b2+K¯2z¯4b4,    (126)


K¯0=K0,    (127)
K¯1=K1+κ(12ω2a2-πGρ0a21)B,    (128)
K¯2=K2-κ(12ω2a2-πGρ0a21)B.    (129)

Because the potential W¯ is to be constant on the level surface (Kopeikin et al., 2015), the numerical coefficients K¯1 and K¯2 must vanish. The first condition, K¯1=0, yields a relation between the angular velocity of rotation, ω, and the eccentricity, e, of the rotating fluid ellipsoid:

ω22πGρ0[1+ω2a22c2+2κ(50401+662303B)] =310κ[6025601+9912+670217612 3003+9013+(3152)A(0181+782 1203+604)B].    (130)

Equation (130) generalizes the famous result that was first obtained by Colin Maclaurin in 1742, from the Newtonian theory of gravity to the realm of general relativity.

The second condition, K¯2=0, yields an algebraic equation for the two coefficients A, B,

12(0101+3532)A+(ω22πGρ0+0211+9021403+704)B=ω48π2G2ρ2+ω24πGρ(50631+1302703)+16[180250(541892+423)+1(543111602+6303)].    (131)

Equation (131) imposes one constraint on two coefficients A, B defining the shape of the PN-ellipsoid and the law of distribution of mass density. One more constraint is required to fix these coefficients. We can use, for example, the Maclaurin relation Equation (130) to set an additional (geodetic) constraint on the parameters by demanding that the post-Newtonian part of Equation (130) vanishes.

We explored another possibility to impose self-consistent constraints on the shape parameters of the PN-ellpsoid directly related to the gravimetric measurements in geodesy. Newtonian theory of rotating reference-ellipsoid connects the shape parameters of the Maclaurin ellipsoid with the measured values of its gravity force at the pole and equator of the ellipsoid through the Pizzetti and Clairaut theorems (Pizzetti, 1913; Torge and Müller, 2012). Let us denote the force of gravity measured at the pole of the ellipsoid by γp ≡ γi(x = 0, y = 0, z = b), and the force of gravity measured on equator by γe ≡ γi(x = 0, y = a, z = 0). Due to the rotational symmetry of the ellipsoid the equatorial point can be, in fact, chosen arbitrary. The classic form of the Pizetti theorem is

2γea+γpb=3GMa2b2ω2,    (132)

while the Clairaut theorem states

γeaγpb=3GM2a2b3ee331e2arcsinee3+ω2.    (133)

These two theorems were crucial in geodesy of XIX-th century for deeper understanding that the geometric shape of Earth's figure can be determined not only from the measurements of the geodetic arcs but, independently, from rendering the intrinsic measurements of the force of gravity on its surface. The gravity-geometry correspondence led to the pioneering idea that the gravity force and geometry of space must be always interrelated—the idea which paved a way in XX-th century to the development of the general theory of relativity by A. Einstein. We derive the post-Newtonian analogs of the Pizzetti and Clairaut theorems in the next two sections and explore what kind of natural limitations (if any) they can set on the shape parameters of the post-Newtonian ellipsoid.

11. Post-Newtonian Pizzetti's Theorem

The force of gravity on the equipotential surface of a massive rotating body is defined by equation (Kopeikin et al., 2011)

γi=[ΛjijW¯]x=x¯,    (134)

where ixi, the post-Newtonian gravity potential W¯ is given in Equation (116), and

Λji=δij(11c2V¯N)12c2vivj,    (135)

is the matrix of transformation from the global coordinates to the local inertial frame of observer, vi = (ω × x)i is velocity of the observer with respect to the global coordinates, and VN is the Newtonian potential Equation (56). We emphasize that we, first, take a derivative in Equation (134), and then, take the spatial coordinates, x, on the equipotential surface described by the equation of the post-Newtonian ellipsoid Equation (22).

Velocity vi is orthogonal to the gradient iW¯ everywhere, that is

viiW¯=0.    (136)

Indeed, it is easy to prove that

viiW¯=ω(xyW¯yxW¯).    (137)

Partial derivatives of W are calculated from Equation (116),

xW¯=dW¯dCxC(x)=dW¯dC2xa2,yW¯=dW¯dCyC(x)=dW¯dC2ya2.    (138)

Replacing the partial derivative from Equations (138) to (137) yields Equation (136). After accounting for Equations (136), (134) is simplified to

γi(x¯)=[(11c2V¯N)iW¯]x=x¯.    (139)

We shall denote the force of gravity on the pole by γp ≡ γi(x = 0, y = 0, z = rp), and the force of gravity on the equator by γe ≡ γi(x = 0, y = re, z = 0) where the equatorial re and polar rp radii are defined in Equations (24) and (25) respectively. Taking the partial derivative from W¯ in Equation (139), yields

γp=2πGρ0a2b(021)+16ω2a2bκ(132+23)          +2πGρ0a2bκ[(0171+602763+324)B          +(041+832)A]+4πGρ0a23bκ[0(271562          +243)21(331742+363)],    (140)
γe=a(21Gπρ0ω2)ω4a32c2+κ[30+181282        +123+(01)]ω2a2κπGρ0a[(182+183        124)B+2A]+2κπGρ0a[0(51112+63)        21(51142+93)].    (141)

We form a linear combination

2γea+γpb=2πGρ0[21+a2b2(021)]2ω21c2ω4a2                  2κ[20171+2821238a2b2(1                  32+23)]ω22κπGρ0[21162+363                  244a2b2(0171+602763+324)]B                  2κπGρ0[22+a2b2(041+832)]A                  +4κπGρ0[0(51112+63)21(51                  142+93)]+4a23b2κπGρ0[0(271562                  +243)21(331742+363)].    (142)

Making use of integrals given in Appendix A in Supplementary Material, we can check that Equation (142) is simplified to

2γea+γpb=4πGρ02ω21c2ω4a2               κ3e7[e(1e2)(105104e2+42e4)               31e2(54e2)(76e2+2e4)arcsine]ω2κπGρ012e9(74e2)[5e(2131e2+10e4)31e2(3540e2+8e4)arcsine]B+κπGρ0e5[e(7193e22e4)1e2(74e2)arcsine]A+κπGρ03e10(1e2)[315e2+621e4250e6+24e8+2e1e2(315516e2+169e418e6)arcsine3(105242e2+178e440e6)arcsin2e].    (143)

In order to compare Equation (143) with its classic counterpart Equation (132), we convert, then, the density ρ0, to the total relativistic mass of the PN-ellipsoid by making use of Equation (103). It allows us to recast Equation (143) to

2γea+γpb=3GMa2b2ω2a2ω4c2+31e216e9a2b2G2M2c2[e1e2(315+621e2250e4+24e6)+2(315831e2+685e4187e6+42e8)arcsine]+120e7bGMω2c2[e(525+1045e2730e4+234e6)+153593e2+92e442e6+8e81e2arcsine]916G2M2c2105347e2+420e4218e6+40e8e10a2b2arcsin2e3320e8a2b2G2M2c2e[36757525e2+4850e41000e648e8151e2e(74e2)(3540e2+8e4)arcsine]B+916G2M2a2b2e5c2[e(7193e22e4)1e2(74e2)arcsine]A.    (144)

This equation represents the post-Newtonian extension of the classical Pizzetti theorem (132).

12. Post-Newtonian Clairaut's Theorem

In order to derive the post-Newtonian analog of the Clairaut theorem Equation (133) we follow its classic derivation given, for example, in Pizzetti (1913). To this end we form a linear difference between the forces of gravity measured at the pole and on the equator of the PN-ellipsoid by making use of EquationS (139) and (140). We get

γeaγpb=2πGρ0[1+a2b2(021)]+ω2+12c2ω4a2+κ[20171+282123+16a2b2(132+23)]ω2κπGρ0[21+162363+2442a2b2(0171+602763+324ρ)]B2κπGρ0[22a2b2(041+832ρ)]A+2κπGρ0[0(51+11263)+21(51142+93)]+4a23b2κπGρ0[0(271562+243)21(331742+363)].    (145)

We use the results of Appendix A in Supplementary Material to replace the integrals entering the right side of Equation (145), with their explicit expressions given in terms of the eccentricity e of the Maclaurin ellipsoid Equation (20). It yields

γeaγpb=6πGρ0e3(e1e2arcsine)2πGρ+ω2+12c2ω4a2+κ6e7(75e5e3+122e542e7+31e2(25+10e234e4+8e6)arcsine)ω2κ24e9πGρ0[e(5251075e2+662e4112e6)31e2(175300e2+144e416e6)arcsine]Bκ6e5πGρ0[e(1529e2+6e4)31e2(58e2)arcsine]Aκ6e9πGρ0(1e2)[e(225351e258e4+24e6)21e2(225276e285e4+18e6)arcsine]κ2e10πGρ0(1e2)(75142e2+38e4+40e6)arcsin2e.    (146)

This form of the Clairaut theorem apparently depends on the shape parameters A, B in different combinations and can be used to impose a constraints on one of them in addition to the constraint given by the level surface Equation (131). For example, we could demand that all post-Newtonian terms in Equation (146) vanish.

The last step is to replace the density ρ0 in Equation (146) with the total mass of the PN-ellipsoid by making use of the gauge-invariant expression Equation (103). We get

γeaγpb=3GM2a2b3ee331e2arcsinee3+ω2+a2ω42c2e(375+25e2682e4+234e6)31e2(125+50e2194e4+40e6)arcsine40be7GMc2ω23640e(26255375e2+3310e4704e6+48e8)1e2(8751500e2+720e4128e6)arcsinea3be9G2M2c2B332e(1529e2+6e4)31e2(58e2)arcsinea2b2e5G2M2c2A3160e1e2(225351e258e4+24e6)2(225501e2+191e4+175e642e8)arcsinea4e9G2M2c2932G2M2c275142e2+38e4+88e6a4e10arcsin2e    (147)

This is the post-Newtonian extension of the classical Clairaut theorem Equation (133).

13. Reference Ellipsoid in Small-Eccentricity Approximation

Let us apply the formalism of previous sections to derive practically useful post-Newtonian relationships used for processing geodetic and gravimetric measurements on the surface of Earth and in space. To this end we notice that the flattening, f of the terrestrial reference ellipsoid is about fe2 ∕ 2 ≃ 1 ∕ 298 = 0.0034 (Petit and Luzum, 2010, Section 1) that can be used as a small parameter for expanding all post-Newtonian terms in the convergent Taylor series with respect to f. We shall keep the Newtonian expressions as they are without expansion, and take into account only terms of the order of e2 in the post-Newtonian parts of equations by systematically discarding terms of the order of ~e4, and higher. Because, according to Maclaurin's relation Equation (130), the square of the angular velocity, ω2e2, (see Pizzetti, 1913; Chandrasekhar, 1969; Torge and Müller, 2012 for more detail), we shall also discard terms of the order of ω4 and ω2e2. This procedure shall dubbed as a small-eccentricity approximation.

Now, we have to make a decision about what kind of constraints on the model parameters of the PN ellipsoid A, B would be more preferable for practical applications. We have considered the case of A = 0, B ≠ 0 in our previous work (Kopeikin et al., 2016). We have proved that BO(e4), and can be ignored in all geodetic equations of the small-eccentricity approximation under assumption of the homogeneous density distribution. In the present paper we explore the complementary case A ≠ 0, B = 0. This constraint makes the level surface of the post-Newtonian ellipsoid coinciding exactly with the classic ellipsoid Equation (20) which means that the ellipsoidal coordinate system, having been ubiquitously used in geodesy, is not deformed by the relativistic corrections to the gravity field at all. This requires to adopt that the mass density inside the ellipsoid has an inhomogeneous ellipsoidal distribution given by Equation (38) depending on the parameter A which is determined from the equation of the level surface Equation (131). We solve Equation (131) under this constraint with respect to A, and find out that it can be approximated with

8315e4A=4724725e4+11152363825e6+O(e8),    (148)

that yields

A=5915+1394e21155,    (149)

Post-Newtonian correction to the mass of the Earth, M, can be evaluated from Equation (102). Because the mass couples with the universal gravitational constant G, it contributes to the numerical value of the geocentric gravitational constant GM=3.986004418×1014 m3s−2. Numerical values for the geopotential and the semi-major axis, a, of reference ellipsoid are given in Petit and Luzum, 2010, Table 1.1.) They yield, GMc2a=6.7×10-10. After expansion of the right side of Equation (102) with respect to the eccentricity e, the relativistic variation in the value of GM is

δ(GM)pNGM2κ03GMc2a2.08×109,    (150)

where the numerical value of κ was calculated on the basis of relationship κ=πGρ0a2(34)GMc2a. The current uncertainty in the numerical value of GM is δGM=8×105 m3s−2 (Petit and Luzum, 2010, Table 1.1) which gives the fractional uncertainty

δ(GM)GM2.0×109.    (151)

This is comparable with the relativistic contribution Equation (150) that must be taken into account in the reduction of precise geodetic data processing which are currently based on the Newtonian theory.

The Maclaurin relation Equation (130) takes on the following approximate form

ω22πGρ0=1e2(32e2)arcsine3e(1e2)e3+99175GMac2e2,    (152)

which, after replacing ρ0 with the help of Equation (103) and expansion, can be recast to

ω2=3GM2a3[(32e2)arcsine3e1e2e3+99175GMac2e2].    (153)

This can be further decomposed and presented as follows

ω=ω(1+297280GMac2),    (154)


ω=(3GM2a3(32e2)arcsine3e1e2e3)1/2,    (155)

is the currently used value for the Earth angular velocity of rotation, ω=729211.5×10-10 rad·s−1 (Petit and Luzum (2010). Contribution from the post-Newtonian term in the right side of Equation (154) amounts to 7.1 × 10−10 which makes it evident that the currently adopted value for ω should be corrected by the IUGG.

The approximate version of the post-Newtonian Pizzetti theorem Equation (144) reads

2γea+γpb=2ω2+3GMa2b+4GMac2[(3+4725e2)GMa3+14ω2].    (156)

The post-Newtonian corrections in the Clairaut theorem Equation (147), after they are expanded with respect to the eccentricity, yield

γeaγpb=ω2+3GM2a2be(3e2)31e2arcsinee3                    +GMac2(5925GMa3e2+1114ω2).    (157)

The post-Newtonian corrections to the gravitational field entering the the Pizzetti and Clairaut theorems Equations (156), (226) are not so negligibly small, amount to the magnitude of approximately 3 μGal (1 Gal = 1 cm/s2), and are to be taken into account in calculation of the parameters of the reference-ellipsoid from astronomical and gravimetric data in a foreseeable future.

Equations (148)–(151) extrapolates the Newtonian relations adopted by IERS and IUGG for definition of the Earth reference ellipsoid and for connecting its parameters to the angular velocity ω and gravity field, to the realm of general relativity.

14. Post-Newtonian Geoid's Undulation

We define the anomalous (disturbing) gravity potential TT(x) as the difference between the physical gravity potential, WW(x), and the background gravity potential, W¯(x), of the PN reference ellipsoid

T=WW¯.    (158)

The gravity potential W is defined by the Einstein equations with the sources that take into account the real distribution of density, pressure, etc. of the real Earth. We don't need to solve the full Einstein equations again in order to find out T. Because the difference between W and W¯ is very small, we can consider merely a linear perturbation of the background metric tensor [see Equation (13)]

ϰαβ=gαβ(x)g¯αβ(x).    (159)

The gravity potential W is constant on the perturbed equipotential surface, and is defined by

W=c2(1dτdt),    (160)

where dτ=c-1-gαβdxαdxβ is the proper time of observer in the real physical spacetime. Making use of Equations (160) and (16) allows us to recast Equation (158) to

T=c2(dτ¯dtdτdt),    (161)

which can be further simplified by noticing that

(dτdt)2=(1u¯αu¯βϰαβ)(dτ¯dt)2,    (162)

where ūα is the unperturbed four-velocity of observer which coincides with the four-velocity of matter of the rotating reference ellipsoid. Accounting for definition Equation (161), we get the anomalous gravity potential in the form,

T=c22u¯αu¯βϰαβ,    (163)

which is fully sufficient for practical applications. We emphasize that ϰαβ is the difference between the actual physical metric, gαβ, and the metric ḡαβ of the background manifold, and as such, should not be confused with the post-Newtonian expansion of the metric gαβ around a flat spacetime with the Minkowski metric ηαβ = diag(−1, 1, 1, 1) like it is shown in Equations (2)–(4). Our next task is to derive the differential equation for the anomalous gravity potential T.

Let us assume that inside the Earth the deviation of the real matter distribution from its unperturbed value is described by the symmetric energy-momentum tensor

c2Tαβ=euαuβ+sαβ,    (164)

where uα is four-velocity, 𝔢 is the energy density, and 𝔰αβ is the symmetric stress tensor of the perturbing matter. The stress tensor includes the isotropic pressure (diagonal components) and shear (off-diagonal components), and is orthogonal to uα, that is 𝔰αβuα=0. The energy density of the matter perturbation

e=ρ(c2+P),    (165)

where ρ is the mass density of the perturbation, and 𝔓 is the internal (compression) energy of the perturbation. For further calculations, a more convenient metric variable is

lαβϰαβ+12g¯αβϰ,    (166)

where ϰαβϰαβ. The dynamic field theory of manifold perturbations leads to the following equation for lαβ (Kopeikin and Petrov, 2013, 2014),

lαβ|μ|μ+g¯αβAμ|μ2Aα|βR¯μαlβμR¯μβlαμ2R¯αμνβlμν+2Fαβm=16πGc4Tαβ,    (167)

where Aαlαβ|β is the gauge vector function, depending on the choice of the coordinates, R¯αμνβ is the Riemann (curvature) tensor of the background manifold depending on the metric tensor ḡαβ, its first and second derivatives, R¯αβ=μνR¯μανβ—the Ricci tensor, and Fαβm is the tensorial perturbation of the background matter induced by the presence of the perturbation 𝔗αβ (see Kopeikin and Petrov, 2013, Equations 148–150 for particular details).

In what follows, we focus on derivation of the master equation for the anomalous gravity potential T in the exterior space that is outside of the background matter of the reference ellipsoid. This assumption is normally used in geodesy (Torge and Müller, 2012). To achieve our goal, we introduce two auxiliary scalars,

qu¯αu¯βlαβ+l2,    (168)
p(g¯αβ+u¯αu¯β)lαβ,    (169)


lg¯aβlαβ=2(pq).    (170)

In terms of the scalar 𝔮 the anomalous gravity potential Equation (163) reads

T=c22q,    (171)

where we have used the property ϰ = l. Taking from both sides of Equation (171) the covariant Laplace operator yields

Tc22q,    (172)

where qq|μ|μ is to be calculated from Equation (167).

We notice that according to Kopeikin and Petrov (2013, Equations 148–150) Fαβm is directly proportional to the thermodynamic quantities of the background matter and, thus, vanishes in the exterior (with respect to the background matter) space. Hence, we can drop off Fαβm in (237) in the exterior-to-matter domain. After contracting Equation (167) with ḡαβ, and accounting for Equation (170) we obtain

qp|α|αAα|α=8πGc2T,    (173)

where all terms depending on the Ricci tensor R¯aβ cancel out, 𝔗αβ𝔗αβ, and we still have terms with the gauge field Aα. Now, we use the gauge freedom of general relativity to simplify (246). More specifically, we impose the gauge condition

Aα=pα,    (174)

where 𝔭ααββ𝔭. This gauge allows us to eliminate function 𝔭 from Equation (173) and, after making use of Equation (172), reduce it to to a simple form of a covariant Poisson equation

T=4πGT.    (175)

In the Newtonian approximation the trace of the energy-momentum tensor is reduced to the negative value of the matter density of the perturbation, 𝔗≃−ρ, while □T≃ΔNT. Hence, Equation (175) matches its Newtonian counterpart. Outside the mass distribution the master equation for the anomalous gravity potential is reduced to the covariant Laplace equation

T=0.    (176)

Equations (175), (176) extend equations of classic geodesy to the realm of the post-Newtonian chronometric geodesy. The main difference is that the covariant Laplace operator in Equations (175), (176) is taken in curved space with the metric ḡαβ. The explicit form of the covariant Laplace operator applied to a scalar T in coordinates xi = {x1, x2, x3}, reads (Lightman et al., 1975, Problem 7.7.)

T1g¯i(g¯g¯ijjT),    (177)

where the repeated indices mean the Einstein summation, ixi is the partial derivative with respect to the spatial coordinates, and we omitted all time derivatives because the background spacetime is stationary. For the geodetic purposes it is convenient to choose a spherical coordinate system, {x1, x2, x3} = {r, θ, λ} co-rotating rigidly along with the reference ellipsoid and with z-axis directed along the axis of the rotation. With a sufficient post-Newtonian accuracy the metric tensor of the background spacetime in this coordinate system reads

g¯00=1+1c2(ω2r2sin2θ+2V¯N)+O(c3),    (178)
g¯0λ=1cωr2sin2θ+O(c3),    (179)
g¯rr=1+2V¯Nc2+O(c3),    (180)
g¯θθ=(1+2V¯Nc2)r2+O(c3),    (181)
g¯λλ=(1+2V¯Nc2)r2sin2θ+O(c3),    (182)

where V¯ is the Newtonian gravitational potential defined in Equation (56). Determinant of the metric

g¯det[g¯αβ]=r4sin2θ(1+2V¯c2)2,     (183)

After calculating the inverse metric ḡαβ, substituting it to Equation (177), and discarding all terms od the order of O(c−3), the equation for the anomalous gravity potential takes on the following form

ΔT=ω2c22Tλ2,    (184)


Δrr+2rr+1r2θθ+1r2tanθθ+1r2sin2θλλ,    (185)

is the Laplace operator in the spherical coordinates. The post-Newtonian equation Equation (184) can be solved by iterations by expanding the distrubing potential in the post-Newtonian series

T=TN+1c2TpN+O(c4),    (186)

where TN is the Newtonian disturbing potential expressed, for example, by the Stokes integral formula (Heiskanen and Moritz, 1967, Equations 2–163a), and TpN is the post-Newtonian correction which we are searching for. Substituting Equation (186) in Equation (184) yields the inhomogeneous differential equation for TpN,

ΔTpN=ω22TNλ2,    (187)

which can be solved by means of standard mathematical techniques for known TN.

We introduce the relativistic geoid height, N, by making use of relativistic generalization of Bruns' formula (Heiskanen and Moritz, 1967, Equations 2–144). We consider a set of equipotential surfaces defined by constant values of the gravity potential W. Let a point Q lie on an equipotential reference surface S1 and has coordinates xQα, and a point P lie on another equipotential surface S2, and has coordinates xPα. The height difference, N, between the two surfaces is defined as the absolute value of the integral taken along the direction of the plumb line passing through the points Q and P,

N=QPnαdxαdd,    (188)

where nαgα ∕ |g| is the unit (co)vector along the plumb line, gα is the relativistic acceleration of gravity Equation (18), |g|(h¯aβgαgβ)12, and ℓ is the proper length defined in space by Zelmanov et al. (2006) and Landau and Lifshitz (1975)

d2=h¯αβdxαdxβ,    (189)

where h¯αβαβ+ūαūβ is the metric tensor in three-dimensional space. In case, when the height difference is small enough, we can use the second mean value theorem for integration (Hobson, 1909) and approximate the integral in Equation (188) as follows

N=QPgα(x)dxα|g(x)|=c2gQQPαln(1Wc2)dxα                           =c2gQln|1W(Q)c21W(P)c2|,    (190)

where gQ = |g(Q)| denotes the magnitude of the relativistic acceleration of gravity taken on the equipotential surface S1. Equation (190) is exact. Separation of the height N in the Newtonian part and the post-Newtonian corrections depends on how we define the reference equipotential surface S1.

Let us choose the reference surface by equation W(Q)=W¯ where W¯ is the normal (relativistic) gravity field of the reference ellipsoid. Then, expanding the logarithm in Equation (190) with respect to the ratio Wc2 and making use of definition Equation (158) of the anomalous gravity potential T, we obtain from (190)

N=|T(P)|γQ,    (191)

where the disturbing potential, T, is measured at the point P on the geoid surface W, and the acceleration of gravity γQgQ is measured at point Q on the reference surface W¯.

Relativistic Bruns' formula Equation (191) yields geoid's undulation with respect to the unperturbed reference level surface, W¯=const., in general relativity. Because we have defined this surface as an equipotential surface W¯ which is the post-Newtonian solution of the Einstein equations, the height N does not represent the undulation of the relativistic geoid W with respect to the Newtonian equipotential surface defined by the surface of the constant Newtonian gravity potential, V¯N=const. Expansion of the height N in Equation (191) in the post-Newtonian series around the value of the surface W¯, yields

N=N+1c2NpN+O(c4),    (192)

where 𝔑 is the classic definition of the geoid height given in terms of the Newtonian disturbing potential, VN-V¯N as explained, for example, in textbook (Vaníček and Krakiwsky, 1986, chapter 21). The post-Newtonian correction 𝔑pN to the height 𝔑 is caused by the difference between the post-Newtonian terms in W and W¯, and has a magnitude of the order 𝔑pN(VNc2)×𝔑, where VN is the Newtonian gravitational potential of the Earth. Because the largest undulation of the Newtonian geoid of the Earth does not exceed 100 m (Torge and Müller, 2012), the explicit post-Newtonian correction 𝔑pN to the undulation is exceedingly small, 𝔑pN7×10-6 cm.

Thus, the main relativistic effects in terrestrial geodesy enter through the difference between the equipotential surface of the relativistic potential of the gravity force, W¯, and one of the Newtonian gravity potential VN. This difference amounts to (VNc2)×R1 cm which looks pretty small but crucially important in the current study of the changes in the global average sea level which is now rising twice as fast as it did over the past century, providing clear evidence of global warming on a short time scale. The measurements are done with the help of the the satellite laser altimetery which is so precise that allows us to measure the change in the global average sea level with the uncertainty ≤ 1 mm/year (Fu and Haines, 2013; Ablain et al., 2015) which is a factor of 10 larger than the post-Newtonian effects in the determination of the equipotential level surface!

Conflict of Interest Statement

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


We thank the referees for valuable comments and suggestions for improving the presentation of the manuscript. The work of SK has been supported by the grant No 14-27-00068 of the Russian Science Foundation (RSF).

Supplementary Material

The Supplementary Material for this article can be found online at:


1. ^Post-newtonian definitions of mass, center of mass, and other multipole moments can be found, for example, in Kopeikin et al. (2011).

2. ^We remind the reader that the coordinates on the surface of the PN-ellipsoid are denoted as x¯, ȳ, z¯.


Ablain, M., Cazenave, A., Larnicol, G., Balmaseda, M., Cipollini, P., Faugère, Y., et al. (2015). Improved sea level record over the satellite altimetry era (1993-2010) from the Climate Change Initiative project. Ocean Sci. 11, 67–82. doi: 10.5194/os-11-67-2015

CrossRef Full Text | Google Scholar

Ansorg, M., Fischer, T., Kleinwächter, A., Meinel, R., Petroff, D., and Schöbel, K. (2004). Equilibrium configurations of homogeneous fluids in general relativity. Mon. Not. Roy. Astron. Soc. 355, 682–688. doi: 10.1111/j.1365-2966.2004.08371.x

CrossRef Full Text | Google Scholar

Bardeen, J. M. (1971). A Reexamination of the Post-Newtonian Maclaurin Spheroids. Astrophys. J. 167, 425. doi: 10.1086/151040

CrossRef Full Text | Google Scholar

Beverini, N., Allegrini, M., Beghi, A., Belfi, J., Bouhadef, B., Calamai, M., et al. (2014). Measuring general relativity effects in a terrestrial lab by means of laser gyroscopes. Laser Phys. 24:074005. doi: 10.1088/1054-660X/24/7/074005

CrossRef Full Text | Google Scholar

Bjerhammar, A. (1985). On a relativistic geodesy. Bull. Géodésique 59, 207–220. doi: 10.1007/BF02520327

CrossRef Full Text | Google Scholar

Bondarenko, N. P., and Pyragas, K. A. (1974). On the Equilibrium figures of an ideal rotating liquid in the post-newtonian approximation of general relativity. II: Maclaurin's P-Ellipsoid. Astrophys. Space Sci. 27, 453–466. doi: 10.1007/BF00643890

CrossRef Full Text | Google Scholar

Bondarescu, R., Bondarescu, M., Hetényi, G., Boschi, L., Jetzer, P., and Balakrishna, J. (2012). Geophysical applicability of atomic clocks: direct continental geoid mapping. Geophys. J. Int. 191, 78–82. doi: 10.1111/j.1365-246X.2012.05636.x

CrossRef Full Text | Google Scholar

Bondarescu, R., Schärer, A., Lundgren, A., Hetényi, G., Houlié, N., Jetzer, P., et al. (2015). Ground-based optical atomic clocks as a tool to monitor vertical surface motion. Geophys. J. Int. 202, 1770–1774. doi: 10.1093/gji/ggv246

CrossRef Full Text | Google Scholar

Brumberg, V. A. (1991). Essential Relativistic Celestial Mechanics. Bristol: Adam Hilger.

Google Scholar

Chandrasekhar, S. (1965). The post-newtonian effects of general relativity on the equilibrium of uniformly rotating bodies. I. The maclaurin spheroids and the virial theorem. Astrophys. J. 142, 1513–1518. doi: 10.1086/148433

CrossRef Full Text | Google Scholar

Chandrasekhar, S. (1967). Ellipsoidal figures of equilibrium - a historical account. Commun. Pure Appl. Math. 20, 251–265. doi: 10.1002/cpa.3160200203

CrossRef Full Text | Google Scholar

Chandrasekhar, S. (1967a). The post-newtonian effects of general relativity on the equilibrium of uniformly rotating bodies. II. The deformed figures of the maclaurin spheroids. Astrophys. J. 147, 334–352. doi: 10.1086/149003

CrossRef Full Text | Google Scholar

Chandrasekhar, S. (1967b). The post-newtonian effects of general relativity on the equilibrium of uniformly rotating bodies. III. The deformed figures of the jacobi ellipsoids. Astrophys. J. 148, 621–644. doi: 10.1086/149183

CrossRef Full Text | Google Scholar

Chandrasekhar, S. (1967c). The post-newtonian effects of general relativity on the equilibrium of uniformly rotating bodies.IV. The roche model. Astrophys. J. 148, 645–649. doi: 10.1086/149184

CrossRef Full Text | Google Scholar

Chandrasekhar, S. (1969). Ellipsoidal Figures of Equilibrium. New Haven, CT: Yale University Press.

Google Scholar

Chandrasekhar, S. (1971a). The post-newtonian effects of general relativity on the equilibrium of uniformaly rotating bodies.VI. The deformed figures of the jacobi ellipsoids. Astrophys. J. 167, 455–463. doi: 10.1086/151042

CrossRef Full Text | Google Scholar

Chandrasekhar, S. (1971b). The post-newtonian effects of general relativity on the equilibrium of uniformly rotating bodies. V. The deformed figures of the maclaurin spheroids. Astrophys. J. 167, 447–453. doi: 10.1086/151041

CrossRef Full Text | Google Scholar

Chandrasekhar, S., and Elbert, D. D. (1974). The deformed figures of the Dedekind ellipsoids in the post-Newtonian approximation to general relativity. Astrophys. J. 192, 731–746. doi: 10.1086/153111

CrossRef Full Text | Google Scholar

Chandrasekhar, S., and Elbert, D. D. (1978). The deformed figures of the Dedekind ellipsoids in the post-Newtonian approximation to general relativity - Corrections and amplifications. Astrophys. J. 220, 303–313. doi: 10.1086/155906

CrossRef Full Text | Google Scholar

Chandrasekhar, S., and Miller, J. C. (1974). On slowly rotating homogeneous masses in general relativity. Mon. Not. Roy. Astron. Soc. 167, 63–80. doi: 10.1093/mnras/167.1.63

CrossRef Full Text | Google Scholar

Ciufolini, I., and Wheeler, J. A. (1995). Gravitation and Inertia. Princeton, NJ: Princeton University Press.

Google Scholar

Drewes, H. (ed.). (2009). Geodetic Reference Frames. Berlin: Springer. doi: 10.1007/978-3-642-00860-3

CrossRef Full Text | Google Scholar

Falke, S., Lemke, N., Grebing, C., Lipphardt, B., Weyers, S., Gerginov, V., et al. (2014). A strontium lattice clock with 3 × 10−17 inaccuracy and its frequency. New J. Phys. 16:073023. doi: 10.1088/1367-2630/16/7/073023

CrossRef Full Text | Google Scholar

Fock, V. A. (1964). The Theory of Space, Time and Gravitation, 2nd Edn. (Trans. N. Kemmer). New York, NY: Macmillan.

Google Scholar

Fu, L.-L., and Haines, B. J. (2013). The challenges in long-term altimetry calibration for addressing the problem of global sea level change. Adv. Space Res. 51, 1284–1300. doi: 10.1016/j.asr.2012.06.005

CrossRef Full Text | Google Scholar

Galtsov, D. V., Tsvetkov, V. P., and Tsirulev, A. N. (1984). The spectrum and polarization of the gravitational radiation of pulsars. Sov. Phys. JETP 59, 472–477.

Google Scholar

Gürlebeck, N., and Petroff, D. (2010). The axisymmetric case for the post-newtonian dedekind ellipsoids. Astrophys. J. 722, 1207–1215. doi: 10.1088/0004-637X/722/2/1207

CrossRef Full Text | Google Scholar

Gürlebeck, N., and Petroff, D. (2013). A generalized family of post-newtonian dedekind ellipsoids. Astrophys. J. 777, 1–16. doi: 10.1088/0004-637X/777/1/60

CrossRef Full Text | Google Scholar

Heiskanen, W. A., and Moritz, H. (1967). Physical Geodesy. San Francisco, CA: W. H. Freeman.

Google Scholar

Hobson, E. W. (1909). On the second mean-value theorem of the integral calculus. Proc. Lond. Math. Soc. s2-7, 14–23. doi: 10.1112/plms/s2-7.1.14

CrossRef Full Text | Google Scholar

Hofmann-Wellenhof, B., and Moritz, H. (2006). Physical Geodesy. Berlin: Springer.

Google Scholar

Hurst, R. B., Rabeendran, N., Wells, J.-P. R., and Schreiber, K. U. (2015). “Large ring laser gyroscopes: towards absolute rotation rate sensing,” in Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, volume 9444 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series (Bellingham, WA), 944407.

Google Scholar

Kómár, P., Kessler, E. M., Bishof, M., Jiang, L., Sórensen, A. S., Ye, J., et al. (2014). A quantum network of clocks. Nat. Phys. 10, 582–587. doi: 10.1038/nphys3000

CrossRef Full Text | Google Scholar

Kopeikin, S., Efroimsky, M., and Kaplan, G. (2011). Relativistic Celestial Mechanics of the Solar System. Berlin: Wiley. doi: 10.1002/9783527634569

CrossRef Full Text | Google Scholar

Kopeikin, S., Han, W., and Mazurova, E. (2016). Post-Newtonian reference-ellipsoid for relativistic geodesy. Phys. Rev. D. eprint arXiv:1510.03131

Google Scholar

Kopeikin, S. M., Mazurova, E. M., and Karpik, A. P. (2015). Towards an exact relativistic theory of Earth's geoid undulation. Phys. Lett. A 379, 1555–1562. doi: 10.1016/j.physleta.2015.02.046

CrossRef Full Text | Google Scholar

Kopeikin, S. M., and Petrov, A. N. (2013). Post-newtonian celestial dynamics in cosmology: field equations. Phys. Rev. D 87:044029. doi: 10.1103/PhysRevD.87.044029

CrossRef Full Text | Google Scholar

Kopeikin, S. M., and Petrov, A. N. (2014). Dynamic field theory and equations of motion in cosmology. Ann. Phys. 350, 379–440. doi: 10.1016/j.aop.2014.07.029

CrossRef Full Text | Google Scholar

Kopejkin, S. M. (1991). Relativistic manifestations of gravitational fields in gravimetry and geodesy. Manuscripta Geodaetica 16, 301–312.

Google Scholar

Landau, L. D., and Lifshitz, E. M. (1975). The Classical Theory of Fields. Oxford: Pergamon Press.

Google Scholar

Lightman, A. P., Press, W. H., Price, R. H., and Teukolsky, S. A. (1975). Problem Book in Relativity and Gravitation. Princeton, NJ: Princeton University Press.

Google Scholar

Mai, E. (2014). Time, atomic clocks, and relativistic geodesy. Report No 124, Deutsche Geod'́atische Kommission der Bayerischen Akademie der Wissenschaften (DGK), 128. Available online at:

Mai, E., and Müller, J. (2013). General remarks on the potential use of atomic clocks in relativistic geodesy. ZFV Zeitschrift Geodasie Geoinformation Landmanagement 138, 257–266.

Google Scholar

Meinel, R., Ansorg, M., Kleinwächter, A., Neugebauer, G., and Petroff, D. (2008). Relativistic Figures of Equilibrium. Cambridge: Cambridge University Press. doi: 10.1017/CBO9780511535154

CrossRef Full Text | Google Scholar

Müller, J., Soffel, M., and Klioner, S. A. (2008). Geodesy and Relativity. J. Geod. 82, 133–145. doi: 10.1007/s00190-007-0168-7

CrossRef Full Text | Google Scholar

Nicholson, T. L., Campbell, S. L., Hutson, R. B., Marti, G. E., Bloom, B. J., McNally, R. L., et al. (2015). Systematic evaluation of an atomic clock at 2 × 10−18 total uncertainty. Nat. Commun. 6, 1–8. doi: 10.1038/ncomms7896

PubMed Abstract | CrossRef Full Text | Google Scholar

Oltean, M., Epp, R. J., McGrath, P. L., and Mann, R. B. (2015). Geoids in general relativity: geoid quasilocal frames. eprint arXiv:1510.02858

Google Scholar

Petit, G., and Luzum, B. (2010). IERS conventions. IERS Tech. Note 36, 179.

Google Scholar

Petit, G., Wolf, P., and Delva, P. (2014). “Atomic time, clocks, and clock comparisons in relativistic spacetime: a review,” in Frontiers in Relativistic Celestial Mechanics. Vol. 2 Applications and Experiments, ed S. M. Kopeikin (Berlin: W. de Gruyter), 249–279.

Google Scholar

Petroff, D. (2003). Post-Newtonian Maclaurin spheroids to arbitrary order. Phys. Rev. D 68, 104029. doi: 10.1103/physrevd.68.104029

CrossRef Full Text | Google Scholar

Pizzetti, P. (1913). Principii Della Teoria Meccanica Della Figura dei Pianeti. Pisa: E. Spoerri.

Google Scholar

Plag, H.-P., and Pearlman, M. (ed.). (2009). Global Geodetic Observing System. Dordrecht: Springer. doi: 10.1007/978-3-642-02687-4

CrossRef Full Text | Google Scholar

Poli, N., Oates, C. W., Gill, P., and Tino, G. M. (2013). Optical atomic clocks. Nuovo Cimento Riv. Ser. 36, 555–624. doi: 10.1393/ncr/i2013-10095-x

CrossRef Full Text | Google Scholar

Pyragas, K. A., Bondarenko, N. P., and Kravtsov, O. V. (1974). On the equilibrium figures of an ideal rotating liquid in the post-newtonian approximation of general relativity. I: equilibrium conditions. Astrophys. Space Sci. 27, 437–452. doi: 10.1007/BF00643889

CrossRef Full Text | Google Scholar

Pyragas, K. A., Bondarenko, N. P., and Kryshtal, A. N. (1975). On the equilibrium figures of an ideal rotating fluid in the post-Newtonian approximation of general relativity. III - Stability of the forms of equilibrium. Astrophys. Space Sci. 33, 75–97. doi: 10.1007/BF00646009

CrossRef Full Text | Google Scholar

Rendall, A. D. (1990). Convergent and divergent perturbation series and the post-minkowskian approximation scheme. Class. Quant. Grav. 7, 803–812. doi: 10.1088/0264-9381/7/5/010

CrossRef Full Text | Google Scholar

Schreiber, K. U. (2013). Variations of earth rotation from ring laser gyroscopes: one hundred years of rotation sensing with optical interferometry (invited). AGU Fall Meet. Abstr. Available online at:

Soffel, M., Klioner, S. A., Petit, G., Wolf, P., Kopeikin, S. M., Bretagnon, P., et al. (2003). The IAU 2000 resolutions for astrometry, celestial mechanics, and metrology in the relativistic framework: explanatory supplement. Astron. J. (USA) 126, 2687–2706. doi: 10.1086/378162

CrossRef Full Text | Google Scholar

Soffel, M. H. (1989). Relativity in Astrometry, Celestial Mechanics and Geodesy. Berlin: Springer-Verlag. doi: 10.1007/978-3-642-73406-9

CrossRef Full Text | Google Scholar

Torge, W., and Müller, J. (2012). Geodesy, 4th Edn. Berlin: De Gruyter. doi: 10.1515/9783110250008

CrossRef Full Text | Google Scholar

Tsirulev, A. N., and Tsvetkov, V. P. (1982a). Rotating post-newtonian near ellipsoidal configurations of a magnetized homogeneous fluid - Part I. Sov. Astron. 26, 289–292.

Google Scholar

Tsirulev, A. N., and Tsvetkov, V. P. (1982b). Rotating post-newtonian near ellipsoidal configurations of a magnetized homogeneous fluid - Part II. Sov. Astron. 26, 407–412.

Google Scholar

Tsvetkov, V. P., and Tsirulev, A. N. (1983). Gravitational waves emitted by a spinning magnetized blob of homogeneous post-newtonian gravitating fluid. Sov. Astron. 27, 66–69.

Google Scholar

Vaníček, P., and Krakiwsky, E. J. (1986). Geodesy, the Concepts, 2nd Edn. Amsterdam: North Holland.

Google Scholar

Weinberg, S. (1972). Gravitation and Cosmology: Principles and Applications of the General Theory of Relativity. New York, NY: John Wiley & Sons, Inc.

Google Scholar

Will, C. M. (1993). Theory and Experiment in Gravitational Physics. Cambridge: Cambridge University Press. doi: 10.1017/CBO9780511564246

CrossRef Full Text | Google Scholar

Zelmanov, A., Rabounski, D., Crothers, S. J., and Borissova, L. (2006). Chronometric Invariants: on Deformations and the Curvature of Accompanying Space. Rehoboth, NM: American Research Press. Available online at:

Keywords: physical geodesy, reference ellipsoid, geoid, relativistic geodesy, reference frames, general relativity and gravitation

Citation: Kopeikin SM (2016) Reference Ellipsoid and Geoid in Chronometric Geodesy. Front. Astron. Space Sci. 3:5. doi: 10.3389/fspas.2016.00005

Received: 22 November 2015; Accepted: 08 February 2016;
Published: 25 February 2016.

Edited by:

Agnes Fienga, Observatoire de la Côte d'Azur, France

Reviewed by:

Alberto Vecchiato, Osservatorio Astrofisico di Torino - INAF, Italy
Ramesh Jaga Govind, University of Cape Town, South Africa

Copyright © 2016 Kopeikin. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Sergei M. Kopeikin,