The waveform

The Teukolsky formalism gives the two components $h_+$ and $h_\times$ of the gravitational waves produced by a test mass in a fixed, circular orbit of radius $r_0$ as

$$h_+(u) - i h_\times(u) = \frac{2 G \mu}{r c^2} \sum_{l\geq ... ...^{-i m \Omega (t - c r^*)} {}_{-2}Y_{lm}(\vartheta, \varphi),$$ (7.1.159)

where $u=t-c r^*$ is the standard retarded time. Here $G$ is Newton's constant and $c$ the speed of light. Here the functions ${}_{-2}Y_{lm}(\vartheta,\varphi)$ are the spherical harmonics of spin weight $-2$, $\Omega$ is the angular velocity and $M=m_1+m_2$ is the total mass of the system. The angles $\vartheta$ and $\varphi$ are the usual spherical coordinates, as defined in Section , Figure ( $\vartheta = \iota$ and $\vartheta = \beta$). Since the motion is symmetric around the $z$-axis $\varphi$ does not have an intrinsic meaning. Throughout this section angles are measured in radians. The amplitudes $A_{lm}$ are constants and have to be calculated numerically by solving the Teukolsky equation. We provide $A_{lm}$ tabulated in a datafile [45] as a function of the orbital velocity $v$ (measured in units of the speed of light $c$)

$$c v = r_0 \Omega = \sqrt{\frac{G M}{r_0}} = (G M \Omega )^\frac{1}{3} = (2 \pi G M f)^{\frac{1}{3}}.$$ (7.1.160)

Here $f$ is the orbital frequency measured in units on Hertz.

To account for the decay of the orbit due to the emission of gravitational waves we use an adiabatic approximation: The energy radiated away as calculated from expression () is used to calculate the change in orbital frequency. By doing so, $\Omega$, and thus also $v$, become functions of time and we have to replace the product $\Omega u$ by an integral $\int \Omega(u) du$.

Following [46] we find that the time evolution of the velocity is governed by

$$\dot{v} = \frac{32}{5} \frac{\mu c^3}{G M^2} v^9 \frac{P(v)}{Q(v)},$$ (7.1.161)

where $\mu = \nicefrac{m_1 m_2}{M}$ is the reduced mass of the system. The function $P(v)$, which is determined by the $A_{lm}$'s, is defined by the equation

$$\dot{E} = \frac{dE}{dt} = P(v) \left(\frac{dE}{dt}\right)_N$$

where $(dE/dt)_N = \nicefrac{-32 \mu^2 v^{10} c^4}{5 G M^2}$ is the quadrupole-formula expression for the gravitational-wave luminosity. We also provide $P(v)$ in tabulated form [45]. The function $Q(v) = (1-6v^2)(1-3v^2)^{-\frac{3}{2}}$ relates the energy and frequency of the orbiting particle:

$$\frac{dE}{df} = E(v) \left(\frac{dE}{df}\right)_N.$$

This can be calculated by solving the geodesic equation for circular orbits in the Schwarzschild spacetime.

Note that $v(t)$ can be calculated from equation () once and for all: First (numerically) calculate the solution $V(t)$ of equation () with the factor $\nicefrac{\mu c^3}{G M^2}$ set to one; then the solution $v(t)$ for general $\mu/M^2$ is simply

$$v(t) = V\left(\frac{\mu c^3}{G M^2} t\right).$$ (7.1.162)

As mentioned before, we have to replace the product $\Omega u$ in equation () by an appropriate integral. We first note that we want to look at waves at a fixed radius $r^* \rightarrow \infty$. We can thus simply ignore the dependence on $r^*$ because it will only contribute a fixed phase $\varphi_0$. Since our data are tabulated as functions of the orbital velocity $v$, it is convenient to use $v$ instead of the time $t$ as an independent variable. This is permissible because $v$ depends monotonically on $t$. Thus the phase becomes:

$$-i m \Omega t \ \ \longrightarrow \ \ -i m \frac{M}{\mu} \f... ...v' \frac{Q(v')}{P(v') v'^6} =: -i m \frac{M}{\mu} \Phi(v_0,v).$$ (7.1.163)

It is important to note that $\Phi$ is not unique, but depends on an arbitrary parameter $v_0$. As with $r^*$, a change in $v_0$ will only affect the phase $\varphi_0$. Since $\Phi$ can become rather large, the freedom in choosing $v_0$ can be used to keep $\Phi$ small in the region of interest. A good choice of $v_0$ will improve the numerical accuracy tremendously. (For example, the standard trigonometric functions in C become hopelessly inaccurate for arguments $> 10^6$ for floats and $> 10^{10}$ for doubles).

The signal is now given by

$$h_+(t) - i h_\times(t) = \frac{2 G \mu}{r c^2} \sum_{l\geq ... ...c{M}{\mu} \Phi(v_0, v(t))} {}_{-2}Y_{lm}(\vartheta, \varphi),$$ (7.1.164)

where $v(t)$ is given by equation ().

To extract $h_+$ and $h_\times$ independently we use the fact that $A_{l-m} = (-1)^l \overline{A}_{lm}$ and, since ${}_{-2}Y_{lm}(\vartheta, \varphi) = {}_{-2}Y_{lm}(\vartheta,0)\, e^{i m \varphi}$ and ${}_{-2}Y_{lm}(\vartheta,0)$ is real, we have ${}_{-2}Y_{l-m}(\vartheta, \varphi)= {}_{-2}Y_{lm}(\vartheta,0)\, e^{-i m \varphi}$ [26]. We can now split the sum () into real and imaginary part. This gives the two components as

$$h_+ = \frac{2 G \mu}{r c^2} \sum_{2 \leq l, 1 \leq m \leq l } \le... ...t) \left({}_{-2}Y_{lm}(\vartheta,0) + (-1)^l {}_{-2}Y_{l-m}(\vartheta,0)\right)$$

$${}$$ (7.1.165)

$$h_{\times} = \frac{2 G \mu}{r c^2} \sum_{2 \leq l, 1 \leq m \leq ... ... \left(-{}_{-2}Y_{lm}(\vartheta,0) + (-1)^l {}_{-2}Y_{l-m}(\vartheta,0)\right).$$

Here we have introduced the notation $\Re_{lm} := {\rm Re}\, A_{lm}$, $\Im_{lm} := {\rm Im}\, A_{lm}$ and $s_m := \sin(-m (\Phi M/\mu - \varphi))$ and $c_m := \cos(-m (\Phi M/\mu - \varphi))$.

The GRASP routine which calculates $h_+$ and $h_\times$ uses expression () truncated to a finite number of terms determined by the user. Finally we note that a change in $\varphi$ has the same effect on the waveforms as a change in $r^*$ and $v_0$.