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

We are searching data for your request:

Forums and discussions:
Manuals and reference books:
Data from registers:
Wait the end of the search in all databases.
Upon completion, a link will appear to access the found materials.

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) ight)^{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}) ight)^{3/8} = left( frac{-1.6,t}{10^4}+ P_o^{8/3} ight)^{3/8} = left( frac{-1.6,t}{10^4}+ 27900^{8/3} ight)^{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. Write the mean anomaly as a Taylor series ($nequiv 2pi/P$): $$egin{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{2pi}{P_0}t - pifrac{dot{P}_0}{P_0^2}t^2+mathcal{O}(t^3) ext{.} end{eqnarray*}$$ The periastron time (epoch) happens when$M(t) = 2pi 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)}{2pi}P_0 approx frac{1}{2}frac{dot{P}_0}{P_0}t^2 ext{.}$$ 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{2dot{P}_0^2}{P_0}-ddot{P}_0 ight]t^3 ext{,}$$ 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)}{2pi}P_0 approx frac{1}{2}frac{dot{P}_0}{P_0}t^2 - frac{11}{18}left(frac{dot{P}_0}{P_0} ight)^2t^3 ext{.}$$ 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 ight)^{5/8} ight] ext{.}$$$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{egin{eqnarray*} P_0 &=& 0.322997448930,(4),mathrm{d} %dot{P}_0 &=& −(2.4184pm0.0009) imes 10^{-12} %dot{P}_0^{scriptsize ext{corr}} &=& -(2.4056pm 0.0051) imes 10^{-12} dot{P}_0^{scriptsize ext{GTR}} &=& −2.40242,(2) imes 10^{-12} end{eqnarray*}}$$ Thus, we can calculate the ratio$x_0 = dot{P}_0/P_0 = -8.60867,(7) imes 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: $$ewcommand{constKeq}{overset{scriptsizedot{K}approx 0}{=}}egin{eqnarray*} P &constKeq& left[P_0^{8/3}+frac{8}{3}dot{P}_0P_0^{5/3}t ight]^{3/8}!!!!! ext{,} x &constKeq& frac{dot{P}}{P} = frac{x_0}{1+frac{8}{3}x_0t} ext{,} end{eqnarray*}$$ and therefore$left.x ight|_{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 imes 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 ight|_{t=0}$. The periastron time is well-defined because it depends on the closest approach, so the magnitude of the positional periastron advance$langledot{omega} angle = 4.226595^circ/mathrm{yr}$is not directly relevant. GTR also predicts$langledot{omega} angle$, but that's not what we're considering here. The mean anomaly is not actually an angle, despite counting off each periastron as$2pi$. It isn't an angle even in exact Keplerian elliptic orbits, even though they have the exact same periastron at each period. Thus,$langledot{omega} angle$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. 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: $${3over8}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({8over3}int K, dt ight)^{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. 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}inom{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_0dot 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 This Site by User73985 $$sum_1^N (1+k)^i qquad = qquad N + sum_{j=2}^{N+1}inom{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+Ndot P) -1] =frac{P_0}{dot P},(1+dot P)(Ndot 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$2pi$. 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$ au$during the orbit ($0leq au leq 1$) the current instantaneous phase is given by $$Phi( au) = int_{t=0}^{t= au} dot{Phi}(t), mathrm{d}t = int_{t=0}^{t= au} F(t), mathrm{d}t = int_{t=0}^{t= au} 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( au) = int_{t=0}^{t= au} frac{1}{P_o + dot P t}, mathrm{d}t qquad mathrm{and} qquad dot{Phi}( au) = F( au) = frac{1}{P( au)} = frac{1}{P_o + dot P au}.$$ Using $$frac{mathrm{d}}{mathrm{d}t} left( frac{1}{f(x)} ight) = frac{-dot f(x)}{f(x)^2} qquad mathrm{and} qquad f(x) = P_o + au dot P qquad mathrm{and} qquad dot f(x) = dot P$$ we get $$ddot Phi( au) = frac{mathrm{d}}{mathrm{d}t} left( frac{1}{P_o + au dot P} ight) = frac{-dot P}{P_o^2 + ( au dot P)^2 + 2 P_o au dot P}$$ and for$ au = 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$ au$as a Maclaurin series (a Taylor series with$a=0$): $$egin{eqnarray*} Phi( au) equiv int_{t=0}^{t= au} dot{Phi}(t),mathrm{d}t &=& Phi(0) + dot{Phi}(0) au + frac{1}{2}ddot{Phi}(0) au^2 + mathcal{O}( au^3) &=& 0 +frac{1}{P_0} au - frac{dot{P}}{2P_o^2} au^2+mathcal{O}( au^3) ext{.} end{eqnarray*}$$ Dropping the third order and residual error terms (Stan Liou's answer explains justification for this) we get $$Phi( au) approx frac{1}{P_0} au- frac{dot{P}}{2P_o^2} au^2.$$ Let$N$count the number of completed orbits at some time$ au_N$. Then the phase at the end of orbit number$N$is given by$Phi( au_N) = N$. In a steady system (A) the orbital period$P_o$stays constant$dot P=0$and so the time$ au_{NA}$at which the$N$th orbit completed is given exactly by$ au_{NA} = NP_o$and thus$N= au_{NA}/P_o$. In a decaying-period system (B) the Nth orbit completes at time$ au_{NB}$when$Phi( au_{NB}) = N$but with$Phi( au_{NB})$depending on a non-zero value of$dot P. egin{align} &Phi( au_{NA}) = N = frac{1}{P_0} au_{NA} &Phi( au_{NB}) = N approx frac{1}{P_0} au_{NB}- frac{dot{P}}{2P_o^2} au_{NB}^2. end{align} We can obtain a formula forDelta t_{N}$which is the difference in time between$ au_{NB}$(the epoch of the decaying system$N$th periastron) and$ au_{NA}$(the epoch of the steady system$N$th periastron) as follows. Note that following Weisberg & Taylor we define$Delta t_{N},(=$CPTS$)= au_{NB} - au_{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{ au_{NA}}{P_o} = N approx frac{ au_{NB}}{P_0}- frac{dot P}{2P_o^2} ,{ au_{NB}}^2$$ $${ au_{NA}} approx { au_{NB}} - frac{dot P}{2P_o} ,{ au_{NB}}^2$$ $$Delta t_N = au_{NB} - au_{NA} approx + frac{dot P}{2P_o} ,{ au_{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 ($ au_{NB}$). Using$ au_{NA}=N P_o$, we can easily calculate$Delta t_N = au_{NB} - au_{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} { au_{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$ au_{NA}=N P_o$we can easily calculate the hypothetical epochs ($ au_{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 ($ au_{NA}$). But Stan Liou's equation for$Delta t_N$uses$ au_{NB}$not$ au_{NA}$. We need to find$Delta t_N$as a function of$ au_{NA}$. Therefore we would need to solve the following quadratic equation for$Delta t_N$in terms of$ au_{NA}$and$K$, where$K=frac{dot P}{2 P_o}$: $$Delta t_{N} = K { au_{NB}}^2 = K ({ au_{NA}}^2 + 2 au_{NA}Delta t_{N} +{Delta t_{N}}^2)$$ $$0 = (K{ au_{NA}}^2 + 2K au_{NA}Delta t_{N} +K{Delta t_{N}}^2) -Delta t_{N}$$ $$0 = (K){Delta t_{N}}^2 +(2K au_{NA} -1)Delta t_{N} +(K{ au_{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$ au_{NA}$). Starting with $$frac{ au_{NA}}{P_o} =frac{1}{P_0} au_{NB}- frac{dot P}{2P_o^2} ,{ au_{NB}}^2$$ $$au_{NA} = au_{NB}- frac{dot P}{2P_o}{ au_{NB}}^2 qquad ightarrow qquad 0 = - au_{NA} + au_{NB} - frac{dot P}{2P_o} ,{ au_{NB}}^2$$ Using$ x = frac{-b pm sqrt{b^2-4ac}}{2a}$with$a=- frac{dot P}{2P_o} , b=1, c=- au_{NA}$we obtain $$au_{NB} = frac{-(1) pm sqrt{(1)^2-4(- frac{dot P}{2P_o}).(- au_{NA})}}{2(- frac{dot P}{2P_o})}$$ $$au_{NB} = frac {1 pm sqrt {1 - 2 frac{ dot P }{P_o}. au_{NA} }} {dot P /{P_o}} qquad longrightarrow qquad au_{NB} = frac{P_o}{ dot P } left( 1 pm sqrt {1 - 2 frac{dot P }{P_o}. au_{NA} },, ight)$$ Now$ au_{NA} = P_o ,N$so we can write $$au_{NB} = frac{P_o}{ dot P } left(1 pm sqrt {1 - 2 dot P N } ight) qquad longrightarrow qquad -Delta t = au_{NA}- au_{NB} = P_o left( N -frac{1}{ dot P} left(1 pm sqrt {1 - 2 dot P N } ight) ight)$$ = = = = = = = We can try a simple approximation of this using$(1+alpha)^{0.5} = (1+ 0.5alpha)$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 ight) ight) = P_o ( N - N ) = 0$$ This is a valid approximation but it goes a bit too far for our purposes! We know that$ au_{NA}- au_{NB} eq 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}(-2dot P N)^{1} -frac{1}{8}(-2 dot P N)^{2} +frac{1}{16}(-2dot P N)^{3}+…$$ $$[1+ (-2dot 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 ight) ight)$$ by inspection the "$pm$" becomes a "$-$" $$-Delta tapprox P_o left( N -frac{1}{ dot P } left( + dot P N +frac{1}{2}dot P ^2 N^2 ight) ight) qquad longrightarrow qquad -Delta tapprox P_o left( N - left( N +frac{1}{2} dot P N^2) ight) ight)$$ $$-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${ au_{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} ,{ au_{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( { au_{NA}}^2 -{ au_{NB}}^2 ight)$$ $$= frac{1}{2} frac {dot P}{P_o} left( {{ au_{NA}}^2 - au_{NA}^2 +2 au_{NA} Delta t + Delta t^2} ight)$$ $$= frac{1}{2} frac {dot P}{P_o} left( {2 au_{NA} Delta t + Delta t^2} ight)$$ and the difference as a fraction of$Delta t$$$= frac{1}{2} frac {dot P}{P_o} left( {2 au_{NA} + Delta t} ight) .$$ 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 ( Cdot{P}t + CP_0 ) ight]_0^{T1}$$ hence $$dot{P} = ln ( Cdot{P}T1 + CP_0 ) - ln ( CP_0 ) = ln left( frac{ C ( dot{P}T1 + P_0 )} { C ( P_0 )} ight) = ln left( frac{ dot{P}T1 + P_0 } {P_0 } ight)$$ $$ightarrow exp(dot{P}) = frac{ dot{P}T1 + P_0 } {P_0 } ightarrow P_0 exp(dot{P}) -P_0 = dot{P}T1$$ giving us $$ightarrow 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}} ight) .$$ 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} ight) ^N = D_0 left(frac{1}{dot P} ight)^N left(exp^{dot P} -1 ight) ^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} ight)^{1} :D_0(1+dot P)^1$$ becoming $$frac{exp(dot P)-1}{dot P} :1+dot P qquad ightarrow exp(dot P)-1 :dot P +dot P^2 qquad ightarrow exp(dot P) : 1+dot P +dot P^2$$ applying the series expansion of$exp(dot P)$we obtain $$qquad ightarrow 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}\$.