Detailed explanation of phase_frequency() routine

0 The phase_frequency() routine starts with inputs describing the physical properties of the system (the masses) and an initial frequency from which to start the evolution. We then compute the orbital frequency evolution [in cycles/second] directly from the formula given in [7]

$$f(t) = {M_\odot \over 16 \pi T_\odot m_{\rm tot}} \biggl\{ \Theta^{-3/8}... ... 2688} + {11\over 32}\eta \right) \Theta^{-5/8} - {3\pi \over 10} \Theta^{-3/4}$$

$$\qquad\qquad\qquad+ \left( {1855099\over 14450688} + {56975\over 258048} \eta + {371\over 2048} \eta^2 \right) \Theta^{-7/8} \biggr\} \;,$$ (6.4.33)

where $m_{\rm tot}$ is the total mass of the binary. The time integral of this equation gives the orbital evolution in cycles. Multiplying by $2\pi$ yields the orbital phase in radians

$$\phi (t) = \phi_c -{1\over\eta} \biggl\{ \Theta^{5/8} +\left({3715\over 8064} + {55\over 96}\eta \right) \Theta^{3/8} - {3\pi \over 4} \Theta^{1/4}$$

$$\qquad\qquad\quad+\left({9275495\over 14450688}+{284875\over 258048} \eta\ +{1855\over 2048} \eta^2 \right) \Theta^{1/8} \biggr\}\;.$$ (6.4.34)

Here $\Theta$ is a dimensionless time variable

$$\Theta={\eta M_\odot \over 5 T_\odot m_{\rm tot}} (t_c-t) \; ,$$ (6.4.35)

$\eta = \mu/m_{\rm tot}$, and $t_c$ is the time of coalescence of the two point masses. Similarly the constant $\phi_c$ is the phase at coalescence, which is arbitrarily set in phase_frequency() so that $\phi=0$ at the initial time. [See the detailed discussion of the phase conventions below.] Also notice that the mass quantities only appear as ratios with the solar Mass $M_\odot$, and the time only appears as a ratio with the quantity $T_\odot = 4.925491\times 10^{-6}\> \rm sec$ in Eq.().

These formulations of the post-Newtonian equations for the phase and frequency are simple to implement: each pass through the loop increments the time by the sample time (Sample_Time in the example) and computes the phase and frequency using Eqs. () and (). However, there is an alternative formulation. In deriving these equations the ``natural'' equation that arises is of the form $\dot f= F(f)$. [See e.g. [12] Eq.(3).] This in turn can be integrated to give an equation of the form $t_c-t = T(f)$. In our formulation this equation has been inverted - throwing away higher-order post-Newtonian terms as you go - to give Eq.(). However the equation in the form $t_c-t = T(f)$ can also be implemented directly. In this type of formulation one would again increment the time, but then use a root-finding routine to find the frequency at each time step. Our chosen method has the advantage of avoiding a time-consuming root-finder at each time step; however the alternative formulation has undergone fewer damaging post-Newtonian transformations, and may therefore be more accurate.

In our formulation we only need to call a root-finding routine at the start of the chirp to find the value of $t_c-t$ when the system is at the initial frequency. In order to insure that we find the correct root for the starting time we begin a search at a time when the leading order prediction of the frequency is well below the desired starting frequency. We step forward in time until we bracket the root; we then call the Numerical Recipes root-finder rtbis() to compute the root precisely. This is depicted in the lower right corner of figure where we show the value of the ``time'' coordinate $X$ that corresponds to an initial frequency of $60$Hz. This method is virtually assured of finding the correct root in that it will find the first solution as we proceed from right to left in figure . The primary problem in finding this root is that there may actually be no meaningful start-time for the specified chirp. For example, if you you were to specify a chirp with two $1.4M_\odot$ objects with an initial frequency of 1000Hz, you can see from the figure that there is no value of $X$ (i.e. $t_c-t$) that corresponds to this frequency. In this case phase_frequency() will search from right to left for the start time. It will notice that it is passing over the peak in the graph and out of the regime of post-Newtonian viability. It will then terminate the search and notify the caller that there is no solution for the requested chirp.

The behavior of the frequency equation is shown in figure . As time increases the frequency rises to a maximum and then begins to decrease dramatically. Notice that the maximum occurs when the dimensionless time parameter $\Theta={\eta(t_c-t) \over 5 T_\odot m_{\rm tot} } = X^8$ is approximately unity; this feature is only weakly dependent on the mass ratio. The fact that $\Theta \approx 1$ means the post-Newtonian corrections in Eq.() are comparable to the leading order term. Therefore, this peak is a natural place to terminate the post-Newtonian chirp approximation. In the example the code terminated the chirp for precisely this reason. [See the warning message.]

Although it is not shown in the figure the behavior of $f$ as $X$ nears zero is very abrupt; the function goes sharply negative and then turns around and diverges to $+\infty$ as $X\rightarrow 0$ (i.e. $t\rightarrow t_c$). This abrupt behavior will happen on a time scale of order $T_\odot$ (a few microseconds). Typical sample times are likely to be on the order of a tenth of a millisecond, and therefore the iterative loop may step right over this maximum-minimum-divergence behavior of the frequency function altogether. Don't worry. The routine phase_frequency() handles this case gracefully. The routine will stop the chirp calculation and warn the caller if the time stepper goes beyond the coalescence time. It will also stop the chirp calculation if it senses that the time has stepped over the dip in frequency and is on the strongly divergent part of the frequency curve near the $X=0$ axis.

Figure: Orbital frequency as a function of the ``time'' coordinate $X=\biggl ( {\eta (t_c-t) M_\odot \over 5 T_\odot m_{\rm tot} } \biggr )^{1/8}$.