This computation is being done with the method used in "Gravitation and Relativity" by M. G. Bowler (1976, Pergamon Press, Elmsford NY) pp. 148 - 150.
The setup begins with the a diffeomophism of the Distorted Schwarzschild Solution derived in Ïmposing a 4-Dimensional background on General Relativity" by Edward M. Schaefer in "Theoretical Physics: MRST2002: A Tribute to George Leibbrabdt / AIP Conference Proceedings Volume 646" (2002, Amreican Institute of Physics, Melville NY) on p. 156 as Eq. 21.
ds^2 = A(rho)dt^2 - [1/A(rho)][d rho^2 + rho^2 (d theta^2 + (sin^2 theta) d phi^2)],
A(rho) = rho^2/[m + sqrt(rho^2 + m^2)]^2,
where rho is a radial coordinate as measured from a given worldline in the spacetime, and m is the geometrized mass of the SSNRMO as measured from the same worldline. The other symbols are as used in the Schwarzschild Solution. The coordinate transformation between the radial coordinates of the distorted Schwarzschild Solution and the above diffeomorphism is rho = sqrt(r[r-2m]).
Since this metric contains no off-diagonal terms, the technique of Bowler's Eqs. (9.10.3) through (9.10.14) can be used to provide constants of motion for the energy and angular momentum of the test object. This includes the simplification technique of using the plane theta = Pi/2, in which case d theta = 0 and sin(theta) = 1. The expressions are:
d tau = K (A(rho)) dt = K (rho^2/[m + sqrt(rho^2 + m^2)]^2) dt and
[m + sqrt (rho^2 + m^2)]^2 d tau/d phi = H
where H and K are the desired constants.
The computations start by defining B = sqrt (1/A(rho)), which appears regularly in this exercise.
> B := (m + sqrt (rho(t)^2 + m^2))/rho(t);
|
This form of the metric equation is that of Bowler's Eq. (9.10.16). It arises by doing substitutions of d tau and d phi in the PGR SSNRMO metric (in the plane theta = Pi/2) using the equations for the constants of motion given above, and then dividing the equation by A(rho)^2 to leave just K^2 on the lhs. h = H*K, and is therefore itself a constant in these calculations.
> metr := K^2 = B^2 - B^6 * diff(rho(t),t)^2 - h^2/(rho(t)^2 * B^2);
|
This is a description of h in the above equation. diff_ph_t = d phi/d t, and will be given an appropriate value by substitution when hh is used.
> hh := B^4 * rho(t)^2 * diff_ph_t;
|
We want to know how the orbiting particle's position varies over time. So the metric is differentiated wrt time to give:
> diff (metr, t);
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Differentail-cubed terms dropped (since the are insignificant for a body like Mercury)
> d_metr := % - (0 = -6*(m+sqrt(rho(t)^2+m^2))^5 * diff(rho(t),t)^3 / (rho(t)^5*sqrt(rho(t)^2+m^2))+6*(m+sqrt(rho(t)^2+m^2))^6*diff(rho(t),t )^3/(rho(t)^7));
| |||||||||||||||||||||||||||||||||||||||||||||||||||||
This equation is now reorganized and normalized to leave (d/d t)^2 rho on the lhs.
> % - (-2*(m+sqrt(rho(t)^2+m^2))^6*diff(rho(t),t)*diff(rho(t),``(t,2))/(rho(t)6) = -2*(m+sqrt(rho(t)2+m2))6*diff(rho(t),t)*diff(rho(t),``(t,2))/(rho( t)^6));
| |||||||||||||||||||||||||||||||||||||||||||||||||||||
> rho_ddiff := expand(%/(B^6 * 2 * diff (rho(t),t)));
| |||||||||||||||||||||||||||||||||||||||||||||||||||
To obtain the period of revolution, a circular orbit is assumed, or (d/dt)^2 rho = 0. h is now replaced with its value, which is where the (d phi/dt)^2 term comes from.
> 0 = subs(h=subs (diff_ph_t = diff (phi(t),t), hh),rhs(rho_ddiff));
| |||||||||||||||||||||||||||||||||||||||||||
Now reorganize so that d phi/dt is on the lhs.
> %-(rho(t)^3*diff(phi(t),t)^2/((m+sqrt(rho(t)^2+m^2))*sqrt(rho(t)^2+m^2)) = rho(t)^3*diff(phi(t),t)^2/((m+sqrt(rho(t)^2+m^2))*sqrt(rho(t)^2+m^2))) ;
| |||||||||||||||||||||||||||||||||||||||||||
> normal (%);
| |||||||||||||||||||||||||||||||
> sq_df_phi := %/(-rho(t)^3/((m+sqrt(rho(t)^2+m^2))*sqrt(rho(t)^2+m^2)));
|
Convert the above into a expression in terms of r. This makes for a more appropriate comparison with the results of the Schwarzschild Solution since the map between rho and r is rho = sqrt(r(r-2m)).
> subs (rho(t) = sqrt (r*(r-2*m)), rhs(sq_df_phi));
|
Cast the recipricol of the rhs of the above as a series. This is more computationally stable and easier to analyze. dt^2 = this * d phi^2. It is notable that this is through the 5th order excatly the same as for Schwarzschild.
> series (1/%, r=infinity);
|
The revolution time is obtained by integration over the range of phi (0..2*Pi) on the square root of the above.
> rev_time := 2 * Pi * sqrt(%):
For the oscillation time (the time it takes to return to the same radial coordinate), one can replace rho by R + sigma, to obtain a result for small quantities of the from (d/dt)^2 sigma = -sigma*(d X(R)/dR), where X(R) is the rhs of rho_ddiff. This produces a solution of sigma = C*sqrt(dX(R)/dR)*sin t, which gives an oscillation period of 2*Pi/sqrt(dX(R)/dR)
Continuing to use rho, sqrt(d X(rho)/d rho)) will now be computed.
> oscl_eqn := diff (subs (rho(t) = rho, rhs (rho_ddiff)), rho);
| |||||||||||||||||||||||||||||||||||||||||||||||||||||||
Replace h in the above by its value. We now use the result obtained above for d phi / dt in the place of diff_ph_t.
> hh2:= subs(rho(t)=rho,subs(diff_ph_t=sqrt(rhs(sq_df_phi)),hh));
|
> subs (h=hh2, oscl_eqn);
| |||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Obtain the series for of the reciporicol (in terms of r), and multiply by 2*Pi to get the oscillation time.
> series (rationalize(-1/subs (rho=sqrt(r*(r-2*m)), %)), r=infinity);
| |||||||||||||||||||||||||
> oscl_time := 2*Pi*sqrt(%):
Now the perihelion advance per radian can be computed. This is the difference between the oscillation time divided by the revolution time and 1. For consistency, this will be determined as a function of r instead of rho.
As a series, the first-order term is 3*m/rho, and so an advance of 6*Pi*m/rho per revolution is obtained, which is the same as for the Schwarzschild Solution. The second order term for the Schwarzschild Solution is (27/2)*m^2/r^2, meaning that the metric used herein causes a difference of 0.5*pi*m^2/rho^2 in the predicted perihelion shift, with PGR prediction the greater shift. (The Schwarzschild shift was computed on a different worksheet: "Schwarzschild perihelion shift.mws".)
> series (oscl_time/rev_time - 1, r=infinity);
|
> perih_adv := 2*Pi*convert(%,polynom):
The computation of the perihelion advance for Mercury begins here.
astronomical unit in cm.
> au := 1.496 * 10^13;
|
the semi-major axis of Mercury
> a := 0.3871 * au;
|
the eccentricity of Mercury's orbit.
> e := .2056;
|
Mercury's semi-latus rectum, the radius of a circular orbit with the same average potential (as a function of time) as Mercury's.
> ll := a * (1 - e^2);
|
The mass of the Sun in grams.
> M_s := 1.989 * 10^33;
|
The speed of light
> c := 29979845800;
|
The gravitational constant
> G := 6.6726 * 10^(-8);
|
The geometrized mass of the Sun. (m_s has units of length).
> m_s := M_s * G /c^2;
|
The PGR contribution to the perihelion advance of Mercury in radians/orbit.
> shift_rate := evalf (subs (r=ll, m=m_s, perih_adv));
|
The length of Mercury's revolution in years.
> merc_rev := 87.97/365.24219;
|
The GR perihelion advance of Mercury in arc-seconds per century.
> evalf((100/%)*shift_rate/Pi*180*3600);
|
The difference between the predictions for the Schwarzschild metric and my own. Note that the difference of ~1*10^(-6) arc-seconds per century is currently way too small to detect.
in radians per revolution,
> evalf(-0.5*Pi*(m_s/ll)^2);
|
in arc-seconds per century,
> evalf ((100/merc_rev)*%/Pi*180*3600);
|