<p>Write the mean anomaly as a Taylor series ($n\equiv 2\pi/P$):
$$\begin{eqnarray*}
M(t) \equiv \int_0^tn\,\mathrm{d}t &amp;=&amp; M_0 + \dot{M}_0t + \frac{1}{2}\ddot{M}_0t^2 + \mathcal{O}(t^3)\\
&amp;=&amp;n_0t + \frac{1}{2}\dot{n}_0t^2 + \mathcal{O}(t^3)\\
&amp;=&amp;\frac{2\pi}{P_0}t - \pi\frac{\dot{P}_0}{P_0^2}t^2+\mathcal{O}(t^3)\text{.}
\end{eqnarray*}$$</p>

<p>The periastron time (epoch) happens when $M(t) = 2\pi N$, where $N$ counts number of orbits, while $NP_0$ <em>would</em> be the periastron time if the orbital period stayed constant. Therefore, treating $N$ as continuous, the cumulative difference is just
$$t - \frac{M(t)}{2\pi}P_0 \approx \frac{1}{2}\frac{\dot{P}_0}{P_0}t^2\text{.}$$
If we want more accuracy, we can calculate the next term in the Taylor series,
$$\frac{1}{6}\ddot{n}_0t^3 = \frac{\pi}{3P_0^2}\left[\frac{2\dot{P}_0^2}{P_0}-\ddot{P}_0\right]t^3\text{,}$$
and if $\dot{P} = KP^{-5/3}$ for $\dot{K}\approx 0$, we can eliminate the $\ddot{P}_0$ term to give a cumulative shift of
$$t - \frac{M(t)}{2\pi}P_0 \approx \frac{1}{2}\frac{\dot{P}_0}{P_0}t^2 - \frac{11}{18}\left(\frac{\dot{P}_0}{P_0}\right)^2t^3\text{.}$$
According to the Weisberg &amp; Taylor paper, $|\dot{P}_0/P_0| \sim 10^{-16}\,\mathrm{s}^{-1}$, so the cubic term is negligible over the thirty-year time scale, and we can safely ignore anything past the quadratic term.</p>

<p>Just for kicks, we can solve $\dot{P} = KP^{-5/3}$ for constant $K$, which yields a cumulative periastron time shift of
$$\frac{3}{5x}\left[1+\frac{5}{3}xt-\left(1-\frac{8}{3}xt\right)^{5/8}\right]\text{.}$$
$K$ is not actually constant, being a function of the eccentricity $e$, although the rate of change of $e$ is tiny, on the order of $10^{-11}\,\mathrm{yr}^{-1}$, and we can ignore it on the scales we're talking about. Additionally, the proportionality formula is based on a orbital average of radiation reaction at 2.5-th post-Newtonian order.</p>

<hr>

<p>The values presented in the paper are for the epoch $T_0 = 52144.90097844$ (MJD):
$$\small{\begin{eqnarray*}
P_0 &amp;=&amp; 0.322997448930\,(4)\,\mathrm{d}\\
%\dot{P}_0 &amp;=&amp; −(2.4184\pm0.0009)\times 10^{-12}\\
%\dot{P}_0^{\scriptsize\text{corr}} &amp;=&amp; -(2.4056\pm 0.0051)\times 10^{-12}\\
\dot{P}_0^{\scriptsize\text{GTR}} &amp;=&amp; −2.40242\,(2)\times 10^{-12}\\
\end{eqnarray*}}$$
Thus, we can calculate the ratio $x_0 = \dot{P}_0/P_0 = -8.60867\,(7)\times 10^{-17}\,\mathrm{s}^{-1}$. The difference between $T_0$ and the beginning of 1975 is $\delta t = 9731\,\mathrm{d}$, which is not long enough for this ratio to change appreciably, since relative first-order correction would be around $|\dot{x}_0(\delta t)/x_0| \sim |x_0(\delta t)| \sim 10^{-7}$.</p>

<p>Alternatively, following <strong>zibadawa timmy</strong>, we can also solve $K = \dot{P}P^{5/3}$ explicitly for constant $K$, which is not exactly true but approximates well:
$$\newcommand{\constKeq}{\overset{\scriptsize\dot{K}\approx 0}{=}}\begin{eqnarray*}
P &amp;\constKeq&amp; \left[P_0^{8/3}+\frac{8}{3}\dot{P}_0P_0^{5/3}t\right]^{3/8}\!\!\!\!\!\text{,}\\
x &amp;\constKeq&amp; \frac{\dot{P}}{P} = \frac{x_0}{1+\frac{8}{3}x_0t}\text{,}
\end{eqnarray*}$$
and therefore $\left.x\right|_{1975}$ is indeed indistinguishable from $x_0$ at this level of accuracy.</p>

<p>I've graphed (in light blue) the result from the quadratic shift formula with this $x$  and superimposed the graph from the Weisberg &amp; Taylor paper (cropped on outer borders and then rescaled to be $500\times 500\,\mathrm{px}$):</p>

<p><img src="https://i.sstatic.net/jteR1.png" alt="Hulse–Taylor graph reconstruction"></p>

<p>It's an exact match except for a constant time shift, which is just a matter of when we set the $t = 0$ time. PSR B1913+16 was actually discovered in 1974, and the first data point is at the end of 1974. That is why my graph is shifted slightly by a fraction of a year.</p>

<hr>

<blockquote>
  <p>I am unclear whether your formula excludes the effect of periastron advance (4.226595 degrees per year):- (A) If not then the periastron advance effect is so large that the change in period due to gravitational wave emission is not visible on the graph and the agreement with theory cannot be assessed ...</p>
</blockquote>

<p>The GTR prediction of the cumulative periastron time shift is made through through equation (1) in the paper, which prescribes $\dot{P}_0 = \left.\dot{P}_b\right|_{t=0}$. The periastron time is well-defined because it depends on the closest approach, so the magnitude of the positional periastron advance $\langle\dot{\omega}\rangle = 4.226595^\circ/\mathrm{yr}$ is not directly relevant. GTR <em>also</em> predicts $\langle\dot{\omega}\rangle$, but that's not what we're considering here.</p>

<p>The mean anomaly is not actually an angle, despite counting off each periastron as $2\pi$. It isn't an angle even in exact Keplerian elliptic orbits, even though they have the exact same periastron at each period. Thus, $\langle\dot{\omega}\rangle$ is a different issue.</p>

<blockquote>
  <p>(B) If so then must we conclude that gravitational wave emission causes change in period of approx 4s in 10 years? </p>
</blockquote>

<p>The period changes by a bit less than $1\,\mathrm{ms}$ over a decade, as you have already calculated above.</p>

CC BY-SA 3.0
2014-11-15T13:56:27.660Z
That gives a curve that looks close to the original. I haven't digested your method yet. I am unclear whether your method excludes the effect of periastron advance (4.226595 degrees per year):- (A) If not then it could be that the periastron advance effect (not on P but on the calculation, e.g. ambiguity of Mean Anomaly) is so large that any change in P due to gravitational wave emission (GWE) is not visible on the graph and so agreement with GWE theory cannot be assessed. (B) If so then must we conclude that GWE alone causes big changes in period e.g. approx 40s in 30 years?
steveOw
CC BY-SA 3.0
2014-11-15T19:22:32.383Z
OK I now get what "Cumulative Shift of Periastron Time" means (Update 2 in my Question). If $P$ varies with $t$ then periastron generally occurs when $t/P(t)$ is not an integer and so if $M(t) = 2\pi t/P(t)$ then $M(t)$ at periastron will not be an integral multiple of $2\pi$. (Or so it seems to me).
steveOw
CC BY-SA 3.0
2014-11-15T19:32:54.610Z
@steveOW $M(t) = \int_0^t\frac{2\pi}{P(t')}\,\mathrm{d}{t'}$ ensures that if the first periastron is when $M(t) = 0$, then the next is when $M(t) = 2\pi$, etc. If $P$ is constant, then $M(t) = 2\pi t/P$; otherwise, no.
Stan Liou
CC BY-SA 3.0
2014-11-16T00:52:56.403Z
I get $M''(0) =-4\pi P'(0)/P(0)^2 \neq d(M'(0))/dt$.
steveOw
CC BY-SA 3.0
2014-11-16T04:15:04.707Z
@steveOw $\dot{M} = 2\pi/P\,\implies\,\ddot{M} = -2\pi\dot{P}/P^2$. You have an extra factor of $2$.
Stan Liou
CC BY-SA 3.0
2014-11-16T13:53:48.253Z
I expect you are right but I will write my derivation out in a separate temporary answser to discover where I went wrong:-)
steveOw
CC BY-SA 3.0
2014-11-16T18:20:54.623Z
Agreed $\dot M(0) = 2\pi/P$ but $\dot M(t) = 2\pi(P-t \dot P)/P^2$ so $\ddot M(0)$ doesn't necessarilly $=d(\dot M(0))/dt$ (I think). Please see my separate answer regarding differentiation of $M$.
steveOw
CC BY-SA 3.0
2014-11-20T17:25:14.957Z
Thanks for another excellent answer. I satisfied myself with your solution using a Binomial formula independent of Mean Anomaly (see my other answer). I'm still not happy with definition of Mean Anomaly in a non-Keplerian orbit. I might ask a separate question about it.
steveOw
CC BY-SA 3.0
2014-11-22T00:56:55.450Z
@steveOw It's not assumption; it's just vocabulary. If $P$ were constant, then $1/P$ is orbital frequency; in time interval $\Delta t$, the number of (or fraction of) orbits completed is $\Delta t/P$. $P$ actually varies continuously, but that's what calculus is for: over infinitesimal $\mathrm{d}t$, the fraction of orbit traversed is $\mathrm{d}t/P$, and integration is adding them. Thus, $\int\mathrm{d}t/P$ is the *exactly correct* orbital count, regardless of whether one wishes to call its $2\pi$ scaling 'mean anomaly' or not. If prefer not, just replace 'mean anomaly' with 'orbital phase'.
Stan Liou
CC BY-SA 3.0
2014-11-23T03:11:18.410Z
OK I get it now $1/P_i$ is like speed around an imaginary phase clock dial which starts at 0 and ends ($360^o$ later) at 1. And then the $\Delta t/P_i$'s are like distances around the circumference of that clock dial and the sum of them all in one orbit period MUST add up to 1.  Thanks for your patience!
steveOw
CC BY-SA 3.0
2014-11-25T14:24:25.370Z
(A bit pedantic but) in deriving your formula for $\Delta t$ I think you make an implicit, unstated, unverified presumption that the times of the Nth periastra in steady-period (S) and decaying-period (D) systems are approximately equal i.e. $t_{NB}\approx t_{NA}$. For the system in question it turns out the assumption is valid as the error (as a fraction of $\Delta t$) is very small $~ 2\Delta t/t_{NB}$. My (Method 3) answer (based on yours) derives a formula for $\Delta t$ independent of $t$.
steveOw
CC BY-SA 3.0
2014-11-25T23:04:57.360Z
@steveOw No, it is based entirely on whether $M(t)$ is well-approximated, and by Taylor's theorem the error has about the magnitude of the the third-order term, which I've explicitly examined to be negligible. (To be fully rigorous, the remainder term needs to be examined explicitly, but it not essentially different.) Since cumulative periastron shift is $t - NP_0$ exactly and $N = M(t)/2\pi$ exactly, all I'm doing is plugging in the approximation for $M(t)$ to substitute for $N$. No other approximations are used, not even implicitly.
Stan Liou
CC BY-SA 3.0
2014-11-26T12:42:47.747Z
The Taylor approximation error isn't the issue. You use one variable $t$. I use two variables:- (for steady system A) $t_A=NP_0$  and (for decaying system B) $t_B$ as used in $M(t)$. OK now I see that your $t$ represents $t_B$ throughout, you don't actually write $t_A$.  Previously I was mis-reading your 4th equation as being of the expected form $t_A$-$t_B$ =$\Delta t$.  Sorry.
steveOw
CC BY-SA 3.0
2014-11-15T10:17:00.620Z
2014-11-16T15:08:46.497Z
2014-11-16T15:08:46.497Z
Stan Liou
<p>As for solving the equation, $dP/dt=K/P^{5/3}$ is a separable differential equation. See the <a href="http://en.wikipedia.org/wiki/Examples_of_differential_equations" rel="nofollow">first section here</a>. Google also brings up several pages describing how to solve such equations.</p>

<p>In particular, we rewrite the equation in differential form as:
$$P^{5/3}\, dP = K\, dt.$$</p>

<p>Loosely speaking, I've cleared the denominators, with a goal to get the $P$ stuff all together on one side, and the $t$ stuff on the other. Technically it's not that at all, but pretending like it is works out. I am implicitly assuming that $K$ does not depend on $P$, which your post supports.</p>

<p>Now we integrate both sides:
$${3\over8}P^{8/3}=\int P^{5/3}\, dP = \int K\, dt,$$
where we let the right hand integral cover the constant of integration.</p>

<p>So our solution is $$P=\left({8\over3}\int K\, dt \right)^{3/8},$$ with initial condition setting the constant of integration.</p>

<p>Note I took an 8th root in finding $P$, so we need the inside function to be nonnegative.</p>

<p>Your post suggests $K$ is independent of time, in which case $\int K\, dt= K t+C$ for some constant $C$ which is determined by the initial conditions. If it does depend on time, you have to do the integration accordingly, or use data to estimate it.</p>

CC BY-SA 3.0
2014-11-14T21:46:16.440Z
Many thanks :-) As you assumed: K is independent of time and P.
steveOw
CC BY-SA 3.0
2014-11-14T21:27:47.047Z
2014-11-14T21:32:56.253Z
2014-11-14T21:32:56.253Z
zibadawa timmy
<p><strong>METHOD 3 Predicting Periastron Times and Cycle Durations from Orbital Phase ($\Phi$)</strong></p>

