The Gravity Train Revisited

Introduction 

The gravity train is a conceptual of design of transportation making use of the gravity of the planet, first proposed by 17th century scientist Robert Hooke in a letter to Issac Newton. It features a tunnel through the Earth.

Suppose a straight tunnel was drill through the Earth. A train operating solely on gravity starts from one end of the tunnel. Assuming there is no friction inside the tunnel, and the density of the Earth is uniform. The problem asks how long it takes for the train to arrive at the other end.

In the classic solution, the spin (rotation) of the Earth is not considered. (Strictly speaking, the rotation effect can be factorized in by using a locally measured gravitational acceleration $g$). The answer then is $t = \pi \sqrt{\frac{R_E}{g}} \simeq 42 \mathrm{min.}$, regardless of the exact location of the tunnel, where $R_E \simeq 6371 \mathrm{km}$ is the radius of the Earth, $g \simeq 9.8 \mathrm{m/s^2}$ is the gravitational acceleration. The rotation of the planet introduces a second time scale: $T_s \simeq 24 \mathrm{hr} \simeq 17 T_G$, where  $T_G = 2\pi \sqrt{\frac{R_E}{g}}$. The separation of the scales validates the neglecting of the Earth spin in the classical solution. Nevertheless, in this essay we want to look at the correction effect introduced by the spin of the Earth.

Fig. 1

General Settings

Let the center of the Earth be the origin of our coordinate system. The tunnel is characterized by the position of the center of the tunnel $\vec d$, and the orientation of the tunnel $\hat \tau$, $\vec d \cdot \hat \tau = 0$. Let $\vec{r} = \vec{d} + x \hat{\tau}$ be the position of the train, where $x$ is the distance of the train from the center of the tunnel, $\vec{v} = \dot x\hat\tau$ the velocity of the train. The spin of the Earth is described by an angular velocity vector $\vec\Omega$. Suppose $\Omega \cos \theta \equiv \vec \Omega \cdot \hat \tau$, $\Omega d \cos\phi = \vec\Omega \cdot \vec d$.

In the reference frame of the Earth (a non-inertial frame), there are four forces acting on the train: the gravity $\vec F_G = -\frac{GM_Em}{R^3_E} \vec r \equiv -mg \frac{\vec r}{R_E}$, the centrifugal force $\vec F_e = - m \vec\Omega \times (\vec\Omega \times \vec r)$, the Coriolis force $\vec F_c = 2 m \vec v \times \vec \Omega$, and the contact force. Since there is no friction, the contact force simply counter all forces in the direction that perpendicular to the tunnel. We only need to consider the force along the tunnel. $\hat \tau \cdot \vec F_G = -m g \frac{x}{R_E}$, $\hat \tau \cdot \vec F_e = m \Omega^2 x - m \Omega^2 (d \cos\phi + x \cos \theta)$, $\hat \tau \cdot \vec F_c = 0$. The total force along the tunnel is $F = -mg\frac{x}{R_E} + m \Omega^2 ( (1-\cos\theta)x - d\cos\phi)$.

The equation of motion thus reads, \[
\ddot x = \left( -\frac{g}{R_E} + \Omega^2 (1-\cos\theta) \right) x - \Omega^2 d \cos\phi
\] The solution is a harmonic oscillator with angular frequency $\omega = \frac{g}{R_E} - \Omega^2 (1-\cos\theta)$ or period $T = 2\pi \sqrt{\frac{R_E}{g - \Omega^2R_E(1-\cos\theta)}}$. The center of the harmonic oscillator is not in the middle of the tunnel. In stead, it is shifted away from there by $\Delta = \frac{\Omega^2dR_E\cos\phi}{g - \Omega^2 R_E(1-\cos\theta)}$.

