PSR1913+16: The Binary Pulsar, Mathematical Reproduction and Extension

Kathleen Hamilton, Howard Community College 2017, University of Maryland College Park
Ye Hie (Joy) Cho, Howard Community College

Mentored by: Alex M. Barr, Ph.D.


The discovery of PSR1913+16, a binary pulsar system, brought light to the first indirect detection of gravitational waves in the early 1970s. This work, conducted by Joseph Taylor’s research team, which included graduate student Russell Hulse, led to the 1993 Nobel Prize in Physics. The research presented here aims to reproduce some of their calculations, then extend beyond those to analyze the eventual collision. The initial prediction for coalescence is given, along with the strength of the gravitational waves produced by the system. Through calculus and astrophysics concepts, the mathematical assumptions and equations underlying Hulse and Taylor’s work are verified and compared with their data.



“Pulsars … yielded two exciting scientific stories which began with serendipity and ended with a Nobel Prize.” – Russell Hulse, Nobel Prize acceptance speech [1]

Since the stunning direct confirmation of the existence of gravitational waves in September 2015 by the Laser Interferometer Gravitational Wave Observatory (LIGO), gravitational waves have been at the forefront of scientific research and discovery. Gravitational waves can best be visualized as ripples in space-time, produced by accelerating objects. Because gravitational waves are not part of the electromagnetic spectrum, they propagate through the universe with exceptionally little degradation and provide information about the vast cosmic phenomena. Due to the extreme distances between us and the gravitational wave sources, the waves are much less powerful when they reach our detectors, which is part of the reason why it took one hundred years after Einstein’s prediction to directly detect them.

The first indirect detection of gravitational waves was made decades earlier, thanks to another remarkable scientific discovery. In 1967, Jocelyn Bell Burnell discovered pulsars by recording incredibly regular radio pulses [Fig. 1]. Further investigation revealed the source to be neutron stars, whose rapid rotation rates create the regular electromagnetic signals. Neutron stars are astonishingly dense remains of larger, gaseous stars that have collapsed and undergone supernova explosions. There is a possibility that all neutron stars are in fact pulsars, just with their electromagnetic jets aligned in such a way that we cannot observe them from our position on Earth. Pulsars are regarded as universal clocks due to their incredible regularity and can even be used to make cosmic maps because of their precise, measurable periods of rotation.


Figure 1: Diagram of a pulsar, showing its radio wave emissions and regular period. Reproduced from [1].

       In 1974, graduate student Russell Hulse, in working with Joseph Taylor’s pulsar research group, discovered an intriguing phenomenon: a pulsar – PSR1913+16 – with a very fast period that changed noticeably over the course of a day! Even better, the shifts were reproducible, meaning that PSR1913+16’s period was a periodic function of time [2]. Hulse suspected that he had found a binary pulsar, with the pulses shifted due to the pulsar and its companion star orbiting around their center of mass. Due to the Doppler effect, the time between pulses decreases when the pulsar approaches Earth and increases when it moves away. As Hulse continued to observe, he found the evidence he was looking for, and the measurements to detect gravitational waves began [Fig. 2].


Figure 2: Diagram of the pulsar (left) and its companion, showing the system releasing gravitational waves as they orbital. Looking at the picture from the perspective of someone on Earth, the pulsar period would appear shortened when moving toward Earth (on the left), then elongated when moving away (on the right). Reproduced from [3].


A basic tenet of general relativity is that any accelerating mass radiates gravitational waves. In a binary system, energy will be lost due to this radiation, causing the orbital radius to shrink over time. To simplify initial calculations, a circular system is considered in determining the luminosity of gravitational wave emissions [2]. Luminosity (L) is the energy emitted by a star per second. As shown in the equation below, luminosity depends on the total mass of the system M, the sum of the masses of the pulsar and its companion; the separation between the two bodies a, twice the orbital radius r; the orbital angular velocity ; the gravitational constant G; and the speed of light c.


