# Hulse Taylor Binary Pulsar - How is “Cumulative Periastron Time shift” calculated?

Following on from a previous question about the interesting Hulse Taylor Binary Pulsar.

BACKGROUD

Weisberg & Taylor, 2004 present a graph showing the change in "Cumulative Periastron Time Shift" over a 30 year period. Here is a later version of that graph

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 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. • 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 Nov 14 '14 at 19:51 • @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 Nov 14 '14 at 21:05 • 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 Nov 14 '14 at 21:50 • 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 Nov 14 '14 at 21:54 • 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 Nov 14 '14 at 22:09 ## 4 Answers 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}$): 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. • 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 Nov 15 '14 at 13:56 • 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 Nov 15 '14 at 19:22 • @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 Nov 15 '14 at 19:32 • I get$M''(0) =-4\pi P'(0)/P(0)^2 \neq d(M'(0))/dt$. – steveOw Nov 16 '14 at 0:52 • @steveOw$\dot{M} = 2\pi/P\,\implies\,\ddot{M} = -2\pi\dot{P}/P^2$. You have an extra factor of$2$. – Stan Liou Nov 16 '14 at 4:15 As for solving the equation,$dP/dt=K/P^{5/3}$is a separable differential equation. See the first section here. 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. • Many thanks :-) As you assumed: K is independent of time and P. – steveOw Nov 14 '14 at 21:46 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 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$. 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 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}\$.

• I do not think you should do this way – Py-ser Nov 24 '14 at 1:55
• @Py-ser. Please could you elaborate? (Note there have been major edits since your comment). – steveOw Nov 26 '14 at 13:28
• I meant, to edit more answers instead of one. I did not criticize your actual answer, but the way of answering. – Py-ser Nov 26 '14 at 13:40
• @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 Nov 26 '14 at 17:00