To calculate the travel time, we need to consider two cases. If the train starts from the high altitude side, it arrives at the low altitude terminal in time $t = \pi \sqrt{\frac{R_E}{g - \Omega^2R_E(1-\cos\theta)}}\left(1 - \frac{1}{\pi}\arccos\frac{\sqrt{R_E^2-d^2}-\Delta}{\sqrt{R_E^2-d^2}+\Delta} \right)$. If the train starts from low altitude terminal (farther from the axis of the Earth) with zero velocity, it would never reach the high altitude side. Some minimum velocity is needed in order to arrive the other side, $v_{\min} = 2 \Omega \sqrt{\left( g - \Omega^2R_E (1-\cos\theta) \right) d\cos\phi} \sqrt[4]{1-\frac{d^2}{R^2_E}}$. by doing so, it also takes time $t$ to arrive at the destination.

Discussions

The spin of the Earth affects the problem in two ways. The correction in the period of the oscillator is about $\sim \frac{\Omega^2 R_E}{g}\sin^2\frac{\theta}{2} = 0.345\% \sin^2\frac{\theta}{2}$. The second correction is through the displacement of the center of the harmonic oscillator. We compare $\Delta$ and $l = \sqrt{R_E^2 - d^2}$ half the length of the tunnel. $\Delta \simeq \frac{\Omega^2 R_E}{g} \cos\phi d = 0.345371\% \cos\phi \sqrt{R_E^2 - l^2}$. Compare $\Delta_{\max}$ (with $\phi = 0$) and $l$ (see Fig. 2). The correction from $\Delta$ is only small for tunnels either of thousands of kilometers long or well oriented (i.e. choosing a good $\phi$). Thousands of kilometers long tunnel would be for transcontinental express.

Fig. 2, the numbers for the tunnel length $l$ and maximum shift of the harmonic oscillator center $\Delta_{\max}$ are in kilometers.

We have only considered the straight tunnels. It turns out, straight tunnels are not the fastest ones for gravity train. The fastest tunnels are the brachistochrone tunnels. The correction due to the Earth rotation can be considered accordingly.

On dynamics of the Kepler problem

The two-body central-force problem is also known as the Kepler problem is a classical problem in the domain of classical mechanics. The textbooks pay much attention on the orbits, at least partly because
the orbits have closed-form. In this post, I put emphasis on the dynamic solutions, i.e. the time evolution of the motion, which are the direct solutions of the dynamic equation.

The notations

$r$: the radial coordinate, the distance to the center-of-mass;
$\theta$: the angular coordinate;

$m$: the reduced mass;
$l$: the module of the angular momentum ($l>0$);
$V(r) = - \frac{k m}{r}$: the potential energy. ($k > 0$);
$E$: total energy. ($ E < 0$);
$\varepsilon$: the eccentricity. ($-1 < \varepsilon < 1$);
$a$: the semi-major axis;
$r_{\min} = a(1-\varepsilon), r_{\max} = a(1+\varepsilon)$: the periapsis and apoapsis distances.
$T=2\pi\sqrt{\frac{a^3}{k}}$: the period;

Here we have listed two sets of parameters: 1. the physical parameters $m, l, E, k$ (more precisely, $\frac{l}{m}, \frac{E}{m}, k$); 2. the orbit parameter (kinematical parameters) $a, \varepsilon, T$. Each set can determine the dynamic solution. We shall see their relation later.

The equations of motion

[ equation 1 ]\[ \ddot r - \frac{l^2}{m^2 r^3} + \frac{1}{m}\frac{\partial }{\partial r} V(r) = 0; \\
\dot \theta = \frac{l}{m r^2}; \\
r(t=0) = r_{\min}, \quad \dot r(t=0) = 0, \quad \theta(t=0) = 0;
\] where $\dot{}$ stands for the time derivative.

The orbits

