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. Many thanks. That appears to confirm my computer results. If possible please can you point me to the general method for your solution? 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. 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. 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!). 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. I have fleshed out the equations and deltaP derivations at the end of the Question. 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. 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. 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. 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. 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 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). • Stan Liou Correct answer 8 years ago 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? 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$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. I get$M''(0) =-4\pi P'(0)/P(0)^2 \neq d(M'(0))/dt$. @steveOw$\dot{M} = 2\pi/P\,\implies\,\ddot{M} = -2\pi\dot{P}/P^2$. You have an extra factor of$2$. I expect you are right but I will write my derivation out in a separate temporary answser to discover where I went wrong:-) 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$. 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 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'. 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! (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 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. 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.