LaTeX: L_{cir}=\frac{2GM^2a^4\omega_{orb}^6}{5c^5}\tag{1}

The luminosity should equal the negative time derivative of the total energy of the system. Summing the kinetic energy (KE) and gravitational potential energy (GPE) yields the total energy of the system. We assume here that both bodies have the same mass so that LaTeX: M=m_1+m_2=2m.

LaTeX: E=KE+GPE=2\left(\frac{1}{2}m\omega_{orb}^2r^2\right)-\frac{Gm^2}{2r}\tag{2}

Newton’s second law states that the net force on a system equals its mass times its acceleration. This is combined with his law of universal gravitation, which states that the product of two masses divided by the square of their separation is proportional to their attractive force.

LaTeX: \frac{Gm^2}{4r^2}=m\omega_{orb}^2r\tag{3}

This allows us to rewrite Equation (2) as

LaTeX: E=\frac{Gm^2r-2Gm^2r}{4r^2}=-\frac{Gm^2}{4r}=-\frac{GM^2}{8a}\tag{4}

As previously mentioned, the orbital separation shrinks in conjunction with the system’s energy loss. Equating the time derivative of the system’s total energy with the negative of luminosity gives an equation for the change in the pulsar and its companion’s separation over time.

LaTeX: \frac{dE}{dt}=\frac{GM^2}{8a^2}\frac{da}{dt}=-L_{cir}
LaTeX: \frac{da}{dt}=-\frac{L_{cir}8a^2}{GM^2}=-\frac{\frac{2GM^2a^4\omega_{orb}^68a^2}{5c^5}}{GM^2}=-\frac{16a^6\omega_{orb}^6}{5c^5}=-\frac{16a^6}{5c^5}\left(\frac{GM}{a^3}\right)^3
LaTeX: \frac{da}{dt}=-\frac{16G^3M^3}{5a^3c^5}\tag{5}

The orbital period is inversely proportional to the orbital angular velocity.

LaTeX: T=\frac{2\pi}{\omega_{orb}}

The fractional rate of change of the period can be written in terms of orbital velocity. Further equivalences can be derived from the above expressions involving Newton’s second law (3) and the rate of change in the orbital separation (4).

LaTeX: \frac{1}{T}\frac{dT}{dt}=-\frac{1}{\omega_{orb}}\frac{d\omega_{orb}}{dt}=-\frac{3}{2a}\frac{da}{dt}=-\frac{24G^3M^3}{5a^4c^5}\tag{6}

Having derived the basic equations, the system’s eccentricity can now be considered. Represented by , the eccentricity accounts for the elliptical orbits of the pulsar and its companion. For circular orbits, the value of eccentricity is 0, and it is between 0 and 1 for elliptical orbits. The eccentricity function multiplies the above equations, making them more accurate for non-circular motion [2]. The previous equations have not explicitly included this function, since an eccentricity of 0 for circular orbits means the function’s value is 1.

LaTeX: f(\epsilon)=(1-\epsilon^2)^{-7/2}\left(1+\frac{73}{24}\epsilon^2+\frac{37}{96}\epsilon^4\right)

Hulse and Taylor found that the pulsar and its companion orbit each other elliptically with an eccentricity value of 0.617. Inputting this value into the eccentricity function gives the modification factor.

LaTeX: f(0.617)=(1-0.617^2)^{-7/2}\left(1+\frac{73}{24}0.617^2+\frac{37}{96}0.617^4\right)\approx11.843

With this, the system’s yearly orbital decay can be found, for both the separation and period.

LaTeX: \frac{da}{dt}=-\frac{16G^3M^3}{5a^3c^5}f(\epsilon)
LaTeX: \frac{da}{dt}=-\frac{16(6.674\times10^{-11}m^3s^{-2}kg^{-1})^3(5.62887\times10^{30}kg)^3}{5(2\times9.762\times10^8m)^3(3\times10^8 ms^{-1})^5}(11.843)(31536000\frac{s}{yr})
LaTeX: \frac{da}{dt}\approx -3.5\frac{m}{yr}