The orbit is an equation of the coordinates $F(r,\theta) =0$. It is convenient to express $r$ as a function of $\theta$: $r = r(\theta)$. The classical way to get the orbit is first to use the chain rule: $\frac{d}{dt} = \dot{\theta}\frac{d}{d\theta} = \frac{l}{mr^2}\frac{d}{d\theta}$. Then the equation of motion becomes
[ equation 2 ] \[
r''_\theta - 2\frac{r^{'2}_\theta}{r} - r + \frac{m}{l^2}r^4\frac{\partial}{\partial r}V(r) = 0.
\] By doing a change of variable $u \equiv \frac{1}{r}$, the orbit equation becomes,
[ equation 3 ] \[
u''_\theta + u + \frac{m}{l^2}\frac{\partial}{\partial u}V(1/u) = 0
\]
For $V = -\frac{k m}{r}$, it becomes $ u''_\theta + u = \frac{km^2}{l^2}$. And the solution is,
[ solution 1 ] \[ u = \frac{1}{r} = \frac{km^2}{l^2}(1+\varepsilon \cos(\theta)) \] The orbit is an ellipse. The semi-major axis $a (1-\varepsilon^2) = \frac{l^2}{km^2}$.

The dynamic solutions

To solve for the dynamics $\vec{r}(t)$ (in polar coordinate $\vec{r}(t) = (r(t), \theta(t))$), we shall return to equation 1. Note that equation 1 does not depend on $t$ explicitly. This type of 2nd order differential equation can always be solved by using the chain rule: $\frac{d}{dt} = \frac{dr}{dt}\frac{d}{dr}$. Equation 1 becomes:
[ equation 4 ] \[ \dot r \frac{d \dot r}{dr} + \frac{d}{dr} \left( \frac{1}{m}V(r) + \frac{l^2}{2m^2 r^2} \right) = 0 \\ \implies \frac{1}{2} m \dot r^2 + V(r) + \frac{l^2}{2mr^2} = E \] Substitute the orbit equation, we get two relations $a(1-\varepsilon^2) = \frac{l^2}{km^2}, E = -\frac{km}{2a}$ thus \[ a = -\frac{km}{2E}, \quad \varepsilon = \sqrt{1 + 2\frac{l^2E}{k^2m^3}}, \quad T = 2\pi \sqrt{-\frac{k^2m^3}{8E^3}}, \] where the expression about $T$ comes from Kepler's third law. We'll derive it later. The second one in equation 4 is nothing but the energy conservation. To solve equation 4, we note $\dot r = \frac{dr}{dt}$. Thus, \[ t = \int_C \frac{dr}{\pm\sqrt{\frac{2}{m} \left(E- \frac{l^2}{2m r^2}-V(r) \right)}} \] where $C$ is an integral path. The path $C$ also encodes the initial condition. In our case, $r(t=0) = a(1-\varepsilon)$. For $V = -\frac{km}{r}$, $E - \frac{l^2}{2m r^2} - V(r)$ in the denominator of equation 4 is a rational polynomial and can be rewritten as $E (r-r_+)(r-r_-)/r^2$, where $r_\pm$ are the roots of the polynomial. On the other hand, $E - \frac{l^2}{2m r^2_\pm} - V(r_\pm) = 0 \Rightarrow \dot r = 0$. So $r_\pm$ are the two extrema, $r_\pm = a(1\pm \varepsilon)$. Equation 5 can be written as
[ equation 5 ] \[ t = \sqrt{-\frac{m}{2E}} \int_C \frac{r dr}{\pm\sqrt{(r- r_-)(r_+-r)}} \] Note that $\sqrt{-\frac{m}{2E}} = \sqrt{\frac{a}{k}}$.

Fig. The plot of the imaginary part of the integrand with $r_- = 1, r_+ = 2$. $r_\pm$ are the two brach points of the integrand.

To evaluate this integral, we need proper parametrization $r = r(\xi)$. There are several viable parametrization schemes.

scheme 1. If we integrate from $r_{\min}$ to $r< r_{\max}$, then we can evaluate the integral directly (numerically for example):
[ equation 6-1 ] \[
t = \sqrt{\frac{a}{k}} \int_{a(1-\varepsilon)}^r \frac{\rho d\rho}{\sqrt{(\rho-a+a\varepsilon)(a+a\varepsilon-\rho)}} \]

scheme 2. $r = a(1-\varepsilon)\cos^2\varphi + a(1+\varepsilon)\sin^2\varphi = a(1-\varepsilon\cos2\varphi), \quad (0 \le \varphi \le \frac{\pi}{2})$. So the integral becomes,
[ equation 6-2 ] \[
t = \sqrt{\frac{a}{k}}\int_0^{\frac{1}{2}\arccos\frac{a-r}{a\varepsilon}} 2a(1-\varepsilon\cos2\varphi) d \varphi \\
= \sqrt{\frac{a^3}{k}} \left[ \arccos\frac{a-r}{a\varepsilon}-\sqrt{\varepsilon^2-\left(\frac{a-r}{a}\right)^2}\right],
\] were $\arccos$ is defined on $[0, \pi]$. $t = \frac{1}{2}T$ when $r = a(1+\varepsilon)$. Then $T = 2\pi \sqrt{\frac{a^3}{k}}$, as we know from Kepler's third law. In order to get the expression $r = r(t)$, we need to invert the function, especially $f_\varepsilon(z) \equiv \arccos(z) - \varepsilon\sqrt{1-z^2}$. $f_\varepsilon(z)$ is monotone and its inverse function exists. Then $r(t) = f^{-1}_\varepsilon\left(2\pi\frac{t}{T}\right)$. Unfortunately, $f^{-1}_\varepsilon$ can not be expressed by elementary function. In numerical calculation, however, this is not as bad as one may think. The inverse function can be computed via Newton's algorithm.

scheme 3. We can also parametrize $r$ with $\theta$ directly, $r = a(1-\varepsilon^2)/(1+\varepsilon\cos\theta)$.
[ equation 6-3 ] \[t = \sqrt{\frac{a^3(1-\varepsilon^2)^3}{k}} \int_0^\theta \frac{\mathrm d \phi}{(1+\varepsilon\cos\phi)^2}.\] We can do a change of variable $s = \tan\frac{\phi}{2}$. Then $d\phi = \frac{2}{1+s^2}ds, \cos\phi = \frac{1-s^2}{1+s^2}$ The integral part becomes, \[
\int  \frac{\mathrm d \phi}{(1+\varepsilon\cos\phi)^2} = \int ds \frac{2(1+s^2)}{\left(1+s^2+\varepsilon(1-s^2)\right)^2} = \frac{1}{(1-\varepsilon)^2}\int \frac{ds}{s^2+\frac{1+\varepsilon}{1-\varepsilon}} - \frac{2\varepsilon}{(1-\varepsilon)^3}\int \frac{ds}{\left(s^2+\frac{1+\varepsilon}{1-\varepsilon}\right)^2}
\] Note that \[
\int \frac{ds}{s^2+a^2} = \frac{1}{a}\arctan \frac{s}{a} \\
\int \frac{ds}{(s^2+a^2)^2} = \frac{s}{2a^2(s^2+a^2)} + \frac{1}{2a^3}\arctan\frac{s}{a}
\]
Therefore,
\[
t = \sqrt{\frac{a^3}{k}} \left( 2 \arctan \left(\sqrt{\frac{1-\varepsilon}{1+\varepsilon}} \tan\frac{\theta}{2}\right)
- \sqrt{1-\varepsilon^2} \frac{\varepsilon\sin\theta}{1+\varepsilon\cos\theta} \right)
\] (note that the $\arctan$ function returns values on $[0, 2\pi]$. Sometimes, the so-called two-argument inverse tangent function $\arctan(x,y)$ is used.)

Again, in this case, we need to invert the function and obtain an expression for $\theta$.

With either of the above method, we can solve for $\vec{r}(t) = (r(t), \theta(t))$. Finally, we can put the orbit and dynamic solution together, in an animation.

Applications

1. the quasi-satellite

A quasi-satellite of a planet is an object that orbiting the same center and with the same period as the planet (co-orbital configuration). Usually, the quasi-satellite has a different eccentricity. The relative motion of the object with respect to the planet form a closed orbit. So, the object looks as if it is orbiting around the planet with irregular orbits (quasi-orbit). Historically, quasi-satellites resulted several false claims of the existence of other moons of the earth.

Fig. 1a, the relative orbit (yellow), the orbit of the planet (blue), the orbit of the object (purple).

Fig. 2, Blue: the relative orbit (yellow), the orbit of the planet (blue), the orbit of the object (purple).

cf. The quantum Hamilton-Jacobi theory of a central-force problem