<p>In the answer by <strong>Stan Liou</strong> he uses a Taylor Series approximation of the <a href="http://en.wikipedia.org/wiki/Mean_anomaly" rel="nofollow">Mean Anomaly</a> to derive a nice formula which determines the CPTS (Cumulative Perihelion Time Shift) value as a function of $t^2$.  This formula produces results very close to those graphed by Weisberg &amp; Taylor.  As it took me a while to understand how the Mean Anomaly can be applied for a decaying-period orbit I thought it useful to note here what I learnt and present a slightly different way of obtaining a formula to predict periastron times.</p>

<p>Mean Anomaly basically indicates the phase of the orbit at a particular epoch, e.g. what time fraction of the orbit period has been completed.  For Mean Anomaly the fraction is scaled by $2\pi$. Let $\Phi(t)$ be the phase of the orbit at time $t$ such that at the start of the orbit $t=0$ and $\Phi(0)$ = 0 and at the end of the orbit when $t=P$ (the period of orbit) $\Phi(P)$ = 1.  I will assume that an orbit starts at one periapsis and finishes at the next periapsis.</p>

<p>In a decaying-period orbital system the value of $P(t)$ changes with $t$. Here $P$ is the period of orbit.  In a decaying-period orbit I will assume that orbital phase $\Phi(t)$ is completely specified by time $t$ and either $P(t)$ or $F(t)$, i.e. aspidial precession, progressive changes in semi-major axis length or eccentricity do not lead to significant additional changes in phase.  </p>

<p>In a decaying-period system the concept of period at a given epoch is somewhat abstract, it might be construed as a hypothetical period that would occur thereafter if the force causing period change ceased at that particular moment. Instead of period we can think in terms of $F(t)=1/P(t)$ where $F$ is the orbital frequency.  </p>

<p>I like to think of $F(t)$ as "the rate of change of phase with time", i.e. $F(t) = \dot{\Phi}(t)$.  Picture an imaginary phase clock comprising a circular dial whose circumference has markings running from $\Phi=0$ around in a full circle to $\Phi=1$. The phase $\Phi(t)$ of our subject system at any given epoch $t$ will be indicated by a marker at a particular point on the circumference of the dial. Then $F(t)$ can be thought of as the speed at which the marker is moving around the circumference of the dial at a given epoch $t$.  In a particular short interval of time $\delta t$ the (system state) marker will move a certain distance around the dial and this distance travelled will indicate the change in phase. The "travel" (change of phase) will depend on the value of $F$ during that time as per $\delta \Phi \approx \delta t * F(t)$. This is only approximate because F(t) changes during the interval $\delta t$.</p>

<p>Now let us assume that the time interval $T$ between succesive periaspes is known and that the time of the first periapsis is at $t=0$. If $P(t)$ (and therefore $F(t)$ too) changed in a step-wise manner between time intervals, then we could write</p>

<p>$$
\sum_{i=1}^{i=N=T/\delta t} \frac{\delta t}{P(t)}
=
\sum_{i=1}^{i=N=T/\delta t}  \delta t F(t)
=
\sum_{i=1}^{i=N=T/\delta t}  \delta \Phi(t)
= 1.
$$</p>

<p>To represent continuously-changing $P$ and $F$ we reduce $\Delta t$ towards zero and obtain these integral functions</p>

<p>$$
\int_{t=0}^{t=T} \frac{1}{P(t)}\, \mathrm{d}t
=
\int_{t=0}^{t=T} F(t)\, \mathrm{d}t
=
\int_{t=0}^{t=T} \dot{\Phi}(t)\, \mathrm{d}t
= 1.
$$</p>

<p>At any time $\tau$ during the orbit ($0\leq \tau \leq 1$) the current instantaneous phase is given by
$$
\Phi(\tau) = 
\int_{t=0}^{t=\tau} \dot{\Phi}(t)\, \mathrm{d}t
=
\int_{t=0}^{t=\tau} F(t)\, \mathrm{d}t
=
\int_{t=0}^{t=\tau} \frac{1}{P(t)}\, \mathrm{d}t.
$$</p>

<p>and given a constant value of $\dot P$ (the rate of change of period), and $P_o=P(0)$ we obtain
$$
\Phi(\tau) 
=
\int_{t=0}^{t=\tau} \frac{1}{P_o + \dot P t}\, \mathrm{d}t
\qquad \mathrm{and} \qquad 
\dot{\Phi}(\tau) = F(\tau) = \frac{1}{P(\tau)} = \frac{1}{P_o + \dot P \tau}.
$$</p>

<p>Using </p>

<p>$$
\frac{\mathrm{d}}{\mathrm{d}t} \left( 
\frac{1}{f(x)} 
\right)
=
\frac{-\dot f(x)}{f(x)^2}
\qquad \mathrm{and} \qquad
f(x) = P_o + \tau \dot P
\qquad \mathrm{and} \qquad
\dot f(x) =  \dot P
$$</p>

<p>we get
$$
\ddot \Phi(\tau)
= 
\frac{\mathrm{d}}{\mathrm{d}t} \left( 
\frac{1}{P_o + \tau \dot P}
\right)
=
\frac{-\dot P}{P_o^2 + (\tau \dot P)^2 + 2 P_o \tau \dot P}
$$</p>

<p>and for $\tau = 0$ we get
$$
\Phi(0) = 0
\qquad \mathrm{and} \qquad
\dot{\Phi}(0)  = \frac{1}{P_o }
\qquad \mathrm{and} \qquad
\ddot \Phi(0)
= \frac{-\dot P}{P_o^2}.
$$</p>

<hr>

<p><strong>PREDICTING PERIASTRON TIMES</strong></p>

<p>We can obtain an expression for CPTS by copying Stan Liou's approach , but using phase ($\Phi$) instead of Mean Anomaly ($M$).</p>

<p>Write the phase at time $\tau$ as a Maclaurin series (a Taylor series with $a=0$):
$$\begin{eqnarray*}
\Phi(\tau) \equiv \int_{t=0}^{t=\tau} \dot{\Phi}(t)\,\mathrm{d}t 
&amp;=&amp; \Phi(0) + \dot{\Phi}(0) \tau + \frac{1}{2}\ddot{\Phi}(0) \tau^2 + \mathcal{O}(\tau^3)\\
&amp;=&amp;  0 +\frac{1}{P_0}\tau - \frac{\dot{P}}{2P_o^2}\tau^2+\mathcal{O}(\tau^3)\text{.}
\end{eqnarray*}$$</p>

