A Complete Analytic Gravitational Wave Model for Undergraduates
AA Complete Analytic Gravitational Wave Model forUndergraduates
Dillon Buskirk a and Maria C. Babiuc Hamilton b Department of Physics, Marshall University, Huntington, WV 25755, USA
Abstract
Gravitational waves are produced by orbiting massive binary objects, such as black holes andneutron stars, and propagate as ripples in the very fabric of spacetime. As the waves carry off orbitalenergy, the two bodies spiral into each other and eventually merge. They are described by Einstein’sequations of General Relativity. For the early phase of the orbit, called the inspiral, Einsteinequations can be linearized and solved through analytical approximations, while for the late phase,near the merger, we need to solve the fully nonlinear Einstein’s equations on supercomputers. Inorder to recover the gravitational wave for the entire evolution of the binary, a match is requiredbetween the inspiral and the merger waveforms. Our objectives are to establish a streamlinedmatching method, that will allow an analytical calculation of the complete gravitational waveform,while developing a gravitational wave modeling tutorial for undergraduate physics students. Weuse post-Newtonian (PN) theory for the inspiral phase, which offers an excellent training ground forstudents, and rely on
Mathematica for our calculations, a tool easily accessible to undergraduates.For the merger phase we bypass Einstein’s equations by using a simple analytic toy model namedthe Implicit Rotating Source (IRS). After building the inspiral and merger waveforms, we constructour matching method and validate it by comparing our results with the waveforms for the firstdetection, GW150914, available as open-source. Several future projects can be developed basedfrom this project: building complete waveforms for all the detected signals, extending the post-Newtonian model to take into account non-zero eccentricity, employing and testing a more realisticanalytic model for the merger, building a separate model for the ringdown, and optimizing thematching technique. a Undergraduate Student, email: [email protected] b email: [email protected] a r X i v : . [ g r- q c ] O c t . INTRODUCTION Why are gravitational waves important? They ripple through the very fabric of space-timeand, like other kind of waves, carry information about their sources: cataclysms of cosmicproportion, such as colliding black holes, exploding supernovae, and even the origins of theuniverse. This information, once decoded, will enable us to answer deep and fundamentalquestions about the Universe and the nature of space and time.The first direct detection of gravitational waves happened on September 14, 2015, thanksto the precise instruments of the Laser Interferometer Gravitational-wave Observatory(LIGO). This discovery, known as the GW150914 event , came from the collision of twoblack holes. The announcement came at the beginning of 2016, almost as if to celebrate the100 th anniversary of Einstein’s Theory of General Relativity. Three scientists who playedan instrumental role in this discovery received the 2017 Nobel prize in physics for theircontribution to one of the most important achievements in the history of science.Several more detections followed suit: GW151226 , GW170104 , GW170608 , GW170814 ,and GW170817 . The last event is also known as the golden binary , because the collisiongave off – besides gravitational waves – electromagnetic radiation across the spectrum, andhundreds of Earth masses of precious and heavy elements.With these recent discoveries, the era of gravitational and multi-messenger wave astron-omy has begun, and crucial to the success of this new science is the development of a reliablepipeline of well-prepared and capable researchers, ready to move this field forward.The Einstein field equations of general relativity are essential for the correct modelingof gravitational waves. These complicated partial differential equations are extremely chal-lenging to solve analytically. For example, in order to obtain the correct gravitational wavesignal from two orbiting stars, the solution would have to contain everything, starting withthe birth of the binary from the interstellar gas, continuing with the evolution to the merger,and finishing with its collapse into a black hole, with an emission of gravitational radiation .So far, no exact solution exists for such complex systems. In fact, only a few exact solutionsof Einstein’s equations are known . Therefore, in order to calculate the gravitational waves,we need to resort either to numerical simulations or to analytical approximations .This research started as senior undergraduate capstone project during the academic yearof 2017-2018, when the analytical models were implemented, and applied to the first 5 grav-2tational wave detections, without being tested. During the summer of 2018, we continuedwith work on the matching technique, in order to obtain complete analytical gravitationalwaveforms for the entire evolution of the binary, and validate it with the template for theGW150914. We will give here an extended report of this work and the results obtained.We start by presenting in Sec. II, the two analytical approaches we are using to calculatethe gravitational waves: the post-Newtonian (PN), and the Implicit Rotating Source (IRS)models, explaining their domains of applicability. In Sec. III, we elaborate the algorithmemployed in numerically determining the evolution variables that enter in the calculationof the gravitational wave amplitude, using
Mathematica , and we calculate the strain fora fiducial waveform, for both the inspiral and the merger. In Sec. IV, we concentrate ourefforts in finding a simple and effective method for matching the end of the inspiral with thebeginning of the merger, and we build a complete analytic waveform for the whole evolutionof a fiducial binary configuration. Furthermore, we investigate whether the model is correctby calculating the gravitational waveform for the binary configuration corresponding tothe GW150914 event, and testing it against the open-source template for the GW150914strain . We obtain a very good overlap between our calculated strain and the template,which proves that our calculations are correct. Finally, in Sec. V, we summarize our workand give a short outline of future undergraduate projects that can be developed.Throughout this work we are using geometric units, such that G = c = 1. Therefore, theequation for the Newtonian gravitational field simplifies to: Φ = − M/r . Those units arecommonly used in general relativity calculations and undergraduates are likely unfamiliarwith them. It is useful to review the relationship between those units and the
InternationalSystem of Units (SI), because it it is not intuitive and its usefulness is not easy to grasp. Wedescribe in Appendix A how to convert between SI and geometric units, giving examples onhow this allows us to measure mass in seconds, distance in mass, etc. Among the symbolsused in the paper are the following: the total mass of the binary M = m + m , with m and m as the individual masses, the symmetric mass ratio η = m m /M , the orbital separationbetween the location of the centers of the stars in the binary r , the overall orbital velocity v , the orbital angular velocity ω , the frequency of the gravitational waves f GW , the phaseΦ, the distance from the detector to the binary R , the polarization modes of the strain( h × , h + ), and the amplitude of the wave A . Appendix B and Appendix C contain the PNand gIRS coefficients employed in our calculations. In Appendix D we present a short step3y step tutorial that if followed, will enable the reader to change the necessary parametersin order to obtain complete gravitational waveforms for different binary configurations. II. THE ANALYTICAL MODELS
Analytical models of compact binary inspiral and merger are commonly used in combi-nation with numerical simulations, in order successfully to build a large bank of accuratewaveform templates . The best-known analytical procedure of calculating gravitationalwaves uses the PN theory, which was originally developed by Einstein in order to findsolutions to his field equations of general relativity. This theory shows how to construct per-turbative solutions of Einstein’s equations as a series of successive approximations in powersof v/c , in the case of slow motion, large separation, and weak gravitational fields . The ratioof the velocity of the source to the speed of light is called the post-Newtonian parameter, x = v /c . Although its validity is limited to weak fields, the post-Newtonian theory offersremarkably accurate predictions for the gravitational radiation emitted by compact binarysystems . We will show below how to use the post-Newtonian approach to calculate thegravitational waves emitted during the inspiral of two black holes. We must keep in mindthat this approach breaks down for large speeds, as x gets close to unity, and it cannot beapplied at the merger, where we have to rely either on numerical relativity or on analyticaltoy models. One of the simpler and well known analytical ansatz used to model the mergercase is the Implicit Rotating Source model, as described in . We will use the generic
IRS toy-model for the merger , which is tuned to numerical relativity, it is easy to use,and gives satisfactory results. Other merger models are presented in .Once each phase is modeled, a match is required between the post-Newtonian waveformdescribing the inspiral and the merger waveform obtained with the gIRS toy-model. Wedevise a simple and efficient stratagem for the matching region between the PN generatedwaveform and the gIRS wave model, very close to the merger, and we prove its validity. A. The Inspiral Model
We know that in Newton’s law of universal gravitation, a binary system is stable, and willorbit indefinitely with constant frequency, without emitting gravitational waves, or shrinking4heir orbit. However, Einstein’s general theory of relativity predicts that the two orbitingstars will gradually lose orbital energy through emission of gravitational waves and will comecloser together until they merge to form one single star of larger mass, or a black hole. Inthis situation, the energy will decrease in time, and the equation for the energy balancedescribing this behavior will be written as: dEdt = −F , (1)where E is the energy of the binary and F is the flux of the emitted gravitational waves.Even the Earth–Sun system emits gravitational waves, but their effect is too weak (10 − loss relative to the total orbital energy) to be taken into account.In the post-Newtonian approximation, the equations of general relativity take the form ofthe familiar Newtonian two–body equations of motion, in the limit v/c →
0, called the weakfield limit. A correction of order ( v/c ) n to the Newtonian equation of motion is counted asan n/ dvdt = − G Mr [1 + 1 P Nc + 1 . P Nc + 2 P Nc + 2 . P Nc + ... ] (2)At each post-Newtonian expansion we unravel new physics beyond the Newtonian realm.For example, the 1 st order recovers orbit precession, the 1 . th order describes spin-orbitinteraction, and the 2 nd order spin-spin coupling dynamics. The orbital decay with emissionof gravitational waves appears from the 2 . th order onward. The current state of the art inthe post–Newtonian expansion is 3 . th order .Using Kepler’s third law of planetary motion ( ω r = GM ) and writing the orbitalvelocity as v = ωr = ( GM ω ) / , we obtain the following important relationship between thepost-Newtonian parameter x and the orbital angular velocity ω : x = v c = ( GM ω ) / c (3)In geometrical units, eq. (3) becomes x = ( M ω ) / . Using the chain rule, we rewrite theenergy balance eq. (1) in terms of x as: dxdt = − F dE/dx . (4)There are several well known ways of solving eq. (4), referred to as the Taylor T1 throughT5 approximants . We will use the Taylor T4 approximant method, which was reported5o give better agreement with numerical relativity than other approximants for binarieswith comparable mass. This method expands (4) in post-Newtonian powers of x : dxdt = dx dt x + dx dt x + dx dt x + dx dt x + dx HT dt . (5)Here ˙ x HT stands for hereditary terms, and represent the higher order post-Newtonian non-linear terms, which depend on the dynamics of the system in its entire past, the so called tails and tails-of-tails terms that account for the nonlinear interaction between the gravitationalwaves and the spacetime itself. Each post-Newtonian term in eq. (5) is further expressed asa power series in x , truncated at the appropriate order. We will use up to the 6 P N -orderterms, which are the highest order calculated for a quasi-circular orbit . The Taylor-T4approximant in the quasi-circular limit has the following expression: M d x d t (cid:12)(cid:12)(cid:12)(cid:12) = ηx (cid:32) (cid:88) k =2 a k x k (cid:33) . (6)This equation, when integrated, gives the evolution of the post–Newtonian variable x . Wegive in Appendix B the formulas for the expansion coefficients a i , and make available the Mathematica script where we implemented those coefficients and integrated eq. (6) numer-ically. Once x is known, we will obtain the orbital phase by integrating the equation: M d Φ orb dt = M ω orb = x / . (7)The binary orbit shrinks during the evolution, therefore the separation depends on timeas well. In order to describe this, we can calculate r directly from the post-Newtonianparameter x = v c = ω r c , using Kepler’s third law ω r = GM , to obtain r ( t ) = M x ( t ) − in geometrical units. We will push the precision in the calculation of the separation evenfurther, by applying post-Newtonian corrections up 3 P N : r = M ( r x − + r + r x + r x ) , (8)with the terms r i PN included in Appendix B.Once the evolution of the orbital phase and the separation are known, we can constructthe amplitude of the gravitational wave as a combination of two independent states ofpolarizations, similar to electromagnetic waves. For this, we use the general formula: h + = − M ηR (cid:26) (cos θ + 1) (cid:20) (cid:18) − ˙ r + r ˙Φ + Mr (cid:19) cos 2Φ + 2 r ˙ r ˙Φ sin 2Φ (cid:21) + (cid:18) − ˙ r + r ˙Φ + Mr (cid:19) sin θ (cid:27) h × = − M ηR cos θ (cid:20) (cid:18) − ˙ r + r ˙Φ + Mr (cid:19) sin 2Φ − r ˙ r ˙Φ cos 2Φ (cid:21) . (9)6he two polarization modes are denoted h + and h × , to emphasize that the angle betweenthem is π/
4, and not π/
2, like with electromagnetic waves. The amplitude of the gravita-tional waves depends on the orientation of the binary with respect to the detector. Whenworking in the detector’s coordinate system, also called the fundamental frame , we have totake into account the inclination angle θ , which is the angle between the binary orbital planeand the fundamental plane. We will assume an optimal orientation of the detector, normalto the orbital plane, so that the orbit’s inclination angle θ is zero. Then eq. (9) will become: h + = − M ηR (cid:20) (cid:18) − ˙ r + r ˙Φ + Mr (cid:19) cos 2Φ + 2 r ˙ r ˙Φ sin 2Φ (cid:21) , (10) h × = − M ηR (cid:20) (cid:18) − ˙ r + r ˙Φ + Mr (cid:19) sin 2Φ − r ˙ r ˙Φ cos 2Φ (cid:21) . (11)The waveform strain is constructed from the plus and cross polarization modes as follows: h ins ( t ) = h + ( t ) − ih × ( t ) . (12)We mention that the variable h (named strain or amplitude of the wave) is dimensionless,and represents the change in length divided by the length. For our convenience, we chose toexpress equation (12) in a more compact form, by transforming the sine and cosine termsinto their exponential forms. The final result is: h inspiral ( t ) = A ( t ) e − i φ ( t ) , and A = A + iA . (13)The amplitudes are given by A = − M ηR (cid:18) ˙ r + r ˙Φ + Mr (cid:19) , and A = − M ηR (cid:16) r ˙ r ˙Φ (cid:17) . (14) B. The Merger Model
The merger starts beyond the innermost stable circular orbit (ISCO), which is definedas the last complete orbit before the binaries plunge and collide. The radius of this orbitis proved to be r ISCO = 6 M in geometrical units (see Appendix A). The highly nonlinearmerger phase is correctly modeled only by General Relativity, and numerical simulationsof Einstein’s equations are necessary to provide accurate results. However, because of thecomplexity and cost of numerical simulations, semi-analytical models were developed forthe merger case, based on the results provided by numerical relativity. One of the most7uccessful techniques is the Implicit Rotating Source model, as is presented in , wherethe merger waveform is calculated by the analytical fit to numerical simulations. We buildup the waveform strain for the merger starting with the equation: h merger ( t ) = A ( t ) e − i Φ gIRS ( t ) . (15)We mention that we must use the formula given in for the amplitude of the strain, becausein the power factor of 1 / A ( t ) = A ω ( t ) (cid:12)(cid:12) ˙ˆ f (cid:12)(cid:12) α (cid:16) ˆ f − ˆ f (cid:17) / , (16)where: ˆ f = c (cid:18) κ (cid:19) κ (cid:34) − (cid:18) κ e − t/b (cid:19) − κ (cid:35) . (17)The angular orbital velocity is calculated with: ω ( t ) = ω QNM (cid:16) − ˆ f (cid:17) . (18)where ω QNM is the fundamental, or least damped frequency of the quasi-normal modes(QNM) emitted by the final black hole as it settles into its spherical shape. We use for itthe relation given in : ω QNM = 1 − .
63 (1 − ˆ s fin ) . (19)where ˆ s fin is the spin of final black hole:ˆ s fin = 2 √ η − η + 2379287 η − η . (20)The phase is obtained by integrating the orbital angular velocity:Φ gIRS ( t ) = (cid:90) tt ω ( t )d t, (21)The quantities: ˆ f , ˙ˆ f = d ˆ f / d t , and ˆ s fin are obtained through an analytic fit to the numericalrelativity results. The coefficients α , b , c and κ are smooth function of the symmetric mass-ratio η and are given in Appendix C. A is a parameter which we can choose to be unity. Thismodel applies to the merger of non-spinning compact binaries of different mass-ratios, andis called generic IRS (gIRS) model. We implement this simple model in
Mathematica anduse it to calculate the strain of the gravitational waves during the merger (see Appendix D).8
II. THE IMPLEMENTATION OF THE MODELSA. The Inspiral Gravitational Waveform
Before starting the implementation of the models presented above, we need to determinethe domain of the integration, ranging from an initial to a final value for the PN-parameter x . The lower boundary x is dictated by the threshold value of the frequency when thesignal enters the Advanced LIGO detection band. We consider this frequency as beingdetermined by the cut-off frequency due to the Earth’s seismic activity: f lowGW = 10 Hz . Itis worth mentioning that, because x is a unitless physical parameter, we will calculate itusing SI units. The orbital velocity corresponding to f lowGW is, from Kepler’s third law, equalto v = (cid:0) GM ω low (cid:1) / , where ω low = πf low , and M is given in units of the solar mass M (cid:12) (see Appendix A). Note that the frequency of the gravitational waves is twice the orbitalfrequency f GW = 2 f orb , a known feature of quadrupole radiation. This means that thegravitational wave signal goes through two maxima and two minima per one orbit of thebinary motion . With this expression for v , we can calculate the initial value for the PNparameter to be: x = (cid:16) v c (cid:17) = (cid:18) GM πf lowGW c (cid:19) / . (22)The upper boundary is determined by the radius of the last stable orbit of the binary: r ISCO = 3 R Sch = 6
GMc . The velocity corresponding to this orbit in the Newtonian approachis: v ISCO = (cid:112) GM/r
ISCO = c/ √
6. With this value for the velocity, we obtain the upperlimit for the PN parameter of 0 th order to be: x P NISCO = 1 /
6. We add a 2 nd order post-Newtonian correction, that introduces a dependence on the symmetric mass of the binary ,such that: x P NISCO = 16 (cid:18) η (cid:19) . (23)Next, using Kepler’s third law we calculate the frequency of the gravitational wave at ISCOfunction the binary mass, by expressing the mass of the back hole in units of time (seeAppendix A): f ISCO = (cid:115) GMπ r ISCO = c / GM π . (24)This shows that the frequency at the end of the inspiral scales inversely proportional tothe mass of the binary, therefore the smaller the mass, the higher the frequency of the9ravitational wave. After we determine the lower and upper bounds for the PN parameter x , we only need to choose a value for the symmetric mass ratio η , and then we can proceedwith the integration of eq. (6).Let’s start with a fiducial binary configuration of total mass M = 40, and equal masses m = m = 20 given in solar masses, corresponding to a symmetric mass ratio η = 0 . Mathematica . We set the lower boundary x t =0 = x . If we don’t give an upper bound-ary for the time, we run into a known issue in numerical analysis, where the step sizebecomes effectively zero , and the equation becomes stiff . Stiffness is a numerical propertyof differential equations, caused by a set of factors, such as the numerical method, initialconditions, or sudden changes in the solution due to singularities or sharp features in thesolution for the differential equation. All this drives the step size to increasingly small val-ues, until eventually it becomes effectively zero, which leads to unstable numerical results.We adopt the decision to go past the upper x P NISCO boundary for x , and evolve to a finaltime very close to the time t S when the equation becomes stiff. Our choice is motivated bythe fact that we included corrections up to 6PN order for x in eq. (6), and supplementedwith corrections up to 3PN for r in eq. (8) so we can explore solutions beyond the laststable orbit. To this extent we estimate the transition time from the stable binary blackhole system to the coalescence into a single black hole, (also known as the time of flight ),which is the time necessary for the two black holes to fall from the last stable orbit (ISCO)to the light ring. This region is roughly located at twice the Schwarzchild radius: r LR = 4 M for a slowly rotating black hole. While the event horizon is the invisible region around theblack hole from which no light can escape, the light ring is made visible by the light forcedto orbit around the black hole due to the lensing effect of the strong gravitational field. Thedistance from ISCO to the light ring is thus: r ISCO − r LR = 2 M , and the corresponding timeshould be in geometrical units t of = 2 M . We pick as the final time for the PN evolution thetime t F = t S − t of , and integrate numerically eq. (6) from t = 0 s to t F . Transformed inseconds by multiplication with M (cid:12) ( s ), the final time is t F = 11 . s , and corresponds to avalue x final = 0 . x P NISCO = 0 . M = 40 M (cid:12) , the inspiral gravitational wave signal ideally stays in the detectorrange for a total time of nearly 12 seconds, and has a frequency range from 10 Hz to about200 Hz . 10fter we obtain x ( t ), we proceed to calculate the evolution of the phase with time, whichis done by numerically integrating eq. (7). Next, we need to determine the evolution of thedistance between the orbiting black holes, known as the binary separation r ( t ). For this,we use eq (8), that gives r ( t ) corrected up the 3 rd PN order in x ( t ). We plot in Fig. 1the evolution in time of the PN parameter x ( t ) and of the separation r ( t ), for the lastsecond before the merger. Before calculating the amplitude of the gravitational wave, let’s t ( sec ) x t ( sec ) r separation FIG. 1: The evolution in time of the PN parameter x and separation r for the inspiral of ablack hole binary of total mass M = 40 M (cid:12) assume a realistic value for the distance R to the binary and take it to be the distance toAndromeda, the closest galaxy: R = 2 . × km , or 2 . r in seconds in eq. (10). We do this by multiplying it with the mass of the sun in seconds,as explained in Appendix A. We also need to express the total mass M in km if we want totake R in km . We calculate the plus and cross polarizations of the strain with eqs. (10), (11)and the amplitude with eq. (12). Our calculations show that the maximum amplitude of thegravitational wave strain is A max = | h inspiral ( t F ) | = 5 . × − , a shockingly small value.This makes sense if we remember that the strain of the detected signal GW150914 was assmall as 10 − , and the event was located at 1.3 billion light-years away! We plot in Fig. 2the amplitude and h + polarization mode of the gravitational wave for the last second beforethe merger, as it would be seen by an optimally oriented detector here on Earth.The maximum amplitude of the strain will be 1 and it’s reached in the immediate vicinityof the binary black hole. We will rescale the strain to unity by dividing it with A max , andplot it only the last 1 / h × mode lags behind the h + mode. Our11 t ( sec ) × - × - × - × - A inspiral t ( sec )- × - h + inspiral FIG. 2: The last second of the gravitational wave amplitude and h + component of thestrain from the inspiral of a black hole binary of total mass M = 40 M (cid:12) , at a distance R = 2 . × km = 2 . ,
000 parsecs (Andromeda Galaxy). t ( sec )- h inspiral FIG. 3: The evolution of the strain for the inspiral of an equal mass binary with total mass M = 40 M (cid:12) . The solid plot represents the h + , and the dashed plot is the h × polarization.next step is the calculation of the h spherical harmonic component of the strain for theinspiral model. Spherical harmonic functions form a complete set of orthogonal functionsdefined on the surface of a sphere: h = − M ηR e − i Φ (cid:114) π (cid:18) ( r ˙Φ + i ˙ r ) + Mr (cid:19) . (25)This is the dominant spherical harmonic mode in the gravitational wave signal. Indeed,we show in our Mathematica script that the difference between the strain calculated witheq. (25) and with the eq. (13) is in the roundoff error, which is a proof to our calculations.12 . The Merger Gravitational Waveform
The calculation of the gravitational wave strain for the merger proceeds in a straightfor-ward way. We pick the same equal mass configuration for the binary, and start by calculatingthe angular frequency with eq. (18), then we integrate it to obtain the phase Φ gIRS . Inspect-ing eq. (17) we see that the time is in geometric units, and we replace it with t → tMM (cid:12) ( s ) in order to revert to time measured in seconds. Next, we rescale the factor A → MM (cid:12) ( s ) torender the strain unitless. With the amplitude given by eq. (16), we compute the gravita-tional wave strain for the merger using eq. (15). We determine the maximum value for theamplitude of the merger strain, and rescale the amplitude to unity by dividing the strainwith the maximum amplitude. We plot in Fig. 4 the evolution of the amplitude with timefor a small time interval of ( − M, +100 M ) around the origin, which expressed in secondsis around ( − . , . s . We see that the peak of the amplitude is not at t = 0, but cor-responds to a retarded time t r = 0 . ms . We shift the merger strain with half that time: t r / . ms , bringing the highest maximum of the h + component to the time origin,and plot it in Fig 4. We see from Fig: 4 that the h + mode lags behind the h × polarization. - - t A merger - - t ( sec )- - h merger FIG. 4: The amplitude and the rescaled strain for the merger of an equal mass binary of M = 40 M (cid:12) . The solid plot represents the h + , and the dashed plot is the h × polarization.13 V. THE MATCHING TECHNIQUEA. Matching for a Template Waveform
Next we will present a technique for constructing a complete gravitational waveform,by fitting together the strains for the inspiral and the merger, described in Sec. III, forthe same binary configuration with total mass M = 40 M (cid:12) . The gIRS waveform is tunedto numerical relativity results to model the dynamics of the merger, but we observe fromFig. 4 that its amplitude diminished rapidly – an indication that its accuracy deterioratesafter only a few cycles, thus it has a very limited range of applicability. In developing ourmatching technique, we are relying on the accuracy of the inspiral evolution, which includedcorrections terms for the energy up to the 6 th order of the post-Newtonian approximation,and is evaluated up to the light ring.We start by determining the best matching interval by comparing the frequency evolutionat the end of the inspiral and at the beginning of the merging phases. We calculate f inspiralGW = ω inspiral /π and f mergerGW = ω merger / (2 π ). In fact, this relationship can be intuitively seenonly comparing eq. (13) with eq. (15) for the strain of the gravitational wave. In orderto synchronize the two models at a time consistent with their common evolution, we firstmatch them in frequency (see Fig. 5). To do this, we shift the time axis of the inspiralfrequency plot by − t F so that the end of the inspiral is at t = 0, and then we shift the timeaxis of the gIRS frequency plot by a time parameter τ , which is adjusted until the mergerfrequency at t = 0 is approximately equal to the last frequency for the inspiral. This makesthe frequency plot continuous between the inspiral and merger phases. We see that thefrequency of the merger phase increases abruptly, reaching more than twice the frequency atthe end of the inspiral. This sudden increase in frequency during the coalescence is knownas the chirp of the gravitational wave. After this the binary enters into the ringdown phase,which ends when the final black hole is formed. Fig. 5 shows the comparison between theshifted gravitational wave frequency at the end of the inspiral, the merger frequency, and thetranslated merger frequency with time τ , until the overlap with the frequency of the inspiralis reached. This is a straightforward technique and can be easily used in lab settings, orin hands-on demonstrations on gravitational waves. We obtain the frequency overlap for atime shift τ = 2 . ms , and apply this time shift to both polarization modes of the merger.14 - t ( sec ) FIG. 5: Comparison of the gravitational wave frequency at the end of inspiral andbeginning of merger for an equal mass binary of M = 40 M (cid:12) . The solid line is frequency atthe end of the inspiral, the long dashed line is the merger frequency, and the short dashedline (red) is the merger frequency shifted to overlap with the end inspiral frequency.Next we proceed to match the inspiral and merger strains. We analyze first the h + polarization and observe that we need to account for the phase difference between theinspiral and the merger, because the inspiral h + mode leads, while the merger h + lags. Wefind that a strain parameter Φ = π will bring the inspiral and merger h + strain in phase.We obtain a clear overlap at the last maximum of the inspiral strain, then adjust the timeaxis in increments of the retarded time t r /
2, which is about a tenth the time shift τ , untilthe peak of the inspiral overlaps with the peak of the merger waveform. Remarkably, withonly a time shift of ∆ t = 3 / t r we obtain a very good fit and we do not have to rescalethe amplitude of the strain at the overlapping point. The h × polarization is not affected bya phase difference, and an excellent overlap is obtained when we adjust the time axis withonly ∆ t = t r . We see from Fig. 6 that the matching interval is optimal, because we canpick other points in the vicinity of the peak and obtain the same high overlap between theinspiral and the merger amplitudes. 15 - - - - t ( sec )- - h + Strain - - - - - t ( sec )- - h x FIG. 6: The overlap of the strain for an equal mass binary of total mass M = 40 M (cid:12) . Thesolid black plot represents the inspiral strain, and the dashed red plot is the merger strain. B. Comparison with the GW150914 Waveform
We test our implementation by calculating the inspiral and merger strain for the binaryconfiguration GW150914, and comparing our results with the gravitational-wave strain tem-plate for this event, released by the Gravitational Wave Open Science Center . This pro-cedure will test both our implementation of the post-Newtonian and gIRS models and ouroverlapping technique, confirming its viability. The mass parameters for this binary configu-ration are m = 36 . M (cid:12) and m = 29 . M (cid:12) , with a symmetric mass ratio η = 0 . M = 62 . M (cid:12) ,16hich gives a Schwarzchild radius of only 1 . × km , while the distance to the event isestimated to be about 1 . × km . Using the technique described in Sec. III, with the valuefor the cutoff detector frequency f lowGW = 10 Hz , we obtain the strain of the gravitational wavefor the inspiral and merger phases of the binary evolution. The final time for the inspiralstrain t F = 5 . s and its maximum amplitude is A max = 1 . × − .The merger strain is calculated within a time interval of ( − . , . s and the retardedtime for the peak in amplitude is t r = 0 . ms . The frequency overlap is obtained for atime shift τ = 4 . ms , which is used to shift the merger strain before matching. Fig. 7 showsthe complete waveforms for the h + and h × polarizations, obtained with ∆ t = t r / h + mode with Φ = π , to obtain a very good overlap between the inspiral and merger.Lastly, we read in Mathematica the data for the open source GW150914 template, thenoverlap it with our calculated strain for the M = 65 . M (cid:12) binary. Fig. 8 shows an excel-lent match between our post-Newtonian waveform and the GW150914 template, for severalpeaks, with a time adjustment of only ∆ t = 3 / t r and an increase in amplitude of 20%. Ourmodel loses accuracy at the last peak, and this is expected, because we pushed our calcula-tion beyond the limit of its applicability. The gIRS strain for the merger is adjusted with∆ t = 3 / t r , which gives a good overlap with the GW150914 template at peak amplitude. V. CONCLUSIONS
Our objective of developing a streamlined matching method for the analytical calcula-tion of the complete gravitational waveform, while keeping it simple enough to be accessibleto undergraduate physics students, was accomplished. We implemented two analytical al-gorithms for calculating gravitational wave templates during the inspiral and merger ofcompact binary systems, and we built a cohesive method of combining them into a completewaveform. We bypass the complicated Einstein’s equations by using the post-Newtonian(PN) theory to model the inspiral phase and the Implicit Rotating Source (IRS) for themerger phase of the binary evolution. After building the inspiral and merger waveforms, wedevised our matching method and validated it by comparing our results with the waveformtemplate for GW150914, the first detection of gravitational waves. This is a rich and timely17 - - - - t ( sec )- - h + Strain - - - - - t ( sec )- - h x FIG. 7: The overlap of the strain for the binary of total mass M = 65 . M (cid:12) . The solidblack plot represents the inspiral strain, and the dashed red plot is the merger strain.topic, and our approach, accessible to undergraduate students, can easily be implementedin a special topics course or research project for a junior or senior physics students. Weprovide the Mathematica scripts and explain in Appendix D the start-up procedure to befollowed by beginners in this field in order to generate a complete waveform. There areseveral future projects that can be developed based on this report, among which are build-ing complete waveforms for all the detected signals, extending the inspiral mode to includenon-zero eccentricity, testing, improving and optimizing the matching technique, employingand testing a more realistic analytic model for the merger, adding a ringdown model, etc.18 - - - - t ( sec )- - h + Strain
FIG. 8: Comparison between the GW150914 template and our complete model for the h + polarization model. The solid black line is the GW150914 template, the long dash blue lineis the inspiral, and the small dash red line is the merger model. Our model was shifted intime by ∆ t = 1 . ms and its amplitude was increased by 20%. Appendix A: Geometrical Units
Throughout this paper we are working in
Geometrical Units (GU), in which calculationsare simplified, because we don’t have to deal with physical constants such as the universalgravitational constant ( G ) or the speed of light ( c ), because they are set to unity ( G = c = 1).Let’s first set the speed of light to unity. Then we can measure time in units of distance: c = 1 = 2 . × m/s → s = 2 . × m , m = 3 . × − s . (A1)Now by setting Newton’s gravitational constant to unity, and taking the unit for time asmeasured in meters, we can measure mass in units of distance as well: G = 1 = 6 . × − m / kg · s → kg = 0 . × − m = 0 . × − km . (A2)Therefore, we measure both mass and time in meters , which is a distance, or a geometrical unit. In order to establish a straightforward correspondence between theoretical geometricalunits and observations, we express all the relevant quantities in units of solar masses and19se this quantity when converting back and forth between those systems of units. To thispurpose, let’s calculate the mass of the Sun, M (cid:12) , in km and in s , by using the tricks weexplained above in eq. (A1) and eq. (A2) to transform kg into km and into s . Thus M (cid:12) canbe used as a universal unit, that measures mass, distance, and even time. M (cid:12) = 1 . × kg = 1 . km = 4 . × − s . (A3)Let’s explain how can we use M (cid:12) as the unit for mass, distance and time, with this con-version, by giving a few examples. The mass of a black hole measured in solar masses issimply mM (cid:12) , where m is a dimensionless multiplication number, and M (cid:12) is the mass of theSun in kg . As another example, let’s recall that the radius of a black hole for which theescape speed is equal to the speed of light, called the Schwarzchild radius, is R Sch = 2 Gc M .The Schwarzchild radius becomes simply R Sch = 2 M in geometrical units. For the massof the Sun, this radius is R Sch, (cid:12) = 2 . km . Now, by only expressing the black hole massin terms of the solar mass M (cid:12) written as unit of distance, we get back to the SI units: R Sch = 2
M M (cid:12) ( km ) = 2 . M km . Appendix B: Post-Newtonian Coefficients
We give below the coefficients used in eq. (6) for the calculation of the post-Newtonianvariable x . a = 170 . − . η + 370 . η − . η − . η + (14 . − . η ) log( x ( t )) a / = 1047 . − . η + 923 . η + 22 . η − .
446 log( x ( t )) a = 714 . − . η + 3058 . η − . η + 29 . η − . η + ( − . . η + 1146 . η ) log( x ( t )) a / = 3622 . − . η + 12973 . η − .η + 25 . η + (83 . − . η ) log( x ( t )) a = 11583 . − . η + 33371 . η − . η + 648 . η − . η − . η + ( − .
61 + 7001 . η − . η − . η ) log( x ( t )) + 33 . x ( t ) )20elow are coefficients used in eq. (8) for the calculation of the separation r ( t ) as an expansionin the post-Newtonian variable x . r P N = 1 r P N = − . ηr P N = 4 . η + 0 . η r P N = − . η − . η + 0 . η Appendix C: gIRS Coefficients
Here are the coefficients used in the calculation of the merger waveform. Q (ˆ s fin ) = 2(1 − ˆ s fin ) . α ( η ) = 1 Q (ˆ s fin ) (cid:18) η (cid:19) b ( η ) = 16014979 − η c ( η ) = 206903 + 1801141 √ η + 4241205 η log ( η ) κ ( η ) = 7131056 − η Appendix D: Procedures
We lay out below the step-by-step procedure to be followed by beginners in this field,including students with little or no knowledge of gravitational waves or
Mathematica , inorder to generate complete waveforms. One word of caution: this model is tailored to workbest for binary configurations of comparable mass ratio. Experiment with it and please donot hesitate to send an email to the first author if you run into a problem.1. Setting up: Install the
Mathematica software, then go to https://github.com/mbabiuc/MathScripts , click on the Clone or Download button on the right, andchoose the Download ZIP option.2. Mass parameters: Open the Match40.nb script on your computer and save it under adifferent name. Then under the section
Setting up the Mass Parameters change the21ass parameters m1 and m2 with values of your choice, and run the script (click onEvaluation and chose Evaluate Notebook ).3. Final integration time: Scroll down to the section titled Setting up the final integra-tion time text line, and note the error given by the function
NDSolve: At t == ...,step size is effectively zero; singularity or stiff system suspected . Copyand paste the numerical value for that time to the tS variable defined immediatelybelow the NDSolve and run the script again.4. Inspiral Waveform: Scroll to the section
Calculation of the Inspiral Waveform and ifyou want, change the variable R to a realistic value for the distance from the detectorthe to source. Now run only that portion of the script again (press Shift + Enter )until the section Calculations for Merger-Ringdown Waveform . Congratulations, youjust generated your first gravitational wave model for the inspiral! If you want tosave the plot, you will type in the notebook, just below the plot, the command
Export["/Path of File/Name of File", %, "PDF"] and run it.5. Merger Waveform: This is generated without any intervention.6. Matching in frequency: In the
Manipulate plot, click the + button, then play thegraph until the plot goes through the origin of the axes. The shift is done in incrementsof the retarded time t r . Divide the value obtained for τ by t r , to obtain the factor χ ,and check if the frequency of the merger at that time f M is nearly equal to the lastfrequency of the inspiral f F . If not, tweak the factor χ to obtain the best concordance.7. Matching in amplitude: Lastly, the matching in amplitude should follow straight fromthe matching in frequency. It might require a slight adjustment of the time axis for themerger model in increments of t r . This is done by tweaking the (cid:15) + and (cid:15) × coefficients.We hope that this procedure will be easy to follow and rewarding, and will be useful inbootstrapping future projects in gravitational waves with undergraduates, increasing theinvolvement of the physics students and faculty in this new and exciting field.22 CKNOWLEDGMENTS
This work was supported by the department of Physics at Marshall University and theNSF EPSCoR Grant OIA-1458952 to the state of West Virginia “Waves of the Future”. B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), “Observation ofGravitational Waves from a Binary Black Hole Merger,” Phys. Rev. Lett. B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), “GW151226:Observation of Gravitational Waves from a 22-Solar-Mass Binary Black Hole Coalescence,”Phys. Rev. Lett. B. P. Abbott et al. (LIGO Scientific and Virgo Collaboration), “GW170104: Observation of a50-Solar-Mass Binary Black Hole Coalescence at Redshift 0.2,” Phys. Rev. Lett. B. P. Abbott et al. (LIGO Scientific and Virgo Collaboration), “GW170608: Observation of a19 Solar-mass Binary Black Hole Coalescence,” ApJ Letters
L35 (2017). B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), “GW170814: AThree-Detector Observation of Gravitational Waves from a Binary Black Hole Coalescence,”Phys. Rev. Lett. B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), “GW170817:Observation of Gravitational Waves from a Binary Neutron Star Inspiral,” Phys. Rev. Lett. W. E. East, S. T. McWilliams, J. Levin, F. Pretorius, “Observing complete gravitational wavesignals from dynamical capture binaries,” Physical Review D. (4) (2013). E. A. Huerta, P. Kumar, B. Agarwal, D. George, H-Y Schive, H. P. Pfeiffer, R. Haas, W. Ren, T.Chu, M. Boyle, D. A. Hemberger, L. E. Kidder, M. A. Scheel, B. Szil´agyi, “ Complete waveformmodel for compact binaries on eccentric orbits,” Physical Review D (2) (2017). J. Blackman, S. E. Field, M. A. Scheel, C. R. Galley, C. D. Ott, M. Boyle, L. E. Kidder, H. P.Pfeiffer, B. Szil´agyi, “Numerical relativity waveform surrogate model for generically precessingbinary black hole mergers,” Physical Review D (2) (2017). E. Poisson, W. M. Clifford,
Gravity: Newtonian, Post-Newtonian, Relativistic (Cambridge Uni- ersity Press, 2014). R. M. Wald,
General relativity (University of Chicago Press, 1984). F. Pretorius, “ Evolution of Binary Black-Hole Spacetimes,” Physical Review Letters (12)(2005). L. Blanchet, “ Gravitational Radiation from Post-Newtonian Sources and Inspiralling CompactBinaries,” Living Reviews in Relativity (1) 187 (2014). A. Einstein, L. Infeld, B. Hoffman, “The gravitational equations and the problem of motion,”Annals of Mathematics (1) 65 (1938). W. M. Clifford, “On the unreasonable effectiveness of the post-Newtonian approximation ingravitational physics,” Proc. Nat. Acad. Sci. (US) A. Buonanno, G. B. Cook, F. Pretorius, “Inspiral, merger, and ring-down of equal-mass black-hole binaries,” Physical Review D (12) (2007). E. A. Huerta, P. Kumar, S. T. McWilliams, R. O’Shaughnessy, N. Yunes, “Accurate and efficientwaveforms for compact binaries on eccentric orbits,” Physical Review D (8) (2014). Y. Itoh, T. Futamase, “New derivation of a third post-Newtonian equation of motion for rela-tivistic compact binaries without ambiguity,” Phys.Rev. D P. Ajith, N. Fotopoulos, S. Privitera, A. Neunzert, N. Mazumder, A. J. Weinstein, “ An effectualtemplate bank for the detection of gravitational waves from inspiralling compact binaries withgeneric spins,” Physical Review D (8) (2014). J. G. Baker, W. D. Boggs, J. Centrella, B. J. Kelly, S. T. McWilliams, J. R. van Meter, “Mergersof nonspinning black-hole binaries: Gravitational radiation characteristics,” Phys.Rev.D I. Hinder, F. Herrmann, P. Laguna, D. Shoemaker, “Comparisons of eccentric binary black holesimulations with post-Newtonian models,” Physical Review D B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), “The basic physicsof the binary black hole merger GW150914,” Annalen der Physik (1-2) 1600209 (2017). Gravitational Wave Open Science Center, . Wolfram Mathematica: Modern Technical Computing, ..