LaTeX: \frac{dT}{dt}=\frac{3T}{2a}\frac{da}{dt}=\frac{3(2.79\times10^4s)}{2(2\times9.762\times10^8m)}\left(-3.5\frac{m}{yr}\right)\approx -80\frac{\mu s}{yr}\tag{7}

Using the known masses for the pulsar and its companion and their current separation, we find that their separation shrinks by approximately 3.5 meters each year and the period of their orbit decreases by approximately 80 microseconds each year.

Kepler’s third law of planetary motion states that the square of the orbital period of a planet is directly proportional to the cube of its orbital semi-major axis. This can be applied to a binary system by taking both of the masses into account [4]. Using Kepler’s third law, the current orbital period can be determined.

LaTeX: T_0=\sqrt{\frac{a^3}{m_1+m_2}}=\sqrt{\frac{(2\times0.006525AU)^3}{1.42M_\oplus +1.41M_\oplus }}=8.8\times10^{-4}yr=2.79\times10^4s

PSR1913+16’s most significant result, its periastron shift, can now be calculated. The periastron is the point in the pulsar’s orbit where it is closest to its companion star. The periastron shift refers to how many seconds ahead of schedule the pulsar arrives at the periastron after several orbits compared to when it would arrive if the orbit was not shrinking due to energy lost through gravitational waves. Let LaTeX: \Delta represent the decrease in orbital period, which can be converted to a pure, dimensionless number.

LaTeX: \Delta\approx 80\frac{\mu s}{yr}\times\frac{1yr}{31536000s}\approx2.3823\times^{-12}

Let T0 represent the initial period, originally observed by Hulse and Taylor in 1975. Since the period is decaying, the next orbital period takes slightly less time than the original.

LaTeX: T_1=T_0(1-\Delta)

This pattern continues with each subsequent period being a fraction LaTeX: (1-\Delta) of the previous period.

LaTeX: T_2=T_1(1-\Delta)=T_0(1-\Delta)(1-\Delta)=T_0(1-\Delta)^2
LaTeX: T_i=T_0(1-\Delta)^i

Since the orbital decay, LaTeX: \Delta, is quite tiny, this expression can be approximated with the first term in a Taylor series expansion.

LaTeX: T_i\approx T_0(1-i\Delta)

This expression allows for the computation of the total time required for the binary system to complete many orbits, which is the accumulated periastron arrival time. Subtracting the product of the number of orbits and the initial period isolates the arrival shifts, which are the decrease in time due to orbital decay.

LaTeX: accumulated~shift=T_0\sum_{i=1}^n(1-i\Delta)-nT_0\approx -T_0\Delta\frac{n(n+1)}{2}\tag{8}

Here n represents the number of orbits. The accumulated shifts in arrival time over the course of several years can be plotted, showing the orbital decay due to energy lost through gravitational radiation [Fig. 3].


Figure 3: Comparison between the observed accumulated shift in periastron arrival time and the predicted shift based on Equation (8). Experimental values are taken from the 2010 report by Weisberg et al. [5]. The excellent agreement between theory and observations provided conclusive – albeit indirect – proof of gravitational waves and led to the 1993 Nobel prize for Hulse and Taylor.

       After studying the orbital decay, the natural question became when the pulsar and its companion would coalesce. To solve for the collision time, the orbital radius decay equation (5) was set up as a differential equation.

LaTeX: \frac{da}{dt}=-\frac{16G^3M^3}{5a^3c^5}f(\epsilon)=-\frac{16\left(6.674\times10^{-11}m^3s^{-2}kg^{-1}\right)^3\left(5.62887\times10^{30}kg\right)^3}{5a^3\left(3\times10^8ms^{-1}\right)^5}(11.843)
LaTeX: \frac{da}{dt}\approx -\frac{8.26837876\times10^{20}m^4s^{-1}}{a^3}