<p>Dropping the third order and residual error terms (Stan Liou's answer explains justification for this) we get</p>

<p>$$
\Phi(\tau) \approx \frac{1}{P_0}\tau- \frac{\dot{P}}{2P_o^2}\tau^2.
$$</p>

<p>Let $N$ count the number of completed orbits at some time $\tau_N$. Then the phase at the end of orbit number $N$ is given by $\Phi(\tau_N) = N$. In a steady system (A) the orbital period $P_o$ stays constant $\dot P=0$ and so the time $\tau_{NA}$ at which the $N$th orbit completed is given exactly by $\tau_{NA} = NP_o$ and thus $N=\tau_{NA}/P_o$.  In a decaying-period system (B) the Nth orbit completes at time $\tau_{NB}$ when $\Phi(\tau_{NB}) =  N$ but with $\Phi(\tau_{NB})$ depending on a non-zero value of $\dot P$. </p>

<p>$$
\begin{align}
&amp;\Phi(\tau_{NA}) = N = \frac{1}{P_0}\tau_{NA}\\
&amp;\Phi(\tau_{NB}) = N \approx \frac{1}{P_0}\tau_{NB}- \frac{\dot{P}}{2P_o^2}\tau_{NB}^2.
\end{align}
$$</p>

<p>We can obtain a formula for $\Delta t_{N}$ which is the difference in time between $\tau_{NB}$ (the epoch of the decaying system $N$th periastron) and $\tau_{NA}$ (the epoch of the steady system $N$th periastron) as follows.  Note that following Weisberg &amp; Taylor we define $\Delta t_{N}\,(=$ CPTS $)=\tau_{NB} - \tau_{NA}$. The $N$th periastron in the decaying period system occurs earlier than the $N$th periastron in the steady system. Therefore   $\Delta t_{N}$ (=$ CPTS $)  will become increasingly negative as time passes (as $t$ increases).
$$
\frac{\tau_{NA}}{P_o} = N \approx  \frac{\tau_{NB}}{P_0}- \frac{\dot P}{2P_o^2} \,{\tau_{NB}}^2
$$</p>

<p>$$
{\tau_{NA}}  \approx  {\tau_{NB}} - \frac{\dot P}{2P_o} \,{\tau_{NB}}^2
$$</p>

<p>$$
\Delta t_N 
= \tau_{NB} - \tau_{NA} \approx + \frac{\dot P}{2P_o} \,{\tau_{NB}}^2
$$</p>

<p>which is Stan Liou's approximation for $\Delta t_N$.  </p>

<p>So what value is this equation?
In general assume that we have determined the time $t_0$ of the $0$th periastron and measured accurately an initial orbital period $P_o$ and we keep track of the observed periastra. </p>

<p>Case 1 - we observe that the $N$th periastron occurs at a particular time ($\tau_{NB}$). Using $\tau_{NA}=N P_o$, we can easily calculate $\Delta t_N = \tau_{NB} - \tau_{NA}$.  Then, using Stan Liou's approximation of $\Delta t_N$ we can obtain an empirical estimate of $\dot P$ from
$$
{\dot P} = \frac{2 P_o \,\Delta t_N} {\tau_{NB}}^2.
$$</p>

<p>Case 2 - we have a theoretical formula which, given our measurement of $P_o$ predicts the value of  $\dot P$.  Given various values of $N$ and using $\tau_{NA}=N P_o$ we can easily calculate the hypothetical epochs ($\tau_{NA}$ for various $N$) of the hypothetical steady-system periastra.  Now we wish to plot a curve showing the theoretical values of $\Delta t_N$ at various hypothetical periastron epochs ($\tau_{NA}$).  But Stan Liou's equation for $\Delta t_N$ uses $\tau_{NB}$ not $\tau_{NA}$.  We need to find $\Delta t_N$ as a function of $\tau_{NA}$.  Therefore we would need to solve the following quadratic equation for $\Delta t_N$ in terms of $\tau_{NA}$ and $K$, where $K=\frac{\dot P}{2 P_o}$:</p>

<p>$$
\Delta t_{N} = K {\tau_{NB}}^2 = K ({\tau_{NA}}^2 + 2 \tau_{NA}\Delta t_{N} +{\Delta t_{N}}^2)
$$</p>

<p>$$
 0 =  (K{\tau_{NA}}^2 + 2K \tau_{NA}\Delta t_{N} +K{\Delta t_{N}}^2) -\Delta t_{N}
$$</p>

<p>$$
 0 =  (K){\Delta t_{N}}^2    +(2K \tau_{NA} -1)\Delta t_{N}      +(K{\tau_{NA}}^2) 
.$$</p>

<hr>

<p>The following way of getting an equation for $\Delta t$ is more protracted but it will provide a formula for $\Delta t$ as a function of $N$ (which can easily be converted into a function of $\tau_{NA}$).</p>

<p>Starting with 
$$
\frac{\tau_{NA}}{P_o} =\frac{1}{P_0}\tau_{NB}- \frac{\dot P}{2P_o^2} \,{\tau_{NB}}^2
$$</p>

<p>$$
\tau_{NA} =\tau_{NB}- \frac{\dot P}{2P_o}{\tau_{NB}}^2
\qquad\rightarrow \qquad
 0 = -\tau_{NA} + \tau_{NB} - \frac{\dot P}{2P_o} \,{\tau_{NB}}^2
$$</p>

<p>Using $ x = \frac{-b \pm \sqrt{b^2-4ac}}{2a}$ with $a=- \frac{\dot P}{2P_o} , b=1, c=-\tau_{NA}$ we obtain
$$
\tau_{NB} = \frac{-(1) \pm \sqrt{(1)^2-4(- \frac{\dot P}{2P_o}).(-\tau_{NA})}}{2(- \frac{\dot P}{2P_o})}
$$</p>

<p>$$
\tau_{NB} = \frac
{1 \pm \sqrt {1 - 2 \frac{ \dot P  }{P_o}.\tau_{NA}     }}
{\dot P  /{P_o}}
\qquad \longrightarrow \qquad
\tau_{NB} = 
 \frac{P_o}{ \dot P } \left( 1 \pm \sqrt {1 - 2 \frac{\dot P  }{P_o}.\tau_{NA}     }\,\, \right)
$$</p>

<p>Now $\tau_{NA} = P_o \,N$ so we can write
$$
\tau_{NB} = 
\frac{P_o}{ \dot P } \left(1 \pm \sqrt {1 - 2  \dot P  N } \right)
\qquad \longrightarrow \qquad
-\Delta t = \tau_{NA}-\tau_{NB} 
= P_o \left( N -\frac{1}{ \dot P}  \left(1 \pm \sqrt {1 - 2  \dot P  N } \right) \right)
$$</p>

<p>= = = = = = =</p>

<p>We can try a simple approximation of this using $(1+\alpha)^{0.5} = (1+ 0.5\alpha)$ so that $\sqrt {1 - 2  \dot P   N } \approx 1- \dot P  N$ and by inspection the "$\pm$" becomes a "$+$" 
$$
-\Delta t 
= P_o \left( N -\frac{1}{  \dot P  } \left(1 \pm 1 - \dot P   N  \right) \right)
= P_o ( N - N ) 
= 0
$$</p>

<p>This is a valid approximation but it goes a bit too far for our purposes! We know that $\tau_{NA}-\tau_{NB} \neq 0$. So let us try expanding the term $[1 -2 \dot P  N]^{0.5}$ in a series expansion,</p>

<p>$$
[1 + (-2 \dot P  N)]^{0.5} 
= 1 +\frac{1}{2}(-2\dot P  N)^{1} -\frac{1}{8}(-2 \dot P  N)^{2} 
+\frac{1}{16}(-2\dot P  N)^{3}+...
$$</p>

<p>$$
[1+ (-2\dot P  N)]^{0.5} 
= 1 - \dot P N                    -\frac{1}{2}(\dot P  N)^{2} 
-\frac{1}{2}( \dot P  N)^{3}+...
$$</p>

<p>Let us limit ourselves to terms with $ \dot P  $ to the power of $2$ or less thus</p>

<p>$$
[1-2 \dot P  N]^{0.5} \approx 
1 -  \dot P  N -\frac{1}{2}\dot P ^2 N^2 .
$$</p>

<p>We can insert this into the time difference equation</p>

<p>$$
-\Delta t
\approx P_o \left( N -\frac{1}{ \dot P  } 
\left(1 \pm (1 - \dot P N -\frac{1}{2} \dot P ^2N^2 \right)  \right) 
$$</p>

<p>by inspection the "$\pm$" becomes a "$-$"</p>

<p>$$-\Delta t\approx P_o \left( N -\frac{1}{ \dot P  } 
\left(  + \dot P  N +\frac{1}{2}\dot P  ^2 N^2 \right) \right) 
\qquad \longrightarrow \qquad
-\Delta t\approx P_o \left( N - \left( N +\frac{1}{2} \dot P N^2) \right) \right)
$$</p>

<p>$$
-\Delta t \approx -\frac{1}{2}P_o \dot P  \, N^2
\qquad \longrightarrow \qquad
\Delta t \approx \frac{1}{2}P_o \dot P  \, N^2
$$</p>

<hr>

<p>From the Weisberg &amp; Taylor paper we are given  $P_o=27906.979587552 s$ and $\dot P  =  -2.40242 *10^{-12} s/s$ 
so $\frac{1}{2}P_o  \dot P  = -0.5 * 27906.979587552 * 2.40242 *10^{-12}  = -3.352216747 * 10^{-08} s$.</p>

<p>With $N= 10,000$ cycles we get $\Delta t = -3.352217 s$.</p>

<p>Note that we can make the substitution ${\tau_{NA}}^2 = N^2 Po^2$ into the equation for $\Delta t$ to give
$$
\Delta t \approx \frac{1}{2}P_o \dot P  \, N^2  = \frac{1}{2} \frac {\dot P}{P_o}  \,{\tau_{NA}}^2
.$$
The difference between this $\Delta t$ and the value of $\Delta t$ from Stan Liou's equation is
$$
\frac{1}{2} \frac {\dot P}{P_o}  \left( {\tau_{NA}}^2 -{\tau_{NB}}^2   \right)
$$
$$
= 
\frac{1}{2} \frac {\dot P}{P_o}  \left( 
{{\tau_{NA}}^2 
-\tau_{NA}^2 +2 \tau_{NA} \Delta t + \Delta t^2}   \right)
$$
$$
= 
\frac{1}{2} \frac {\dot P}{P_o}  \left( 
{2 \tau_{NA} \Delta t + \Delta t^2}   \right)
$$
and the difference as a fraction of $\Delta t$
$$
= 
\frac{1}{2} \frac {\dot P}{P_o}  \left( 
{2 \tau_{NA}  + \Delta t}   \right)
.$$</p>

<p>For the present case $(\dot P/P_o \approx 8.6 * 10^{-17})$ and so, after 34,000 periastron cycles (948,838,000 s, about 30 years) the difference between the two approximations of $\Delta t$ would be only $0.8 * 10^{-8} s$ which is much smaller than the time resolution of the observations.</p>

<hr>

<p><strong>PREDICTING CYCLE DURATIONS</strong></p>

<p>We can also derive an expression for the cycle duration $D$, as follows.</p>

<p>The start periaston is at $t=0=T0$ and the first following periaston is at $t=T1$ then
$$
\Phi(T1)=\int_{t=0}^{t=T1}\frac{1}{P(t)}\,\mathrm{d}t = 1
$$
this can be expressed as 
$$
1 =\int_{t=0}^{t=T1}\frac{1}{ \dot{P} t + P_0 }\,\mathrm{d}t
= \left[ \frac{1}{\dot{P}} ln  ( C\dot{P}t + CP_0 )\right]_0^{T1}
$$
hence
$$
\dot{P} = 
ln ( C\dot{P}T1 + CP_0 ) - ln  ( CP_0 )
=
ln \left(   \frac{ C ( \dot{P}T1 + P_0 )} { C ( P_0 )}   \right)
=
ln \left(   \frac{ \dot{P}T1 + P_0 } {P_0 }       \right)
$$</p>

<p>$$
\rightarrow exp(\dot{P}) 
= \frac{ \dot{P}T1 + P_0 } {P_0 }
\rightarrow
P_0 exp(\dot{P}) -P_0
= \dot{P}T1 
$$
giving us
$$
\rightarrow
T1 =
P_0 \frac{exp(\dot{P}) -1} {\dot{P}}
$$</p>

<p>so the duration of the first orbit (from start periastron to first following periastron) is $D_1 = T1-T0 = T1$ thus
$$
D_1 =
D_0 \frac{exp(\dot{P}) -1} {\dot{P}}
$$</p>

<p>and we can generalize this to
$$
D_{N} =
D_{N-1} \left( \frac{exp(\dot{P}) -1} {\dot{P}} \right)
.$$</p>

<p>For example using the made-up value $\dot P = -2.34\,10^{-8}$ we obtain
$$
D_{N+1} =
D_{N} \frac{exp(-2.34*10^{-5}) -1} {-2.34*10^{-8}}
= 0.999 999 988 D_{N}
$$</p>

<p>However, using standard arithmetical software (such as Excel) when we try to calculate using $\dot P = -2.34\,10^{-10}$ we get nonsense results because of truncation errors. The approximate value of $\dot P$ reported for the Hulse-Taylor system is about $-2.34\,10^{-12}$.</p>

<p>We can analyze the formula for $D_{N+1}$ using series expansions.  A Taylor Series expansion of $exp(x)/x -1/x$  is given by WolframAlpha as</p>

<p>$$
exp(x)/x -1/x= 1 + 
  \frac{x}{2}
+ \frac{x^2}{6}
+ \frac{x^3}{24}
+ \frac{x^4}{120}
+ \frac{x^5}{720}
+ \frac{x^6}{5040}
+ \frac{x^7}{40320}
+ \frac{x^8}{362880}
+ \frac{x^9}{3628800}
+ \frac{x^{10}}{39916800}
+O(x^{11})
$$
or 
$$1 + \frac{x^{2-1}}{2} +  \frac{x^{3-1}}{3*2} + ... \frac{x^{i-1}}{i!} + ...$$
No easilly computable expression or approximation is obvious at present.
Proceeding anyway, it is clear that the duration of any subsequent orbit $N$ can be computed from </p>

<p>$$
D_N 
= 
D_0 \left(\frac{exp^{\dot P} -1} {\dot P} \right) ^N 
=
D_0 \left(\frac{1}{\dot P} \right)^N \left(exp^{\dot P} -1 \right) ^N 
$$</p>

<p>It is interesting to compare this expression for cycle duration with that which was used as the basis of the coarse binomial solution (see other answer: Method 4)</p>

<p>$$D_1=D_0(1+\dot P)^1$$</p>

<p>The ratio of durations Phase-based to Coarse binomial based is
$$
D_0 \left( \frac{exp(\dot P)-1}{\dot P} \right)^{1}
:D_0(1+\dot P)^1
$$</p>

<p>becoming
$$
\frac{exp(\dot P)-1}{\dot P} 
:1+\dot P
\qquad \rightarrow
exp(\dot P)-1
:\dot P +\dot P^2
\qquad \rightarrow
exp(\dot P)
: 1+\dot P +\dot P^2
$$</p>

<p>applying the series expansion of $exp(\dot P)$ we obtain
$$
\qquad \rightarrow
1+\dot P + \frac{\dot P^2}{2} + \frac{\dot P^3}{6} + \frac{\dot P^4}{24} +\,...\,
: 1+\dot P +\dot P^2
$$</p>

<p>Clearly the two expressions give different values for cycle duration $D$ and the difference appears at the term in $\dot P^2$ and seems to be no bigger than $\frac{\dot P^2}{2}$.</p>

CC BY-SA 3.0
2014-11-24T01:55:43.877Z
I do not think you should do this way
Py-ser
CC BY-SA 3.0
2014-11-26T13:28:21.930Z
@Py-ser. Please could you elaborate? (Note there have been major edits since your comment).
steveOw
CC BY-SA 3.0
2014-11-26T13:40:28.663Z
I meant, to edit more answers instead of one. I did not criticize your actual answer, but the way of answering.
Py-ser
CC BY-SA 3.0
2014-11-26T17:00:14.883Z
@Py-ser. Thanks for the feedback. It was certainly messy yesterday.  Hopefully now it is a bit more simple. I have changed  the question to signpost the different methods.  I am reluctant to merge my two answers into one because they are different methods with different terms employed (even though they produce very similar results).
steveOw
CC BY-SA 3.0
2014-11-24T00:15:39.560Z
2014-11-27T13:46:16.533Z
2014-11-27T13:46:16.533Z
steveOw
<p><strong>METHOD 4:Predicting CPTS using a Coarse Binomial Series approximation.</strong></p>

<p>This method does not make any assumptions about Mean Anomaly. Orbit is defined as the cycle from one periapsis to the next periapsis. </p>

<hr>

<p><strong>SUMMARY</strong></p>

<p>The time from the $0$th to the $N$th periastron can be obtained from the following sum, from which a binomial expression can be derived</p>

<p>$$\frac{T_N}{P_0} 
\qquad = \sum_{i=1}^{i=N} (1+\dot P)^i
\qquad = N + \sum_{j=2}^{N+1}\binom{N+1}{j}\dot P^{j-1}.$$</p>

<p>For our purposes terms in $\dot P$ greater than first order can be dropped leading to </p>

<p>$$ T_{N} = N{P_0} + P_0 \frac{N^2+N}{2} {\dot P} .$$</p>

<p>The time difference between the $N$th periapses of a steady period system and a decaying period system is then</p>

<p>$$ \Delta t_{N} =  0.5 P_0 {\dot P}  (N^2+N) $$</p>

<p>This gives results which agree closely with the Wesiberg &amp; Taylor graph and the formula in Stan Liou's answer.</p>

<hr>

<p><strong>Definitions and Assumptions</strong></p>

<p>Let us consider the variation over time in the <strong>period</strong> of orbit. Weisberg &amp; Taylor provide a formula for the variation in "period" with time.  However the meaning of "period at a point in time" is ambiguous when the period is changing with time. We therefore have to define how "period at a point in time" translates into a time interval between the start and end of a particular orbital cycle (e.g. from periapsis to periapsis).  In this method we will "fudge" the translation by assuming (A1) that the time interval between successive periapses is equal to the "hypothetical period at the point in time of the end of the cycle". Also to "boot" the model we will assume (A2) that at time zero ($t=0$) the hypothetical period is $H_o$ and thus the hypothetical period at the end of the first cycle is given by $H_1 = H_o + H_o \dot H$.  Therefore, by (A1) we have the following definition of the duration of the first cycle $D_1 = H_1 = H_o + H_o \dot H$.    </p>

<p>We further assume (A3) that, over the period of consideration (~ 40 years) the rate of change of period with time is effectively constant ($\ddot H = 0$ and  $\dot H(t) = constant = \dot H$ ). (See Stan Liou's answer for justification). Thus we have
$$
H(t) = H(0) + t \dot H .
$$
It is emphasized that these assumptions are not strictly mappable to a consistent physical model.  However given that we are dealing with very small changes in period over a very large number of orbits the ruleset may be useful.  We shall proceed and see what comes out.  Henceforth I will not distinguish between hypothetical period $H$ and cycle duration $D$.  I will assume that they are both represented by the same term period $P$ and that the duration $P$ of a particular periapsis to periapsis cycle is equal to the hypothetical period $P(t)$ where $t$ is the epoch of the end of the cycle.  Also $\dot P$ will represent both the rate of change of hypothetical period and also the rate of change of cycle duration between cycles.</p>

<p>$P_N$ is the period of the $N$th orbit (orbit $O_N$). Orbit number $0$ terminates at $t=T_0=0$ at perihelion number $0$ and has period $Po = P(0)$.   The first considered orbit ($O_1$) starts at time $t=0$ at perihelion number $0$ and finishes at perihelion number $1$ at time $t=T_1$ and has period duration $P_1 = P_0+ P_0\dot P =P(T_1)$. Therefore we can write</p>

<p>$$P_1 = P_0(1+\dot P)^1$$
$$P_2 = P_1(1+\dot P) = P_0(1+\dot P)(1+\dot P)= P_0(1+\dot P)^2$$
$$P_3 = P_2(1+\dot P) = P_0(1+\dot P)(1+\dot P)(1+\dot P)= P_0(1+\dot P)^3$$
In general then, the period of any orbit is given by the formula:
$$P_N = P_0(1+\dot P)^N.$$</p>

<p>Now let us consider the <strong>time</strong> $T$ at which each orbit ends. We have 
$$T_1 = T_0+P_1 = P_1            = P_0(1+\dot P)$$ 
$$T_2=T_1+P_2 = P_1+P_2  = P_0(1+\dot P)+P_0(1+\dot P)^2$$
$$T_3=T_2+P_3 = P1+P2+P3 = P_0(1+\dot P)+P_0(1+\dot P)^2+P_0(1+\dot P)^3$$
$$T_N=P_1+P_2+...+P_N = P_0(1+\dot P)+P_0(1+\dot P)^2+...+P_0(1+\dot P)^N.$$</p>

<p>Thus
$$T_N = {P_0}\sum_{i=1}^{i=N} (1+\dot P)^i\, .$$</p>

<p>(The Note at the end of this answer explores a way of proceeding which uses an exact equivalent of this expression in a non-serial form.  However, that route is not productive due to difficulties with exponentials.)</p>

<hr>

<p>We shall proceed as follows,  (using $k$ in place of $\dot P$)
$$ \frac{T_N}{P_0}= S_N = \sum_1^N (1+k)^i $$</p>

<p>We can multiply out the first few terms and examine the patterns in the following <strong>pyramid</strong>...
$$(1+k)^1 =               k +1  $$
$$(1+k)^2 =        k^2  +2k +1  $$
$$(1+k)^3 = k^3  +3k^2  +3k +1  $$
$$(1+k)^4 =k^4+4 k^3+6 k^2+4 k+1$$
$$(1+k)^5 =k^5+5 k^4+10 k^3+10 k^2+5 k+1$$
$$(1+k)^6 =k^6+6 k^5+15 k^4+20 k^3+15 k^2+6 k+1$$</p>

<p>So for example, for $N=3$ we would obtain the sum
$$
S_3= k^3  +4k^2 +6k +1
$$</p>

<p>The results suggest a solution with a pattern of the form
$$
S_N = a + bk^1 +ck^2+dk^3...
$$</p>

<p>We can see straightaway that $a=N$.  The other coefficients increase with $N$. It might be possible to determine a formula for the coefficients from the pattern. Although the general pattern is not convergent, it is possible that for certain restricted ranges of $N$ and $k$ a convergent formula could be obtained.  If so then it might be possible to obtain a useful approximation of $S_N$  by evaluating the first few terms in the series.</p>

<hr>

<p>The following binomial equation was derived in an answer to <a href="https://math.stackexchange.com/questions/1029172/how-do-i-express-the-sum-1k1k21k3-1kn-for-k1-as-a-se">this question</a> at Maths Stack Exchange by User73985 </p>

<p>$$ \sum_1^N (1+k)^i \qquad = \qquad N + \sum_{j=2}^{N+1}\binom{N+1}{j}k^{j-1}$$
So
$$ S_N = N + \sum_{j=2}^{N+1}\frac{(N+1)!}{(N+1-j)!j!}k^{j-1}$$</p>

<p>then
$$ S_N = N 
+  \frac{(N+1)!}{(N-1)!2!} k^{1} 
+  \frac{(N+1)!}{(N-2)!3!} k^{2}
+  \frac{(N+1)!}{(N-3)!4!} k^{3} +...
$$</p>

<p>giving
$$ S_N = N 
+  \frac{(N+1)(N)}{2!} k^{1} 
+  \frac{(N+1)(N)(N-1)}{3!} k^{2}
+  \frac{(N+1)(N)(N-1)(N-2)}{4!} k^{3} +...
$$</p>

<p>thus
$$ S_N= N 
+  \frac{N^2+N}{2} k^{1} 
+  \frac{N^3-N}{6} k^{2}
+  \frac{N^4-2 N^3-N^2+2 N}{24} k^{3} +...
$$
(Switching back from using $k$ to using $\dot P$...)
For $N= 10,000$ and ${k=\dot P}= 2.40242 * 10^{-12}$ the ratio of each term to the next term is greater than $1 : N{\dot P} \approx 1:2.4 * 10^{-8}$ so for our purposes this formula can be approximated by truncating to the first order power of ${\dot P}$ to give
$$ S_{N} = \frac{T_{N}}{P_0} = N + \frac{N^2+N}{2} {\dot P}^{1} $$
So
$$ T_{N} = N{P_0} + P_0 \frac{N^2+N}{2} {\dot P} $$
and so the time difference between the Nth periapses (fixed $P$ - decaying $P$) is</p>

<p>$$ \Delta t_{N} = P_0{\dot P} \frac{N^2+N}{2}  $$
and for $N=10,000$ we obtain (from Excel) $\Delta t = 3.35255 s$ which is very close to the value $(3.35221 s)$ obtained from Stan Liou's formula ($\Delta t= 0.5 {\dot P}/P_0 t^2$).</p>

<p>Over 30 years (33,924.42 cycles) $\Delta t$ = 38.5806 $s$ which is very close to the amount indicated by the graph and by Stan Liou's formula $(38.5794 s)$.</p>

<hr>

<p><strong>Note</strong></p>

<p>WolframAlpha indicates the equality</p>

<p>$$T_N = {P_0}\sum_{i=1}^{i=N} (1+\dot P)^i
\qquad =
\frac{P_0}{\dot P}
\,
(1+\dot P)[(1+\dot P)^N -1]
$$</p>

<p>For small $\dot P$ this can be approximated by:-
$$T_N =\frac{P_0}{\dot P}\,(1+\dot P)[(1+N\dot P) -1]
=\frac{P_0}{\dot P}\,(1+\dot P)(N\dot P) 
=N{P_0}\,(1+\dot P). 
$$
The end time of $N$ cycles of period $P_{constant}=P_0$ will be $T_{NC}=N P_0$. And so the time difference between the end times of  (i) N cycles of constant $P$ and (ii) $N$ cycles of decaying $P$ is approximated by 
$$\Delta t = NP_0 - NP_0(1-\dot P) = NP_0 \dot P.
$$
For  $P_0$ = 0.322997448930 d = 27906.979587552 s and $\dot P$ = -2.40242 x 10$^{-12}$ we obtain $(P_0\,dot P) =$ -6.70442859007267$^{-08}$. For  $N=10,000$ cycles (approximately 10 years) this indicates a time difference of 6.70442859007267 $^{-03}$ s, about $7 ms$.  This time difference is very small compared to the value roughly indicated by the original graph and the value calculated by Stan Liou's formula, from which we would expect a time difference of about 3$s$.  The explanantion is that the first-order approximation is not suitable in this case because the true value of $\Delta t$ is so small compared to $t$.</p>

CC BY-SA 3.0
2014-11-20T17:17:37.053Z
2014-11-25T14:13:31.780Z
2017-04-13T12:19:35.350Z
steveOw
<p>Following on from a <a href="https://astronomy.stackexchange.com/questions/6051/how-did-the-hulse-taylor-pulsar-provide-confirming-evidence-of-general-relativit">previous question</a> about the interesting <a href="https://en.wikipedia.org/wiki/Hulse%E2%80%93Taylor_binary" rel="nofollow noreferrer">Hulse Taylor Binary Pulsar</a>.</p>

<p><strong>BACKGROUD</strong></p>

<p><a href="http://arxiv.org/pdf/astro-ph/0407149v1.pdf" rel="nofollow noreferrer">Weisberg &amp; Taylor, 2004</a> present a graph showing the change in "Cumulative Periastron Time Shift" over a 30 year period.  Here is a later version of that graph
<img src="https://i.sstatic.net/992mC.gif" alt="orbit period decay"></p>

<p>The graph indicates excellent agreement between the observed data points and the theoretical curve predicted by the General Theory of Relativity (GTR).  </p>

<p>The authors determined the theoretical change in orbit period <span class="math-container">$P$</span> over time (due to orbital energy loss from emission of gravitational waves) from a formula (Peters
&amp; Matthews (1963)) of the form:-
<span class="math-container">$$ dP/dt = K/P^{(5/3)}$$</span> 
where K is an (effectively constant) function of the masses of the two stars and the orbit eccentricity. For the Hulse-Taylor Binary, the period <span class="math-container">$P$</span> is approximately 7.75 hours and the value of K is approximately <span class="math-container">$(-6 * 10^{-5})$</span>, giving an approximate value of <span class="math-container">$$ dP/dt = 2.34 × 10^{−12} s/s. $$</span></p>

<p>If <span class="math-container">$dP/dt$</span> is constant then obviously we can determine the period <span class="math-container">$P$</span> at a future time by applying the formula
<span class="math-container">$$P_t = P_0 + t\,(dP/dt).$$</span>
As a rough initial estimate over one year the approximate change in period is: 365.25 x 24 x 3600 x -2.34 x 10^{-12} = 0.074 ms. (Wikipedia article indicates 0.076 ms.)</p>

<p>But <span class="math-container">$dP/dT$</span> is not constant because it depends on the value of <span class="math-container">$P$</span> which obviously changes.  So we would expect the period to decay over time and the rate of decay to increase with time. However, the change in P is currently so very small that over the period of 30 years the change in dP/dt (which is a function of <span class="math-container">$P$</span>) is extremely small also. (It is important to appreciate that the graph shown above does <strong>not</strong> plot <span class="math-container">$dP/dt$</span> against time).</p>

<p><strong>zibadawa timmy</strong> (thanks) provided the formula <span class="math-container">$P=(8/3 \int K dt)^{3/8}$</span> giving:- </p>

<p><span class="math-container">$$
P_t = \left( \frac{8}{3} (Kt+C)\right)^{3/8} \qquad \mathrm{with}\qquad C=\frac{3}{8} P_o^{8/3}
$$</span></p>

<p>hence</p>

<p><span class="math-container">$$
P_t = 
 \left( \frac{8}{3} Kt+ P_o^{8/3})\right)^{3/8} 
=
 \left(  \frac{-1.6\,t}{10^4}+ P_o^{8/3}\right)^{3/8} 
=
 \left(  \frac{-1.6\,t}{10^4}+ 27900^{8/3}\right)^{3/8} 
$$</span>
for 1 year , t = 31,557,600 s., so -1.6 t/10^4 = -5,049.12 = c. 5e3</p>

<p>for 10 years t = 310,557,600 s., so -1.6 t/10^4 = -50,491.2 = c. 5e4</p>

<p>while 27900^{8/3} = 716,051,966,959.181 = c. 7e11.</p>

<p>So 
<span class="math-container">$$
P_{1 year} = (-5049.12 + 716,051,966,959.181)^{3/8} = 27,899.999 926 s,\qquad deltaP_{1 year} = -0.0000738 s.
$$</span>
<span class="math-container">$$
P_{10 years} = (-50491.2 + 716,051,966,959.181)^{3/8} =  27,899.999 262 s,\qquad deltaP_{10 years} = -0.000738 s.
$$</span>
These <span class="math-container">$deltaP$</span> values indicate that the value of <span class="math-container">$dP/dt$</span> is effectively constant ($=-0.000 073 8 s/s) over 10years.</p>

<hr>

<p><strong>QUESTION 1 - What is plotted on the graph? (<span class="math-container">$\equiv$</span> What is Cumulative Perihelion Time Shift?) (ANSWERED).</strong> </p>

<p>Initially I was confused about what the Weisberg &amp; Taylor graph was actually plotting.  The 2004 paper has "The observed and theoretical decays are compared graphically in Figure 1".  The wikipedia article labels the y-axis "Cumulative Period Shift" but the caption reads "observed change in the epoch of periastron with date".</p>

<p>An earlier <a href="http://www.nature.com/scientificamerican/journal/v245/n4/pdf/scientificamerican1081-74.pdf" rel="nofollow noreferrer">1981 Scientific American article</a> indicates that
the vertical axis of the graph indicates the difference (round data points) between (A) the calculated epoch of a periastron event in a hypothetical "steady" system whose orbital period remains constant and (B) the observed epoch of the corresponding actual periastron event in the "decaying"-period system. Also shown is the continuous curve showing the difference between (A) steady system periastron epochs and (C) theoretical predictions of decayed periastron epochs.   The periastron's of the actual and hypothetical orbits are synchronized at <span class="math-container">$t=0$</span> (in late 1974).</p>

<hr>

<p><strong>QUESTION 2 - How are the theoretical values of CPTS Calculated? (ANSWERED).</strong></p>

<p>Several (similar) methods have now been identified.</p>

<p>Method 1 The method(s) used by Weisberg &amp; Taylor to calculate theoretical values of CPTS have not been identified. These following methods seem reasonable candidates, based on the agreement with their results.</p>

<p>Method 2 (see Stan Liou's answer)
<strong>Stan Liou</strong> has provided a formula <span class="math-container">$\Delta t = 0.5(\dot P/P_0)t^2$</span> that produces values which closely match the Weisberg &amp; Taylor graph.  </p>

<p>Method 3 (see my answer Method 3)
I have derived a formula identical to Stan Liou's but derived using phase instead of mean anomaly.  Another formula is derived from the same trunk expressing <span class="math-container">$\Delta t$</span> in terms of <span class="math-container">$\dot P, P and N$</span> <span class="math-container">$\Delta t = -0.5 P_o \dot P N^2$</span>.</p>

<p>Method 4 (see my answer Method 4)
I have derived a formula derived using a coarsely approximated  binomial series which gives similar results 
<span class="math-container">$ \Delta t_{N} = P_0{\dot P} \frac{N^2+N}{2}  $</span>.  </p>

<p>Method 5 (not shown)
A finely approximated series expansion approach has been assessed. For the given binary system it produces a formula that is effectively indistinguishable from that of the other methods for the given time range.</p>

CC BY-SA 3.0
2014-11-14T19:51:59.850Z
As for solving the equation, $dP/dt=K/P^{5/3}$ is a separable differential equation with solution $$P=\left({8\over3}\int K\, dt \right)^{3/8},$$ with initial condition setting the constant of integration.
zibadawa timmy
CC BY-SA 3.0
2014-11-14T21:05:50.223Z
@zibadawa timmy. Many thanks. That appears to confirm my computer results. If possible please can you point me to the general method for your solution?
steveOw
CC BY-SA 3.0
2014-11-14T21:50:01.790Z
For the remainder of your issue, is $t$ measured in seconds or years? It sounds like seconds. But the data points are in years. 1975 to 1985 is a lot more than 10 seconds. In particular it looking linear for 20 seconds, or even a year, is no indication it is still linear an entire decade of seconds later.
zibadawa timmy
CC BY-SA 3.0
2014-11-14T21:54:58.297Z
I'm on a mobile, and the backspace is really close to the post button. So I hit the wrong one. I've edited it to finish the thought.
zibadawa timmy
CC BY-SA 3.0
2014-11-14T22:09:32.377Z
I have been using seconds.  Over 10 years 1975-1985 my (+your) formulae indicate relatively linear - from end Year 1 to end Year 10 deltaP increases by x10.  But the source graph indicates much stronger curve (you might not see it on a mobile!).
steveOw
CC BY-SA 3.0
2014-11-14T22:09:44.583Z
There are over 315 million seconds in a decade. Just eyeballing the equation, and assuming $t=0$ is 1975, tells me it's curving over that huge of an interval.
zibadawa timmy
CC BY-SA 3.0
2014-11-14T22:48:07.277Z
I have fleshed out the equations and deltaP derivations at the end of the Question.
steveOw
CC BY-SA 3.0
2014-11-15T01:00:30.173Z
You seem to be supposing that $P$ and the cumulative periastron time shift are the same. The periastron is the point of closest approach in an orbit (says Google). $P$ is the orbital period. These should be related (the "cumulative" part is especially important), but it seems unlikely that they are the same. So I think you need a different formula to get that graph.
zibadawa timmy
CC BY-SA 3.0
2014-11-15T01:44:51.133Z
Agreed. I'm assuming the graph should plot deltaP(t) vs. t, where deltaP(t) = P(t)-P(0) and t=0 in year 1975 AD. I cant presently think what "Cumulative Periastron time shift" should be if it is not same as deltaP.
steveOw
CC BY-SA 3.0
2014-11-15T01:52:21.847Z
Ultimately they are different things by their definition. $P$ is orbital period. Periastron is a point of closest approach in one orbit. The periastron time is how long into an orbit until you are at periastron (requires you pick a start point). And so the periastron time shift is how much the periastron time has changed between complete orbits (this does not actually depend on the  start time). Clearly a change in orbital period induces a periastron time shift, but it seems improbable that the shift equals the period change.
zibadawa timmy
CC BY-SA 3.0
2014-11-15T01:58:10.937Z
As for what the actual equation is, I do not know. One would hope it is mentioned in the paper. I am not familiar with this topic so I do not know the equation myself. I'm making educated guesses on the exact definitions of the terms, as it is.
zibadawa timmy
CC BY-SA 3.0
2014-11-15T02:15:52.647Z
Let us [continue this discussion in chat](http://chat.stackexchange.com/rooms/18669/discussion-between-steveow-and-zibadawa-timmy).
steveOw
CC BY-SA 3.0
2014-11-16T15:11:42.560Z
I've followed up on the PPN calculations for the $\dot{P} = KP^{-5/3}$ formula, and while your $K$ is a function of $e$, $e$ is itself not constant (should have checked sooner). But this makes very little difference because the rate of change of $e$ turns out to be quite negligible.
Stan Liou
CC BY-SA 3.0
2014-11-16T21:55:16.470Z
@Stan Liou
That is why I wrote "K is an (effectively constant) function". I didn't work out how $e$ varies, I just looked at the Peters & Matthews equation and estimated that the "leverage" of a (probably) very small e˙ was so small that it could be ignored.  Thanks for confirming it.  I am also assuming that masses don't change significantly (I think that is an important assumption of the author's model).
steveOw
<general-relativity><binary-star><pulsar>
Hulse Taylor Binary Pulsar - How is "Cumulative Periastron Time shift" calculated?
steveOw
# Hulse Taylor Binary Pulsar - How is "Cumulative Periastron Time shift" calculated?

(Asked by: username_0 on 2014-11-14)

Following on from a [previous question](https://astronomy.stackexchange.com/questions/6051/how-did-the-hulse-taylor-pulsar-provide-confirming-evidence-of-general-relativit) about the interesting [Hulse Taylor Binary Pulsar](https://en.wikipedia.org/wiki/Hulse%E2%80%93Taylor_binary).

**BACKGROUD**

[Weisberg & Taylor, 2004](http://arxiv.org/pdf/astro-ph/0407149v1.pdf) present a graph showing the change in "Cumulative Periastron Time Shift" over a 30 year period. Here is a later version of that graph
![orbit period decay](https://i.sstatic.net/992mC.gif)

The graph indicates excellent agreement between the observed data points and the theoretical curve predicted by the General Theory of Relativity (GTR).

The authors determined the theoretical change in orbit period $P$ over time (due to orbital energy loss from emission of gravitational waves) from a formula (Peters
& Matthews (1963)) of the form:-
$$ dP/dt = K/P^{(5/3)}$$
where K is an (effectively constant) function of the masses of the two stars and the orbit eccentricity. For the Hulse-Taylor Binary, the period $P$ is approximately 7.75 hours and the value of K is approximately $(-6 \* 10^{-5})$, giving an approximate value of $$ dP/dt = 2.34 × 10^{−12} s/s. $$

If $dP/dt$ is constant then obviously we can determine the period $P$ at a future time by applying the formula
$$P\_t = P\_0 + t\,(dP/dt).$$
As a rough initial estimate over one year the approximate change in period is: 365.25 x 24 x 3600 x -2.34 x 10^{-12} = 0.074 ms. (Wikipedia article indicates 0.076 ms.)

But $dP/dT$ is not constant because it depends on the value of $P$ which obviously changes. So we would expect the period to decay over time and the rate of decay to increase with time. However, the change in P is currently so very small that over the period of 30 years the change in dP/dt (which is a function of $P$) is extremely small also. (It is important to appreciate that the graph shown above does **not** plot $dP/dt$ against time).

**zibadawa timmy** (thanks) provided the formula $P=(8/3 \int K dt)^{3/8}$ giving:-

$$
P\_t = \left( \frac{8}{3} (Kt+C)\right)^{3/8} \qquad \mathrm{with}\qquad C=\frac{3}{8} P\_o^{8/3}
$$

hence

$$
P\_t =
\left( \frac{8}{3} Kt+ P\_o^{8/3})\right)^{3/8}
=
\left( \frac{-1.6\,t}{10^4}+ P\_o^{8/3}\right)^{3/8}
=
\left( \frac{-1.6\,t}{10^4}+ 27900^{8/3}\right)^{3/8}
$$
for 1 year , t = 31,557,600 s., so -1.6 t/10^4 = -5,049.12 = c. 5e3

for 10 years t = 310,557,600 s., so -1.6 t/10^4 = -50,491.2 = c. 5e4

while 27900^{8/3} = 716,051,966,959.181 = c. 7e11.

So
$$
P\_{1 year} = (-5049.12 + 716,051,966,959.181)^{3/8} = 27,899.999 926 s,\qquad deltaP\_{1 year} = -0.0000738 s.
$$
$$
P\_{10 years} = (-50491.2 + 716,051,966,959.181)^{3/8} = 27,899.999 262 s,\qquad deltaP\_{10 years} = -0.000738 s.
$$
These $deltaP$ values indicate that the value of $dP/dt$ is effectively constant ($=-0.000 073 8 s/s) over 10years.

---

**QUESTION 1 - What is plotted on the graph? ($\equiv$ What is Cumulative Perihelion Time Shift?) (ANSWERED).**

Initially I was confused about what the Weisberg & Taylor graph was actually plotting. The 2004 paper has "The observed and theoretical decays are compared graphically in Figure 1". The wikipedia article labels the y-axis "Cumulative Period Shift" but the caption reads "observed change in the epoch of periastron with date".

An earlier [1981 Scientific American article](http://www.nature.com/scientificamerican/journal/v245/n4/pdf/scientificamerican1081-74.pdf) indicates that
the vertical axis of the graph indicates the difference (round data points) between (A) the calculated epoch of a periastron event in a hypothetical "steady" system whose orbital period remains constant and (B) the observed epoch of the corresponding actual periastron event in the "decaying"-period system. Also shown is the continuous curve showing the difference between (A) steady system periastron epochs and (C) theoretical predictions of decayed periastron epochs. The periastron's of the actual and hypothetical orbits are synchronized at $t=0$ (in late 1974).

---

**QUESTION 2 - How are the theoretical values of CPTS Calculated? (ANSWERED).**

Several (similar) methods have now been identified.

Method 1 The method(s) used by Weisberg & Taylor to calculate theoretical values of CPTS have not been identified. These following methods seem reasonable candidates, based on the agreement with their results.

Method 2 (see Stan Liou's answer)
**Stan Liou** has provided a formula $\Delta t = 0.5(\dot P/P\_0)t^2$ that produces values which closely match the Weisberg & Taylor graph.

Method 3 (see my answer Method 3)
I have derived a formula identical to Stan Liou's but derived using phase instead of mean anomaly. Another formula is derived from the same trunk expressing $\Delta t$ in terms of $\dot P, P and N$ $\Delta t = -0.5 P\_o \dot P N^2$.

Method 4 (see my answer Method 4)
I have derived a formula derived using a coarsely approximated binomial series which gives similar results
$ \Delta t\_{N} = P\_0{\dot P} \frac{N^2+N}{2} $.

Method 5 (not shown)
A finely approximated series expansion approach has been assessed. For the given binary system it produces a formula that is effectively indistinguishable from that of the other methods for the given time range.

**Tags:** general-relativity, binary-star, pulsar
> **Comments:**
> **username_1**: As for solving the equation, $dP/dt=K/P^{5/3}$ is a separable differential equation with solution $$P=\left({8\over3}\int K\, dt \right)^{3/8},$$ with initial condition setting the constant of integration.  
> **username_0**: @zibadawa timmy. Many thanks. That appears to confirm my computer results. If possible please can you point me to the general method for your solution?  
> **username_1**: For the remainder of your issue, is $t$ measured in seconds or years? It sounds like seconds. But the data points are in years. 1975 to 1985 is a lot more than 10 seconds. In particular it looking linear for 20 seconds, or even a year, is no indication it is still linear an entire decade of seconds later.  
> **username_1**: I'm on a mobile, and the backspace is really close to the post button. So I hit the wrong one. I've edited it to finish the thought.  
> **username_0**: I have been using seconds.  Over 10 years 1975-1985 my (+your) formulae indicate relatively linear - from end Year 1 to end Year 10 deltaP increases by x10.  But the source graph indicates much stronger curve (you might not see it on a mobile!).  
> **username_1**: There are over 315 million seconds in a decade. Just eyeballing the equation, and assuming $t=0$ is 1975, tells me it's curving over that huge of an interval.  
> **username_0**: I have fleshed out the equations and deltaP derivations at the end of the Question.  
> **username_1**: You seem to be supposing that $P$ and the cumulative periastron time shift are the same. The periastron is the point of closest approach in an orbit (says Google). $P$ is the orbital period. These should be related (the "cumulative" part is especially important), but it seems unlikely that they are the same. So I think you need a different formula to get that graph.  
> **username_0**: Agreed. I'm assuming the graph should plot deltaP(t) vs. t, where deltaP(t) = P(t)-P(0) and t=0 in year 1975 AD. I cant presently think what "Cumulative Periastron time shift" should be if it is not same as deltaP.  
> **username_1**: Ultimately they are different things by their definition. $P$ is orbital period. Periastron is a point of closest approach in one orbit. The periastron time is how long into an orbit until you are at periastron (requires you pick a start point). And so the periastron time shift is how much the periastron time has changed between complete orbits (this does not actually depend on the  start time). Clearly a change in orbital period induces a periastron time shift, but it seems improbable that the shift equals the period change.  
> **username_1**: As for what the actual equation is, I do not know. One would hope it is mentioned in the paper. I am not familiar with this topic so I do not know the equation myself. I'm making educated guesses on the exact definitions of the terms, as it is.  
> **username_0**: Let us [continue this discussion in chat](http://chat.stackexchange.com/rooms/18669/discussion-between-steveow-and-zibadawa-timmy).  
> **username_2**: I've followed up on the PPN calculations for the $\dot{P} = KP^{-5/3}$ formula, and while your $K$ is a function of $e$, $e$ is itself not constant (should have checked sooner). But this makes very little difference because the rate of change of $e$ turns out to be quite negligible.  
> **username_0**: @Stan Liou That is why I wrote "K is an (effectively constant) function". I didn't work out how $e$ varies, I just looked at the Peters & Matthews equation and estimated that the "leverage" of a (probably) very small e˙ was so small that it could be ignored.  Thanks for confirming it.  I am also assuming that masses don't change significantly (I think that is an important assumption of the author's model).  

---

## 4 answers

### Answer by username_2 on 2014-11-15. Score: 2 | Accepted answer

Write the mean anomaly as a Taylor series ($n\equiv 2\pi/P$):
$$\begin{eqnarray\*}
M(t) \equiv \int\_0^tn\,\mathrm{d}t &=& M\_0 + \dot{M}\_0t + \frac{1}{2}\ddot{M}\_0t^2 + \mathcal{O}(t^3)\\
&=&n\_0t + \frac{1}{2}\dot{n}\_0t^2 + \mathcal{O}(t^3)\\
&=&\frac{2\pi}{P\_0}t - \pi\frac{\dot{P}\_0}{P\_0^2}t^2+\mathcal{O}(t^3)\text{.}
\end{eqnarray\*}$$

The periastron time (epoch) happens when $M(t) = 2\pi N$, where $N$ counts number of orbits, while $NP\_0$ *would* be the periastron time if the orbital period stayed constant. Therefore, treating $N$ as continuous, the cumulative difference is just
$$t - \frac{M(t)}{2\pi}P\_0 \approx \frac{1}{2}\frac{\dot{P}\_0}{P\_0}t^2\text{.}$$
If we want more accuracy, we can calculate the next term in the Taylor series,
$$\frac{1}{6}\ddot{n}\_0t^3 = \frac{\pi}{3P\_0^2}\left[\frac{2\dot{P}\_0^2}{P\_0}-\ddot{P}\_0\right]t^3\text{,}$$
and if $\dot{P} = KP^{-5/3}$ for $\dot{K}\approx 0$, we can eliminate the $\ddot{P}\_0$ term to give a cumulative shift of
$$t - \frac{M(t)}{2\pi}P\_0 \approx \frac{1}{2}\frac{\dot{P}\_0}{P\_0}t^2 - \frac{11}{18}\left(\frac{\dot{P}\_0}{P\_0}\right)^2t^3\text{.}$$
According to the Weisberg & Taylor paper, $|\dot{P}\_0/P\_0| \sim 10^{-16}\,\mathrm{s}^{-1}$, so the cubic term is negligible over the thirty-year time scale, and we can safely ignore anything past the quadratic term.

Just for kicks, we can solve $\dot{P} = KP^{-5/3}$ for constant $K$, which yields a cumulative periastron time shift of
$$\frac{3}{5x}\left[1+\frac{5}{3}xt-\left(1-\frac{8}{3}xt\right)^{5/8}\right]\text{.}$$
$K$ is not actually constant, being a function of the eccentricity $e$, although the rate of change of $e$ is tiny, on the order of $10^{-11}\,\mathrm{yr}^{-1}$, and we can ignore it on the scales we're talking about. Additionally, the proportionality formula is based on a orbital average of radiation reaction at 2.5-th post-Newtonian order.

---

The values presented in the paper are for the epoch $T\_0 = 52144.90097844$ (MJD):
$$\small{\begin{eqnarray\*}
P\_0 &=& 0.322997448930\,(4)\,\mathrm{d}\\
%\dot{P}\_0 &=& −(2.4184\pm0.0009)\times 10^{-12}\\
%\dot{P}\_0^{\scriptsize\text{corr}} &=& -(2.4056\pm 0.0051)\times 10^{-12}\\
\dot{P}\_0^{\scriptsize\text{GTR}} &=& −2.40242\,(2)\times 10^{-12}\\
\end{eqnarray\*}}$$
Thus, we can calculate the ratio $x\_0 = \dot{P}\_0/P\_0 = -8.60867\,(7)\times 10^{-17}\,\mathrm{s}^{-1}$. The difference between $T\_0$ and the beginning of 1975 is $\delta t = 9731\,\mathrm{d}$, which is not long enough for this ratio to change appreciably, since relative first-order correction would be around $|\dot{x}\_0(\delta t)/x\_0| \sim |x\_0(\delta t)| \sim 10^{-7}$.

Alternatively, following **zibadawa timmy**, we can also solve $K = \dot{P}P^{5/3}$ explicitly for constant $K$, which is not exactly true but approximates well:
$$\newcommand{\constKeq}{\overset{\scriptsize\dot{K}\approx 0}{=}}\begin{eqnarray\*}
P &\constKeq& \left[P\_0^{8/3}+\frac{8}{3}\dot{P}\_0P\_0^{5/3}t\right]^{3/8}\!\!\!\!\!\text{,}\\
x &\constKeq& \frac{\dot{P}}{P} = \frac{x\_0}{1+\frac{8}{3}x\_0t}\text{,}
\end{eqnarray\*}$$
and therefore $\left.x\right|\_{1975}$ is indeed indistinguishable from $x\_0$ at this level of accuracy.

I've graphed (in light blue) the result from the quadratic shift formula with this $x$ and superimposed the graph from the Weisberg & Taylor paper (cropped on outer borders and then rescaled to be $500\times 500\,\mathrm{px}$):

![Hulse–Taylor graph reconstruction](https://i.sstatic.net/jteR1.png)

It's an exact match except for a constant time shift, which is just a matter of when we set the $t = 0$ time. PSR B1913+16 was actually discovered in 1974, and the first data point is at the end of 1974. That is why my graph is shifted slightly by a fraction of a year.

---

> I am unclear whether your formula excludes the effect of periastron advance (4.226595 degrees per year):- (A) If not then the periastron advance effect is so large that the change in period due to gravitational wave emission is not visible on the graph and the agreement with theory cannot be assessed ...

The GTR prediction of the cumulative periastron time shift is made through through equation (1) in the paper, which prescribes $\dot{P}\_0 = \left.\dot{P}\_b\right|\_{t=0}$. The periastron time is well-defined because it depends on the closest approach, so the magnitude of the positional periastron advance $\langle\dot{\omega}\rangle = 4.226595^\circ/\mathrm{yr}$ is not directly relevant. GTR *also* predicts $\langle\dot{\omega}\rangle$, but that's not what we're considering here.

The mean anomaly is not actually an angle, despite counting off each periastron as $2\pi$. It isn't an angle even in exact Keplerian elliptic orbits, even though they have the exact same periastron at each period. Thus, $\langle\dot{\omega}\rangle$ is a different issue.

> (B) If so then must we conclude that gravitational wave emission causes change in period of approx 4s in 10 years?

The period changes by a bit less than $1\,\mathrm{ms}$ over a decade, as you have already calculated above.

> **Comments:**
> **username_0**: That gives a curve that looks close to the original. I haven't digested your method yet. I am unclear whether your method excludes the effect of periastron advance (4.226595 degrees per year):- (A) If not then it could be that the periastron advance effect (not on P but on the calculation, e.g. ambiguity of Mean Anomaly) is so large that any change in P due to gravitational wave emission (GWE) is not visible on the graph and so agreement with GWE theory cannot be assessed. (B) If so then must we conclude that GWE alone causes big changes in period e.g. approx 40s in 30 years?  
> **username_0**: OK I now get what "Cumulative Shift of Periastron Time" means (Update 2 in my Question). If $P$ varies with $t$ then periastron generally occurs when $t/P(t)$ is not an integer and so if $M(t) = 2\pi t/P(t)$ then $M(t)$ at periastron will not be an integral multiple of $2\pi$. (Or so it seems to me).  
> **username_2**: @steveOW $M(t) = \int_0^t\frac{2\pi}{P(t')}\,\mathrm{d}{t'}$ ensures that if the first periastron is when $M(t) = 0$, then the next is when $M(t) = 2\pi$, etc. If $P$ is constant, then $M(t) = 2\pi t/P$; otherwise, no.  
> **username_0**: I get $M''(0) =-4\pi P'(0)/P(0)^2 \neq d(M'(0))/dt$.  
> **username_2**: @steveOw $\dot{M} = 2\pi/P\,\implies\,\ddot{M} = -2\pi\dot{P}/P^2$. You have an extra factor of $2$.  
> **username_0**: I expect you are right but I will write my derivation out in a separate temporary answser to discover where I went wrong:-)  
> **username_0**: Agreed $\dot M(0) = 2\pi/P$ but $\dot M(t) = 2\pi(P-t \dot P)/P^2$ so $\ddot M(0)$ doesn't necessarilly $=d(\dot M(0))/dt$ (I think). Please see my separate answer regarding differentiation of $M$.  
> **username_0**: Thanks for another excellent answer. I satisfied myself with your solution using a Binomial formula independent of Mean Anomaly (see my other answer). I'm still not happy with definition of Mean Anomaly in a non-Keplerian orbit. I might ask a separate question about it.  
> **username_2**: @steveOw It's not assumption; it's just vocabulary. If $P$ were constant, then $1/P$ is orbital frequency; in time interval $\Delta t$, the number of (or fraction of) orbits completed is $\Delta t/P$. $P$ actually varies continuously, but that's what calculus is for: over infinitesimal $\mathrm{d}t$, the fraction of orbit traversed is $\mathrm{d}t/P$, and integration is adding them. Thus, $\int\mathrm{d}t/P$ is the *exactly correct* orbital count, regardless of whether one wishes to call its $2\pi$ scaling 'mean anomaly' or not. If prefer not, just replace 'mean anomaly' with 'orbital phase'.  
> **username_0**: OK I get it now $1/P_i$ is like speed around an imaginary phase clock dial which starts at 0 and ends ($360^o$ later) at 1. And then the $\Delta t/P_i$'s are like distances around the circumference of that clock dial and the sum of them all in one orbit period MUST add up to 1.  Thanks for your patience!  
> **username_0**: (A bit pedantic but) in deriving your formula for $\Delta t$ I think you make an implicit, unstated, unverified presumption that the times of the Nth periastra in steady-period (S) and decaying-period (D) systems are approximately equal i.e. $t_{NB}\approx t_{NA}$. For the system in question it turns out the assumption is valid as the error (as a fraction of $\Delta t$) is very small $~ 2\Delta t/t_{NB}$. My (Method 3) answer (based on yours) derives a formula for $\Delta t$ independent of $t$.  
> **username_2**: @steveOw No, it is based entirely on whether $M(t)$ is well-approximated, and by Taylor's theorem the error has about the magnitude of the the third-order term, which I've explicitly examined to be negligible. (To be fully rigorous, the remainder term needs to be examined explicitly, but it not essentially different.) Since cumulative periastron shift is $t - NP_0$ exactly and $N = M(t)/2\pi$ exactly, all I'm doing is plugging in the approximation for $M(t)$ to substitute for $N$. No other approximations are used, not even implicitly.  
> **username_0**: The Taylor approximation error isn't the issue. You use one variable $t$. I use two variables:- (for steady system A) $t_A=NP_0$  and (for decaying system B) $t_B$ as used in $M(t)$. OK now I see that your $t$ represents $t_B$ throughout, you don't actually write $t_A$.  Previously I was mis-reading your 4th equation as being of the expected form $t_A$-$t_B$ =$\Delta t$.  Sorry.  

---

### Answer by username_1 on 2014-11-14. Score: 2

As for solving the equation, $dP/dt=K/P^{5/3}$ is a separable differential equation. See the [first section here](http://en.wikipedia.org/wiki/Examples_of_differential_equations). Google also brings up several pages describing how to solve such equations.

In particular, we rewrite the equation in differential form as:
$$P^{5/3}\, dP = K\, dt.$$

Loosely speaking, I've cleared the denominators, with a goal to get the $P$ stuff all together on one side, and the $t$ stuff on the other. Technically it's not that at all, but pretending like it is works out. I am implicitly assuming that $K$ does not depend on $P$, which your post supports.

Now we integrate both sides:
$${3\over8}P^{8/3}=\int P^{5/3}\, dP = \int K\, dt,$$
where we let the right hand integral cover the constant of integration.

So our solution is $$P=\left({8\over3}\int K\, dt \right)^{3/8},$$ with initial condition setting the constant of integration.

Note I took an 8th root in finding $P$, so we need the inside function to be nonnegative.

Your post suggests $K$ is independent of time, in which case $\int K\, dt= K t+C$ for some constant $C$ which is determined by the initial conditions. If it does depend on time, you have to do the integration accordingly, or use data to estimate it.

> **Comments:**
> **username_0**: Many thanks :-) As you assumed: K is independent of time and P.  

---

### Answer by username_0 on 2014-11-24. Score: 0

**METHOD 3 Predicting Periastron Times and Cycle Durations from Orbital Phase ($\Phi$)**

In the answer by **Stan Liou** he uses a Taylor Series approximation of the [Mean Anomaly](http://en.wikipedia.org/wiki/Mean_anomaly) to derive a nice formula which determines the CPTS (Cumulative Perihelion Time Shift) value as a function of $t^2$. This formula produces results very close to those graphed by Weisberg & Taylor. As it took me a while to understand how the Mean Anomaly can be applied for a decaying-period orbit I thought it useful to note here what I learnt and present a slightly different way of obtaining a formula to predict periastron times.

Mean Anomaly basically indicates the phase of the orbit at a particular epoch, e.g. what time fraction of the orbit period has been completed. For Mean Anomaly the fraction is scaled by $2\pi$. Let $\Phi(t)$ be the phase of the orbit at time $t$ such that at the start of the orbit $t=0$ and $\Phi(0)$ = 0 and at the end of the orbit when $t=P$ (the period of orbit) $\Phi(P)$ = 1. I will assume that an orbit starts at one periapsis and finishes at the next periapsis.

In a decaying-period orbital system the value of $P(t)$ changes with $t$. Here $P$ is the period of orbit. In a decaying-period orbit I will assume that orbital phase $\Phi(t)$ is completely specified by time $t$ and either $P(t)$ or $F(t)$, i.e. aspidial precession, progressive changes in semi-major axis length or eccentricity do not lead to significant additional changes in phase.

In a decaying-period system the concept of period at a given epoch is somewhat abstract, it might be construed as a hypothetical period that would occur thereafter if the force causing period change ceased at that particular moment. Instead of period we can think in terms of $F(t)=1/P(t)$ where $F$ is the orbital frequency.

I like to think of $F(t)$ as "the rate of change of phase with time", i.e. $F(t) = \dot{\Phi}(t)$. Picture an imaginary phase clock comprising a circular dial whose circumference has markings running from $\Phi=0$ around in a full circle to $\Phi=1$. The phase $\Phi(t)$ of our subject system at any given epoch $t$ will be indicated by a marker at a particular point on the circumference of the dial. Then $F(t)$ can be thought of as the speed at which the marker is moving around the circumference of the dial at a given epoch $t$. In a particular short interval of time $\delta t$ the (system state) marker will move a certain distance around the dial and this distance travelled will indicate the change in phase. The "travel" (change of phase) will depend on the value of $F$ during that time as per $\delta \Phi \approx \delta t \* F(t)$. This is only approximate because F(t) changes during the interval $\delta t$.

Now let us assume that the time interval $T$ between succesive periaspes is known and that the time of the first periapsis is at $t=0$. If $P(t)$ (and therefore $F(t)$ too) changed in a step-wise manner between time intervals, then we could write

$$
\sum\_{i=1}^{i=N=T/\delta t} \frac{\delta t}{P(t)}
=
\sum\_{i=1}^{i=N=T/\delta t} \delta t F(t)
=
\sum\_{i=1}^{i=N=T/\delta t} \delta \Phi(t)
= 1.
$$

To represent continuously-changing $P$ and $F$ we reduce $\Delta t$ towards zero and obtain these integral functions

$$
\int\_{t=0}^{t=T} \frac{1}{P(t)}\, \mathrm{d}t
=
\int\_{t=0}^{t=T} F(t)\, \mathrm{d}t
=
\int\_{t=0}^{t=T} \dot{\Phi}(t)\, \mathrm{d}t
= 1.
$$

At any time $\tau$ during the orbit ($0\leq \tau \leq 1$) the current instantaneous phase is given by
$$
\Phi(\tau) =
\int\_{t=0}^{t=\tau} \dot{\Phi}(t)\, \mathrm{d}t
=
\int\_{t=0}^{t=\tau} F(t)\, \mathrm{d}t
=
\int\_{t=0}^{t=\tau} \frac{1}{P(t)}\, \mathrm{d}t.
$$

and given a constant value of $\dot P$ (the rate of change of period), and $P\_o=P(0)$ we obtain
$$
\Phi(\tau)
=
\int\_{t=0}^{t=\tau} \frac{1}{P\_o + \dot P t}\, \mathrm{d}t
\qquad \mathrm{and} \qquad
\dot{\Phi}(\tau) = F(\tau) = \frac{1}{P(\tau)} = \frac{1}{P\_o + \dot P \tau}.
$$

Using

$$
\frac{\mathrm{d}}{\mathrm{d}t} \left(
\frac{1}{f(x)}
\right)
=
\frac{-\dot f(x)}{f(x)^2}
\qquad \mathrm{and} \qquad
f(x) = P\_o + \tau \dot P
\qquad \mathrm{and} \qquad
\dot f(x) = \dot P
$$

we get
$$
\ddot \Phi(\tau)
=
\frac{\mathrm{d}}{\mathrm{d}t} \left(
\frac{1}{P\_o + \tau \dot P}
\right)
=
\frac{-\dot P}{P\_o^2 + (\tau \dot P)^2 + 2 P\_o \tau \dot P}
$$

and for $\tau = 0$ we get
$$
\Phi(0) = 0
\qquad \mathrm{and} \qquad
\dot{\Phi}(0) = \frac{1}{P\_o }
\qquad \mathrm{and} \qquad
\ddot \Phi(0)
= \frac{-\dot P}{P\_o^2}.
$$

---

**PREDICTING PERIASTRON TIMES**

We can obtain an expression for CPTS by copying Stan Liou's approach , but using phase ($\Phi$) instead of Mean Anomaly ($M$).

Write the phase at time $\tau$ as a Maclaurin series (a Taylor series with $a=0$):
$$\begin{eqnarray\*}
\Phi(\tau) \equiv \int\_{t=0}^{t=\tau} \dot{\Phi}(t)\,\mathrm{d}t
&=& \Phi(0) + \dot{\Phi}(0) \tau + \frac{1}{2}\ddot{\Phi}(0) \tau^2 + \mathcal{O}(\tau^3)\\
&=& 0 +\frac{1}{P\_0}\tau - \frac{\dot{P}}{2P\_o^2}\tau^2+\mathcal{O}(\tau^3)\text{.}
\end{eqnarray\*}$$

Dropping the third order and residual error terms (Stan Liou's answer explains justification for this) we get

$$
\Phi(\tau) \approx \frac{1}{P\_0}\tau- \frac{\dot{P}}{2P\_o^2}\tau^2.
$$

Let $N$ count the number of completed orbits at some time $\tau\_N$. Then the phase at the end of orbit number $N$ is given by $\Phi(\tau\_N) = N$. In a steady system (A) the orbital period $P\_o$ stays constant $\dot P=0$ and so the time $\tau\_{NA}$ at which the $N$th orbit completed is given exactly by $\tau\_{NA} = NP\_o$ and thus $N=\tau\_{NA}/P\_o$. In a decaying-period system (B) the Nth orbit completes at time $\tau\_{NB}$ when $\Phi(\tau\_{NB}) = N$ but with $\Phi(\tau\_{NB})$ depending on a non-zero value of $\dot P$.

$$
\begin{align}
&\Phi(\tau\_{NA}) = N = \frac{1}{P\_0}\tau\_{NA}\\
&\Phi(\tau\_{NB}) = N \approx \frac{1}{P\_0}\tau\_{NB}- \frac{\dot{P}}{2P\_o^2}\tau\_{NB}^2.
\end{align}
$$

We can obtain a formula for $\Delta t\_{N}$ which is the difference in time between $\tau\_{NB}$ (the epoch of the decaying system $N$th periastron) and $\tau\_{NA}$ (the epoch of the steady system $N$th periastron) as follows. Note that following Weisberg & Taylor we define $\Delta t\_{N}\,(=$ CPTS $)=\tau\_{NB} - \tau\_{NA}$. The $N$th periastron in the decaying period system occurs earlier than the $N$th periastron in the steady system. Therefore $\Delta t\_{N}$ (=$ CPTS $) will become increasingly negative as time passes (as $t$ increases).
$$
\frac{\tau\_{NA}}{P\_o} = N \approx \frac{\tau\_{NB}}{P\_0}- \frac{\dot P}{2P\_o^2} \,{\tau\_{NB}}^2
$$

$$
{\tau\_{NA}} \approx {\tau\_{NB}} - \frac{\dot P}{2P\_o} \,{\tau\_{NB}}^2
$$

$$
\Delta t\_N
= \tau\_{NB} - \tau\_{NA} \approx + \frac{\dot P}{2P\_o} \,{\tau\_{NB}}^2
$$

which is Stan Liou's approximation for $\Delta t\_N$.

So what value is this equation?
In general assume that we have determined the time $t\_0$ of the $0$th periastron and measured accurately an initial orbital period $P\_o$ and we keep track of the observed periastra.

Case 1 - we observe that the $N$th periastron occurs at a particular time ($\tau\_{NB}$). Using $\tau\_{NA}=N P\_o$, we can easily calculate $\Delta t\_N = \tau\_{NB} - \tau\_{NA}$. Then, using Stan Liou's approximation of $\Delta t\_N$ we can obtain an empirical estimate of $\dot P$ from
$$
{\dot P} = \frac{2 P\_o \,\Delta t\_N} {\tau\_{NB}}^2.
$$

Case 2 - we have a theoretical formula which, given our measurement of $P\_o$ predicts the value of $\dot P$. Given various values of $N$ and using $\tau\_{NA}=N P\_o$ we can easily calculate the hypothetical epochs ($\tau\_{NA}$ for various $N$) of the hypothetical steady-system periastra. Now we wish to plot a curve showing the theoretical values of $\Delta t\_N$ at various hypothetical periastron epochs ($\tau\_{NA}$). But Stan Liou's equation for $\Delta t\_N$ uses $\tau\_{NB}$ not $\tau\_{NA}$. We need to find $\Delta t\_N$ as a function of $\tau\_{NA}$. Therefore we would need to solve the following quadratic equation for $\Delta t\_N$ in terms of $\tau\_{NA}$ and $K$, where $K=\frac{\dot P}{2 P\_o}$:

$$
\Delta t\_{N} = K {\tau\_{NB}}^2 = K ({\tau\_{NA}}^2 + 2 \tau\_{NA}\Delta t\_{N} +{\Delta t\_{N}}^2)
$$

$$
0 = (K{\tau\_{NA}}^2 + 2K \tau\_{NA}\Delta t\_{N} +K{\Delta t\_{N}}^2) -\Delta t\_{N}
$$

$$
0 = (K){\Delta t\_{N}}^2 +(2K \tau\_{NA} -1)\Delta t\_{N} +(K{\tau\_{NA}}^2)
.$$

---

The following way of getting an equation for $\Delta t$ is more protracted but it will provide a formula for $\Delta t$ as a function of $N$ (which can easily be converted into a function of $\tau\_{NA}$).

Starting with
$$
\frac{\tau\_{NA}}{P\_o} =\frac{1}{P\_0}\tau\_{NB}- \frac{\dot P}{2P\_o^2} \,{\tau\_{NB}}^2
$$

$$
\tau\_{NA} =\tau\_{NB}- \frac{\dot P}{2P\_o}{\tau\_{NB}}^2
\qquad\rightarrow \qquad
0 = -\tau\_{NA} + \tau\_{NB} - \frac{\dot P}{2P\_o} \,{\tau\_{NB}}^2
$$

Using $ x = \frac{-b \pm \sqrt{b^2-4ac}}{2a}$ with $a=- \frac{\dot P}{2P\_o} , b=1, c=-\tau\_{NA}$ we obtain
$$
\tau\_{NB} = \frac{-(1) \pm \sqrt{(1)^2-4(- \frac{\dot P}{2P\_o}).(-\tau\_{NA})}}{2(- \frac{\dot P}{2P\_o})}
$$

$$
\tau\_{NB} = \frac
{1 \pm \sqrt {1 - 2 \frac{ \dot P }{P\_o}.\tau\_{NA} }}
{\dot P /{P\_o}}
\qquad \longrightarrow \qquad
\tau\_{NB} =
\frac{P\_o}{ \dot P } \left( 1 \pm \sqrt {1 - 2 \frac{\dot P }{P\_o}.\tau\_{NA} }\,\, \right)
$$

Now $\tau\_{NA} = P\_o \,N$ so we can write
$$
\tau\_{NB} =
\frac{P\_o}{ \dot P } \left(1 \pm \sqrt {1 - 2 \dot P N } \right)
\qquad \longrightarrow \qquad
-\Delta t = \tau\_{NA}-\tau\_{NB}
= P\_o \left( N -\frac{1}{ \dot P} \left(1 \pm \sqrt {1 - 2 \dot P N } \right) \right)
$$

= = = = = = =

We can try a simple approximation of this using $(1+\alpha)^{0.5} = (1+ 0.5\alpha)$ so that $\sqrt {1 - 2 \dot P N } \approx 1- \dot P N$ and by inspection the "$\pm$" becomes a "$+$"
$$
-\Delta t
= P\_o \left( N -\frac{1}{ \dot P } \left(1 \pm 1 - \dot P N \right) \right)
= P\_o ( N - N )
= 0
$$

This is a valid approximation but it goes a bit too far for our purposes! We know that $\tau\_{NA}-\tau\_{NB} \neq 0$. So let us try expanding the term $[1 -2 \dot P N]^{0.5}$ in a series expansion,

$$
[1 + (-2 \dot P N)]^{0.5}
= 1 +\frac{1}{2}(-2\dot P N)^{1} -\frac{1}{8}(-2 \dot P N)^{2}
+\frac{1}{16}(-2\dot P N)^{3}+...
$$

$$
[1+ (-2\dot P N)]^{0.5}
= 1 - \dot P N -\frac{1}{2}(\dot P N)^{2}
-\frac{1}{2}( \dot P N)^{3}+...
$$

Let us limit ourselves to terms with $ \dot P $ to the power of $2$ or less thus

$$
[1-2 \dot P N]^{0.5} \approx
1 - \dot P N -\frac{1}{2}\dot P ^2 N^2 .
$$

We can insert this into the time difference equation

$$
-\Delta t
\approx P\_o \left( N -\frac{1}{ \dot P }
\left(1 \pm (1 - \dot P N -\frac{1}{2} \dot P ^2N^2 \right) \right)
$$

by inspection the "$\pm$" becomes a "$-$"

$$-\Delta t\approx P\_o \left( N -\frac{1}{ \dot P }
\left( + \dot P N +\frac{1}{2}\dot P ^2 N^2 \right) \right)
\qquad \longrightarrow \qquad
-\Delta t\approx P\_o \left( N - \left( N +\frac{1}{2} \dot P N^2) \right) \right)
$$

$$
-\Delta t \approx -\frac{1}{2}P\_o \dot P \, N^2
\qquad \longrightarrow \qquad
\Delta t \approx \frac{1}{2}P\_o \dot P \, N^2
$$

---

From the Weisberg & Taylor paper we are given $P\_o=27906.979587552 s$ and $\dot P = -2.40242 \*10^{-12} s/s$
so $\frac{1}{2}P\_o \dot P = -0.5 \* 27906.979587552 \* 2.40242 \*10^{-12} = -3.352216747 \* 10^{-08} s$.

With $N= 10,000$ cycles we get $\Delta t = -3.352217 s$.

Note that we can make the substitution ${\tau\_{NA}}^2 = N^2 Po^2$ into the equation for $\Delta t$ to give
$$
\Delta t \approx \frac{1}{2}P\_o \dot P \, N^2 = \frac{1}{2} \frac {\dot P}{P\_o} \,{\tau\_{NA}}^2
.$$
The difference between this $\Delta t$ and the value of $\Delta t$ from Stan Liou's equation is
$$
\frac{1}{2} \frac {\dot P}{P\_o} \left( {\tau\_{NA}}^2 -{\tau\_{NB}}^2 \right)
$$
$$
=
\frac{1}{2} \frac {\dot P}{P\_o} \left(
{{\tau\_{NA}}^2
-\tau\_{NA}^2 +2 \tau\_{NA} \Delta t + \Delta t^2} \right)
$$
$$
=
\frac{1}{2} \frac {\dot P}{P\_o} \left(
{2 \tau\_{NA} \Delta t + \Delta t^2} \right)
$$
and the difference as a fraction of $\Delta t$
$$
=
\frac{1}{2} \frac {\dot P}{P\_o} \left(
{2 \tau\_{NA} + \Delta t} \right)
.$$

For the present case $(\dot P/P\_o \approx 8.6 \* 10^{-17})$ and so, after 34,000 periastron cycles (948,838,000 s, about 30 years) the difference between the two approximations of $\Delta t$ would be only $0.8 \* 10^{-8} s$ which is much smaller than the time resolution of the observations.

---

**PREDICTING CYCLE DURATIONS**

We can also derive an expression for the cycle duration $D$, as follows.

The start periaston is at $t=0=T0$ and the first following periaston is at $t=T1$ then
$$
\Phi(T1)=\int\_{t=0}^{t=T1}\frac{1}{P(t)}\,\mathrm{d}t = 1
$$
this can be expressed as
$$
1 =\int\_{t=0}^{t=T1}\frac{1}{ \dot{P} t + P\_0 }\,\mathrm{d}t
= \left[ \frac{1}{\dot{P}} ln ( C\dot{P}t + CP\_0 )\right]\_0^{T1}
$$
hence
$$
\dot{P} =
ln ( C\dot{P}T1 + CP\_0 ) - ln ( CP\_0 )
=
ln \left( \frac{ C ( \dot{P}T1 + P\_0 )} { C ( P\_0 )} \right)
=
ln \left( \frac{ \dot{P}T1 + P\_0 } {P\_0 } \right)
$$

$$
\rightarrow exp(\dot{P})
= \frac{ \dot{P}T1 + P\_0 } {P\_0 }
\rightarrow
P\_0 exp(\dot{P}) -P\_0
= \dot{P}T1
$$
giving us
$$
\rightarrow
T1 =
P\_0 \frac{exp(\dot{P}) -1} {\dot{P}}
$$

so the duration of the first orbit (from start periastron to first following periastron) is $D\_1 = T1-T0 = T1$ thus
$$
D\_1 =
D\_0 \frac{exp(\dot{P}) -1} {\dot{P}}
$$

and we can generalize this to
$$
D\_{N} =
D\_{N-1} \left( \frac{exp(\dot{P}) -1} {\dot{P}} \right)
.$$

For example using the made-up value $\dot P = -2.34\,10^{-8}$ we obtain
$$
D\_{N+1} =
D\_{N} \frac{exp(-2.34\*10^{-5}) -1} {-2.34\*10^{-8}}
= 0.999 999 988 D\_{N}
$$

However, using standard arithmetical software (such as Excel) when we try to calculate using $\dot P = -2.34\,10^{-10}$ we get nonsense results because of truncation errors. The approximate value of $\dot P$ reported for the Hulse-Taylor system is about $-2.34\,10^{-12}$.

We can analyze the formula for $D\_{N+1}$ using series expansions. A Taylor Series expansion of $exp(x)/x -1/x$ is given by WolframAlpha as

$$
exp(x)/x -1/x= 1 +
\frac{x}{2}
+ \frac{x^2}{6}
+ \frac{x^3}{24}
+ \frac{x^4}{120}
+ \frac{x^5}{720}
+ \frac{x^6}{5040}
+ \frac{x^7}{40320}
+ \frac{x^8}{362880}
+ \frac{x^9}{3628800}
+ \frac{x^{10}}{39916800}
+O(x^{11})
$$
or
$$1 + \frac{x^{2-1}}{2} + \frac{x^{3-1}}{3\*2} + ... \frac{x^{i-1}}{i!} + ...$$
No easilly computable expression or approximation is obvious at present.
Proceeding anyway, it is clear that the duration of any subsequent orbit $N$ can be computed from

$$
D\_N
=
D\_0 \left(\frac{exp^{\dot P} -1} {\dot P} \right) ^N
=
D\_0 \left(\frac{1}{\dot P} \right)^N \left(exp^{\dot P} -1 \right) ^N
$$

It is interesting to compare this expression for cycle duration with that which was used as the basis of the coarse binomial solution (see other answer: Method 4)

$$D\_1=D\_0(1+\dot P)^1$$

The ratio of durations Phase-based to Coarse binomial based is
$$
D\_0 \left( \frac{exp(\dot P)-1}{\dot P} \right)^{1}
:D\_0(1+\dot P)^1
$$

becoming
$$
\frac{exp(\dot P)-1}{\dot P}
:1+\dot P
\qquad \rightarrow
exp(\dot P)-1
:\dot P +\dot P^2
\qquad \rightarrow
exp(\dot P)
: 1+\dot P +\dot P^2
$$

applying the series expansion of $exp(\dot P)$ we obtain
$$
\qquad \rightarrow
1+\dot P + \frac{\dot P^2}{2} + \frac{\dot P^3}{6} + \frac{\dot P^4}{24} +\,...\,
: 1+\dot P +\dot P^2
$$

Clearly the two expressions give different values for cycle duration $D$ and the difference appears at the term in $\dot P^2$ and seems to be no bigger than $\frac{\dot P^2}{2}$.

> **Comments:**
> **username_3**: I do not think you should do this way  
> **username_0**: @Py-ser. Please could you elaborate? (Note there have been major edits since your comment).  
> **username_3**: I meant, to edit more answers instead of one. I did not criticize your actual answer, but the way of answering.  
> **username_0**: @Py-ser. Thanks for the feedback. It was certainly messy yesterday.  Hopefully now it is a bit more simple. I have changed  the question to signpost the different methods.  I am reluctant to merge my two answers into one because they are different methods with different terms employed (even though they produce very similar results).  

---

### Answer by username_0 on 2014-11-20. Score: 0

**METHOD 4:Predicting CPTS using a Coarse Binomial Series approximation.**

This method does not make any assumptions about Mean Anomaly. Orbit is defined as the cycle from one periapsis to the next periapsis.

---

**SUMMARY**

The time from the $0$th to the $N$th periastron can be obtained from the following sum, from which a binomial expression can be derived

$$\frac{T\_N}{P\_0}
\qquad = \sum\_{i=1}^{i=N} (1+\dot P)^i
\qquad = N + \sum\_{j=2}^{N+1}\binom{N+1}{j}\dot P^{j-1}.$$

For our purposes terms in $\dot P$ greater than first order can be dropped leading to

$$ T\_{N} = N{P\_0} + P\_0 \frac{N^2+N}{2} {\dot P} .$$

The time difference between the $N$th periapses of a steady period system and a decaying period system is then

$$ \Delta t\_{N} = 0.5 P\_0 {\dot P} (N^2+N) $$

This gives results which agree closely with the Wesiberg & Taylor graph and the formula in Stan Liou's answer.

---

**Definitions and Assumptions**

Let us consider the variation over time in the **period** of orbit. Weisberg & Taylor provide a formula for the variation in "period" with time. However the meaning of "period at a point in time" is ambiguous when the period is changing with time. We therefore have to define how "period at a point in time" translates into a time interval between the start and end of a particular orbital cycle (e.g. from periapsis to periapsis). In this method we will "fudge" the translation by assuming (A1) that the time interval between successive periapses is equal to the "hypothetical period at the point in time of the end of the cycle". Also to "boot" the model we will assume (A2) that at time zero ($t=0$) the hypothetical period is $H\_o$ and thus the hypothetical period at the end of the first cycle is given by $H\_1 = H\_o + H\_o \dot H$. Therefore, by (A1) we have the following definition of the duration of the first cycle $D\_1 = H\_1 = H\_o + H\_o \dot H$.

We further assume (A3) that, over the period of consideration (~ 40 years) the rate of change of period with time is effectively constant ($\ddot H = 0$ and $\dot H(t) = constant = \dot H$ ). (See Stan Liou's answer for justification). Thus we have
$$
H(t) = H(0) + t \dot H .
$$
It is emphasized that these assumptions are not strictly mappable to a consistent physical model. However given that we are dealing with very small changes in period over a very large number of orbits the ruleset may be useful. We shall proceed and see what comes out. Henceforth I will not distinguish between hypothetical period $H$ and cycle duration $D$. I will assume that they are both represented by the same term period $P$ and that the duration $P$ of a particular periapsis to periapsis cycle is equal to the hypothetical period $P(t)$ where $t$ is the epoch of the end of the cycle. Also $\dot P$ will represent both the rate of change of hypothetical period and also the rate of change of cycle duration between cycles.

$P\_N$ is the period of the $N$th orbit (orbit $O\_N$). Orbit number $0$ terminates at $t=T\_0=0$ at perihelion number $0$ and has period $Po = P(0)$. The first considered orbit ($O\_1$) starts at time $t=0$ at perihelion number $0$ and finishes at perihelion number $1$ at time $t=T\_1$ and has period duration $P\_1 = P\_0+ P\_0\dot P =P(T\_1)$. Therefore we can write

$$P\_1 = P\_0(1+\dot P)^1$$
$$P\_2 = P\_1(1+\dot P) = P\_0(1+\dot P)(1+\dot P)= P\_0(1+\dot P)^2$$
$$P\_3 = P\_2(1+\dot P) = P\_0(1+\dot P)(1+\dot P)(1+\dot P)= P\_0(1+\dot P)^3$$
In general then, the period of any orbit is given by the formula:
$$P\_N = P\_0(1+\dot P)^N.$$

Now let us consider the **time** $T$ at which each orbit ends. We have
$$T\_1 = T\_0+P\_1 = P\_1 = P\_0(1+\dot P)$$
$$T\_2=T\_1+P\_2 = P\_1+P\_2 = P\_0(1+\dot P)+P\_0(1+\dot P)^2$$
$$T\_3=T\_2+P\_3 = P1+P2+P3 = P\_0(1+\dot P)+P\_0(1+\dot P)^2+P\_0(1+\dot P)^3$$
$$T\_N=P\_1+P\_2+...+P\_N = P\_0(1+\dot P)+P\_0(1+\dot P)^2+...+P\_0(1+\dot P)^N.$$

Thus
$$T\_N = {P\_0}\sum\_{i=1}^{i=N} (1+\dot P)^i\, .$$

(The Note at the end of this answer explores a way of proceeding which uses an exact equivalent of this expression in a non-serial form. However, that route is not productive due to difficulties with exponentials.)

---

We shall proceed as follows, (using $k$ in place of $\dot P$)
$$ \frac{T\_N}{P\_0}= S\_N = \sum\_1^N (1+k)^i $$

We can multiply out the first few terms and examine the patterns in the following **pyramid**...
$$(1+k)^1 = k +1 $$
$$(1+k)^2 = k^2 +2k +1 $$
$$(1+k)^3 = k^3 +3k^2 +3k +1 $$
$$(1+k)^4 =k^4+4 k^3+6 k^2+4 k+1$$
$$(1+k)^5 =k^5+5 k^4+10 k^3+10 k^2+5 k+1$$
$$(1+k)^6 =k^6+6 k^5+15 k^4+20 k^3+15 k^2+6 k+1$$

So for example, for $N=3$ we would obtain the sum
$$
S\_3= k^3 +4k^2 +6k +1
$$

The results suggest a solution with a pattern of the form
$$
S\_N = a + bk^1 +ck^2+dk^3...
$$

We can see straightaway that $a=N$. The other coefficients increase with $N$. It might be possible to determine a formula for the coefficients from the pattern. Although the general pattern is not convergent, it is possible that for certain restricted ranges of $N$ and $k$ a convergent formula could be obtained. If so then it might be possible to obtain a useful approximation of $S\_N$ by evaluating the first few terms in the series.

---

The following binomial equation was derived in an answer to [this question](https://math.stackexchange.com/questions/1029172/how-do-i-express-the-sum-1k1k21k3-1kn-for-k1-as-a-se) at Maths Stack Exchange by User73985

$$ \sum\_1^N (1+k)^i \qquad = \qquad N + \sum\_{j=2}^{N+1}\binom{N+1}{j}k^{j-1}$$
So
$$ S\_N = N + \sum\_{j=2}^{N+1}\frac{(N+1)!}{(N+1-j)!j!}k^{j-1}$$

then
$$ S\_N = N
+ \frac{(N+1)!}{(N-1)!2!} k^{1}
+ \frac{(N+1)!}{(N-2)!3!} k^{2}
+ \frac{(N+1)!}{(N-3)!4!} k^{3} +...
$$

giving
$$ S\_N = N
+ \frac{(N+1)(N)}{2!} k^{1}
+ \frac{(N+1)(N)(N-1)}{3!} k^{2}
+ \frac{(N+1)(N)(N-1)(N-2)}{4!} k^{3} +...
$$

thus
$$ S\_N= N
+ \frac{N^2+N}{2} k^{1}
+ \frac{N^3-N}{6} k^{2}
+ \frac{N^4-2 N^3-N^2+2 N}{24} k^{3} +...
$$
(Switching back from using $k$ to using $\dot P$...)
For $N= 10,000$ and ${k=\dot P}= 2.40242 \* 10^{-12}$ the ratio of each term to the next term is greater than $1 : N{\dot P} \approx 1:2.4 \* 10^{-8}$ so for our purposes this formula can be approximated by truncating to the first order power of ${\dot P}$ to give
$$ S\_{N} = \frac{T\_{N}}{P\_0} = N + \frac{N^2+N}{2} {\dot P}^{1} $$
So
$$ T\_{N} = N{P\_0} + P\_0 \frac{N^2+N}{2} {\dot P} $$
and so the time difference between the Nth periapses (fixed $P$ - decaying $P$) is

$$ \Delta t\_{N} = P\_0{\dot P} \frac{N^2+N}{2} $$
and for $N=10,000$ we obtain (from Excel) $\Delta t = 3.35255 s$ which is very close to the value $(3.35221 s)$ obtained from Stan Liou's formula ($\Delta t= 0.5 {\dot P}/P\_0 t^2$).

Over 30 years (33,924.42 cycles) $\Delta t$ = 38.5806 $s$ which is very close to the amount indicated by the graph and by Stan Liou's formula $(38.5794 s)$.

---

**Note**

WolframAlpha indicates the equality

$$T\_N = {P\_0}\sum\_{i=1}^{i=N} (1+\dot P)^i
\qquad =
\frac{P\_0}{\dot P}
\,
(1+\dot P)[(1+\dot P)^N -1]
$$

For small $\dot P$ this can be approximated by:-
$$T\_N =\frac{P\_0}{\dot P}\,(1+\dot P)[(1+N\dot P) -1]
=\frac{P\_0}{\dot P}\,(1+\dot P)(N\dot P)
=N{P\_0}\,(1+\dot P).
$$
The end time of $N$ cycles of period $P\_{constant}=P\_0$ will be $T\_{NC}=N P\_0$. And so the time difference between the end times of (i) N cycles of constant $P$ and (ii) $N$ cycles of decaying $P$ is approximated by
$$\Delta t = NP\_0 - NP\_0(1-\dot P) = NP\_0 \dot P.
$$
For $P\_0$ = 0.322997448930 d = 27906.979587552 s and $\dot P$ = -2.40242 x 10$^{-12}$ we obtain $(P\_0\,dot P) =$ -6.70442859007267$^{-08}$. For $N=10,000$ cycles (approximately 10 years) this indicates a time difference of 6.70442859007267 $^{-03}$ s, about $7 ms$. This time difference is very small compared to the value roughly indicated by the original graph and the value calculated by Stan Liou's formula, from which we would expect a time difference of about 3$s$. The explanantion is that the first-order approximation is not suitable in this case because the true value of $\Delta t$ is so small compared to $t$.


---