After separating the variables, both sides can be integrated.

LaTeX: \int_{2\times9.762\times10^8}^0a^3da=\int_0^{collide}-8.26837876\times10^{20}dt
LaTeX: -8.26837876\times10^{20}collide=-\frac{\left(2\times9.762\times10^8\right)^4}{4}
LaTeX: collide=\frac{-3.63258014\times10^{36}m^4}{-8.26837876\times10^{20}m^4s^{-1}}=4.393340274\times10^{15}s\times\frac{1yr}{31536000s}
LaTeX: collide\approx 140~million~years

Gravitational waves from PSR1913+16 have already been indirectly detected, which leads to the question of direct detection. How does the strain, h, of its gravitational waves – the fractional stretching and compressing of space-time – compare to recent events detected at LIGO?

LaTeX: h\approx 4\frac{GM_c}{c^2D_L}\left(\frac{G\pi fM_c}{c^3}\right)^{2/3}
LaTeX: M_c=\frac{(m_1m_2)^{3/5}}{(m_1+m_2)^{1/5}}=\frac{\left(\left(1.42\cdot1.989\times10^{30}kg\right)\left(1.41\cdot1.989\times10^{30}kg\right)\right)^{3/5}}{\left((1.42+1.41)\cdot1.989\times10^{30}kg\right)^{1/5}}
LaTeX: h\approx \frac{\left(4\cdot6.674\times10^{-11}m^3s^{-2}kg^{-1}\right) M_c}{\left(3\times10^8ms^{-1}\right)^2\left(1.98675\times10^{20}m\right)}\left(\frac{\left(6.674\times10^{-11}m^3s^{-2}kg^{-1}\right)\cdot\pi\cdot2\cdot M_c}{\left(3\times10^8ms^{-1}\right)^3\left(2.79\times10^4s\right)}\right)^{2/3}
LaTeX: h\approx 4.4965\times10^{-23}

Such a minute measure of strain is still beyond our current capabilities for direct detection, which are on the order of 10-21 [6]. Further explanation of the strain equation is out of the scope of this paper. For a more detailed discussion of strain (h), chirp mass (LaTeX: M_c), and how distance (LaTeX: D_L) see the paper by Tzortzakakis and Marvi in this journal.

Going Further

This research provides several avenues for continued exploration. A major future goal is to produce simulations of PSR1913+16, showing its orbital decay and eventual merger. Additionally, PSR1913+16’s eccentricity function can be more refined. While the calculations matched well with the observed data, the eccentricity is not likely to stay constant as the pulsar and its companion move closer together. Even though the estimate for the collision time seems to be on the correct order of magnitude, a more dynamic eccentricity function will produce an improved figure, consistent with other estimates [6].



The authors gratefully acknowledge the support and resources of the Howard Community College Science Division generally, and the Gravitational Waves Honors Seminar of Fall 2017 specifically.


Contacts: and


[1] Hulse, R. A. (1993 December 8). The Discovery of the Binary Pulsar. Nobel lecture. Retrieved from

[2] Amato, J. C. (2016 April 19). Gravitational radiation 1. Physics from Planet Earth. Retrieved from

[3] The Nobel Prize in Physics 1993 – Illustrated Information. Nobel Media AB 2014. Retrieved from

[4] Astronomy Department at Cornell University. Binary star motions. Retrieved from

[5] Weisberg J. M., Nice D. J., and Taylor J. H. (2010). Timing measurements of the relativistic binary pulsar PSRB1913+16. The Astrophysical Journal, 722, 1030-1034.

[6] Wheeler J. (2013). Gravitational waves. Retrieved from


Icon for the Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License

Journal of Research in Progress Vol. 1 by Howard Community College is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, except where otherwise noted.

Share This Book