 |
LALInspiral 5.0.3.1-ea7c608
|
|
|
LALInspiral 5.0.3.1-ea7c608
|
|
Provides routines to generate restricted parametrized post \({}^{5/2}\)-Newtonian inspiral waveforms.
This header covers routines to generate a "restricted" parametrized post \({}^{5/2}\)-Newtonian binary inspiral waveform in the time domain. That is, the calculation of the wave phase is accurate to post \({}^{5/2}\)-Newtonian order (including corrections up to order \(v^5/c^5\), where \(v\) is the orbital speed), but the wave amplitudes are accurate only to leading (post \({}^0\)-Newtonian) order. Furthermore, at each order the post \({}^{n/2}\)-Newtonian correction can be turned on, off, or set to an unphysical value, by adjusting a parameter \(p_n\).
The post-Newtonian expansion implicitly assumes an adiabatic inspiral, where one can represent the waveform by an "instantaneous" amplitude and frequency that vary over timescales longer than one wave period. The orbital frequency of the system to post \({}^{5/2}\)-Newtonian order is given in Eqs. 6.4.1 and 6.9.1 of [1]; here we work entirely in terms of the gravitational-wave frequency, which is twice the orbital frequency:
\begin{eqnarray} \label{eq_ppn_freq} f(t) & = & \frac{M_\odot}{8\pi T_\odot m_\mathrm{tot}}\left\{ p_0\Theta^{-3/8}+ p_1\Theta^{-1/2}+ p_2\left(\frac{743}{2688}+\frac{11}{32}\eta\right)\Theta^{-5/8}- p_3\frac{3\pi}{10}\Theta^{-3/4} \right. \\ & & \left.+ p_4\left(\frac{1855099}{14450688}+\frac{56975}{258048}\eta+ \frac{371}{2048}\eta^2\right)\Theta^{-7/8}- p_5\left(\frac{7729}{21504}+\frac{3}{256}\eta\right)\pi\Theta^{-1} \right\} \; , \end{eqnarray}
where \(M_\odot\) is the mass of the Sun, \(T_\odot=GM_\odot/c^3=4.925491\times10^{-6}\)s is the "geometrized" solar mass in time units, \(m_\mathrm{tot}=m_1+m_2\) is the total mass of the binary, \(\eta=m_1m_2/m_\mathrm{tot}^2\) is the (symmetric) mass ratio parameter, and \(\Theta\) is a dimensionless time parameter:
\begin{equation} \label{eq_ppn_theta} \Theta(t) = \frac{\eta M_\odot}{5T_\odot m_\mathrm{tot}}(t_c-t) \; . \end{equation}
Here \(t_c\) is the time of coalescence of the two masses in the point-mass approximation. The post-Newtonian parameters \(p_k\) are defined such that in a normal (physical) post \({}^{n/2}\)-Newtonian expansion, one sets \(p_1=0\) and \(p_{k>n}=0\), and \(p_k=1\) for all other \(k\). However, changing this convention can be used to model in an approximate way things such as spin, eccentricity, or non-GR theories of gravity. We also note that while most terms are normalized to their normal post-Newtonian values, the normalization on the \(p_1\) term is completely arbitrary, since it is zero in a normal post-Newtonian expansion.
The wave phase as a function of time can be computed analytically from Eq. \eqref{eq_ppn_freq} as \(\phi_\mathrm{orb}=2\pi\int f\,dt\):
\begin{eqnarray} \label{eq_ppn_phi} \phi(t) & = & \phi_c - \frac{2}{\eta}\left\{ p_0\Theta^{5/8}+ p_1\frac{5}{4}\Theta^{1/2}+ p_2\left(\frac{3715}{8064}+\frac{55}{96}\eta\right)\Theta^{3/8}- p_3\frac{3\pi}{4}\Theta^{1/4} \right. \\ & & \left.+ p_4\left(\frac{9275495}{14450688}+\frac{284875}{258048}\eta+ \frac{1855}{2048}\eta^2\right)\Theta^{1/8}- p_5\left(\frac{38645}{172032}+\frac{15}{2048}\eta\right)\pi \log\left(\frac{\Theta}{\Theta_0}\right)\right\} \; . \end{eqnarray}
Here \(\Theta_0\) is an arbitrary constant; changing it is equivalent to changing \(\phi_c\). We note that the post \({}^{5/2}\)-Newtonian term introduces a late-time divergence in phase which renders meaningless the interpretation of \(\phi_c\) as "phase at coalescence"; in our convention we define \(\phi_c\) to correspond to the case \(\Theta_0=1\).
We refer the interested reader to Sec. 6.6 of\ [1] for a discussion of how propagation effects shift the phase of the waveform relative to the orbital phase. To summarize, though: A changing propagation delay does introduce a time-dependent phase shift in the waveform, but the dependence on \(t\) is weak except at very late times; although it looks like a post \({}^{3/2}\)-Newtonian phase correction, it can in fact be represented as a post \({}^{3}\)-Newtonian phase correction combined with a post \({}^{3/2}\)-Newtonian amplitude correction. Since we are concerned with restricted post \({}^{5/2}\)-Newtonian waveforms, which model the amplitude only to leading (post \({}^0\)-Newtonian) order, we can ignore these propagation effects.
To leading order, then, the amplitude of the + and \(\times\) polarizations of the wave are given by Eqs. 6.6.1–6.6.4 of\ [1] as:
\begin{eqnarray} \label{eq_ppn_aplus} A_+(t) & = & -\frac{2T_\odot c}{D}(1+\cos^2 i) \left(\frac{\eta m_\mathrm{tot}}{M_\odot}\right) \left[\frac{\pi T_\odot m_\mathrm{tot}f(t)}{M_\odot} \right]^{2/3} \; , \\ \label{eq_ppn_across} A_\times(t) & = & -\frac{2T_\odot c}{D}(2\cos i) \left(\frac{\eta m_\mathrm{tot}}{M_\odot}\right) \left[\frac{\pi T_\odot m_\mathrm{tot}f(t)}{M_\odot} \right]^{2/3} \; , \end{eqnarray}
where \(D\) is the distance to the source and \(i\) is the inclination of the axis of the source to the line of sight. The normal polarization convention in\ [24] is used, where the reference \(x\)-coordinate axis for the + and \(\times\) polarization tensors is the ascending node of the rotational plane as it crosses the plane transverse to the propagation direction. This convention implies that the + and \(\times\) waveforms are elliptically polarized as follows:
\begin{eqnarray} \label{eq_ppn_hplus} h_+(t) & = & A_+(t)\cos\phi(t) \; , \\ \label{eq_ppn_hcross} h_\times(t) & = & A_\times(t)\sin\phi(t) \; . \end{eqnarray}
Prototypes | |
| void | LALGeneratePPNInspiral (LALStatus *, CoherentGW *output, PPNParamStruc *params) |
| Computes a parametrized post-Newtonian inspiral waveform. More... | |
| void | LALGeneratePPNAmpCorInspiral (LALStatus *, CoherentGW *output, PPNParamStruc *params) |
| Computes a parametrized post-Newtonian inspiral waveform with ampltidude corrections. More... | |
| void | LALGetInspiralParams (LALStatus *, PPNParamStruc *output, GalacticInspiralParamStruc *input, RandomParams *params) |
| Computes the input parameters for a PPN inspiral. More... | |
Data Structures | |
| struct | PPNParamStruc |
| This structure stores the parameters for constructing a restricted post-Newtonian waveform. More... | |
| struct | GalacticInspiralParamStruc |
| This structure stores the position and mass parameters of a galactic inspiral event. More... | |
Files | |
| file | GeneratePPNAmpCorInspiralTest.c |
| Generates a parametrized post-Newtonian inspiral waveform. | |
| file | GeneratePPNInspiralTest.c |
| Generates a parametrized post-Newtonian inspiral waveform. | |
Error Codes | |
| #define | GENERATEPPNINSPIRALH_ENUL 1 |
| Unexpected null pointer in arguments. More... | |
| #define | GENERATEPPNINSPIRALH_EOUT 2 |
| output field a, f, phi, or shift already exists More... | |
| #define | GENERATEPPNINSPIRALH_ETBAD 3 |
| Bad sampling interval. More... | |
| #define | GENERATEPPNINSPIRALH_EFBAD 4 |
| Bad starting frequency; could not get valid start time. More... | |
| #define | GENERATEPPNINSPIRALH_EPBAD 5 |
| Bad post-Newtonian parameters. More... | |
| #define | GENERATEPPNINSPIRALH_EMBAD 6 |
| Bad masses. More... | |
| #define | GENERATEPPNINSPIRALH_EDBAD 7 |
| Bad distance. More... | |
| #define | GENERATEPPNINSPIRALH_EMEM 8 |
| Out of memory. More... | |
More Termination conditions | |
In addition to the error conditions above, there are a number of ways that the signal generation routine can terminate gracefully while still returning a valid waveform. In many cases one wants to continue generating a waveform "until things fall apart"; the following codes, returned in the For the sake of LAL namespace conventions, these termination codes are | |
| #define | GENERATEPPNINSPIRALH_EFSTOP 0 |
| Reached requested termination frequency. More... | |
| #define | GENERATEPPNINSPIRALH_ELENGTH 1 |
| Reached maximum length, or end of provided time series vector. More... | |
| #define | GENERATEPPNINSPIRALH_EFNOTMON 2 |
| Frequency no longer increasing monotonically. More... | |
| #define | GENERATEPPNINSPIRALH_EPNFAIL 3 |
| Evolution dominated by higher-order PN terms. More... | |
| #define | GENERATEPPNINSPIRALH_ERTOOSMALL 4 |
| Orbital radius too small for PN approximation. More... | |
| void LALGeneratePPNInspiral | ( | LALStatus * | stat, |
| CoherentGW * | output, | ||
| PPNParamStruc * | params | ||
| ) |
Computes a parametrized post-Newtonian inspiral waveform.
This function computes an inspiral waveform using the parameters in *params, storing the result in *output.
In the *params structure, the routine uses all the "input" fields specified in GeneratePPNInspiral.h, and sets all of the "output" fields. If params->ppn=NULL, a normal post \({}^2\)-Newtonian waveform is generated; i.e. \(p_0=1\), \(p_1=0\), \(p_2=1\), \(p_3=1\), \(p_4=1\), \(p_{5+}=0\).
In the *output structure, the field output->h is ignored, but all other pointer fields must be set to NULL. The function will create and allocate space for output->a, output->f, and output->phi as necessary. The output->shift field will remain set to NULL, as it is not required to describe a nonprecessing binary.
This function is a fairly straightforward calculation of Eq. \eqref{eq_ppn_freq}– Eq. \eqref{eq_ppn_across} in GeneratePPNInspiral.h. However, there are some nontrivial issues involved, which are discussed in some depth in Secs. 6.4, 6.6, and 6.9.2 of [1] . What follows is a brief discussion of these issues and how this routine deals with them.
Computing the start time:
When building a waveform for data analysis, one would generally like to start the waveform at some well-defined frequency where it first enters the band of interest; one then defines the start time of the integration by inverting Eq. \eqref{eq_ppn_freq} if GeneratePPNInspiral.h. The current algorithm follows this standard approach by requiring the calling routine to specify params->fStartIn, which is then inverted to find \(\Theta_\mathrm{start}\). This inversion is in fact the most algorithmically complicated part of the routine, so we will discuss it in depth.
To help clarify the problem, let us rewrite the equation in dimensionless parameters \(y=8\pi fT_\odot m_\mathrm{tot}/M_\odot\) and \(x=\Theta^{-1/8}\):
\begin{equation} \label{eq_ppn_fnorm} y = \sum_{k=0}^{5} C_k x^{k+3} \; , \end{equation}
where:
\begin{eqnarray} C_0 & = & p_0 \;,\\ C_1 & = & p_1 \;,\\ C_2 & = & p_2\left(\frac{743}{2688}+\frac{11}{32}\eta\right) \;,\\ C_3 & = & p_3\frac{3\pi}{10} \;,\\ C_4 & = & p_4\left(\frac{1855099}{14450688}+\frac{56975}{258048}\eta+ \frac{371}{2048}\eta^2\right) \;,\\ C_5 & = & p_5\left(\frac{7729}{21504}+\frac{3}{256}\eta\right)\pi \;. \end{eqnarray}
We note that \(x\) is a time parameter mapping the range \(t=-\infty\rightarrow t_c\) to \(x=0\rightarrow\infty\).
In a normal post-Newtonian expansion it is possible to characterize the general behaviour of this equation quite accurately, since the values of \(p_k\) are known and since \(\eta\) varies only over the range \([0,1/4]\). In a parametrized post-Newtonian expansion, however, even the relative orders of magnitude of the coefficients can vary significantly, making a robust generic root finder impractical. However, we are saved by the fact that we can restrict our search to domains where the post-Newtonian expansion is a valid approximation. We define the post-Newtonian expansion not to be valid if any of the following conditions occur:
We can further require as a matter of convention that the lowest-order nonzero coefficient in the frequency expansion be positive; this is simply a sign convention stating that the frequency of a system be positive at large radii.
The first two conditions above allow us to set firm limits on the range of the initial \(x_\mathrm{start}\). Let \(C_j\) be the lowest-order nonzero coefficient; then for every nonzero \(C_{k>j}\) we can define a point \(x_k=|C_j/C_k|^{1/(k-j)}\) where that term exceeds the leading-order term in magnitude. We can therefore limit the range of \(x\) to values less than \(x_\mathrm{max}\), which is the minimum of \(\sqrt{2}\) and all \(x_k\). We note that even if we were to extend the post-Newtonian expansion in Eq. \eqref{eq_ppn_fnorm} to an infinite number of terms, this definition of \(x_\mathrm{max}\) implies that the frequency is guaranteed to be monotonic up to \(x_\mathrm{max}(5-\sqrt{7})/6\), and positive up to \(x_\mathrm{max}/2\). Thus we can confidently begin our search for \(x_\mathrm{start}\) in the domain \((0,0.39x_\mathrm{max})\), where the leading-order term dominates, and end it if we ever exceed \(x_\mathrm{max}\).
We therefore bracket our value of \(x_\mathrm{start}\) as follows: We start with an initial guess \(x_\mathrm{guess}=(y_\mathrm{start}/C_j)^{1/(j+3)}\), or \(0.39x_\mathrm{max}\), whichever is less. If \(y(x_\mathrm{guess})>y_\mathrm{start}\), we iteratively decrease \(x\) by factors of 0.95 until \(y(x)<y_\mathrm{start}\); this is guaranteed to occur within a few iterations, since we are moving into a regime where the leading-order behaviour dominates more and more. If \(y(x_\mathrm{guess})<y_\mathrm{start}\), we iteratively increase \(x\) by factors of 1.05 until \(y(x)>y_\mathrm{start}\), or until \(x>x_\mathrm{max}\); this is also guaranteed to occur quickly because, in the worst case, it only takes about 20 iterations to step from \(0.39x_\mathrm{max}\) to \(x_\mathrm{max}\), and if \(x_\mathrm{guess}\) were much lower than \(0.39x_\mathrm{max}\) it would have been a pretty good guess to begin with. If at any point while increasing \(x\) we find that \(y\) is decreasing, we determine that the starting frequency is already in a regime where the post-Newtonian approximation is invalid, and we return an error. Otherwise, once we have bracketed the value of \(x_\mathrm{start}\), we use LALSBisectionFindRoot() to pin down the value to an accuracy of a part in \(10^6\).
Computing the phase and amplitudes:
Once we have \(x_\mathrm{start}\), we can find \(\Theta_\mathrm{start}\), and begin incrementing it; at each timestep we compute \(x\) and hence \(f\), \(\phi\), \(A_+\), and \(A_\times\) according to Eq. \eqref{eq_ppn_freq}, Eq. \eqref{eq_ppn_phi}, Eq. \eqref{eq_ppn_aplus}, and Eq. \eqref{eq_ppn_across}. The routine progressively creates a list of length-1024 arrays and fills them. The process stops when any of the following occurs:
*params. In the last case an error is returned; otherwise the waveform is deemed "complete". Output arrays are created of the appropriate length and are filled with the data.
Internally, the routine keeps a list of all coefficients, as well as a list of booleans indicating which terms are nonzero. The latter allows the code to avoid lengthy floating-point operations (especially the logarithm in the post \({}^{5/2}\)-Newtonian phase term) when these are not required.
When generating the waveform, we note that the sign of \(\dot{f}\) is the same as the sign of \(y'/x^2 = \sum (k+3)C_k x^k\), and use this series to test whether the frequency has stopped increasing. (The reason is that for waveforms far from coalescence \(f\) is nearly constant: numerical errors can cause positive and negative fluctuations \(\Delta f\) bewteen timesteps. The analytic formulae for \(\dot{f}\) or \(y'\) are less susceptible to this.) The coefficients \((k+3)C_k\) are also precomputed for added efficiency.
Warnings and suggestions:
If no post-Newtonian parameters are provided (i.e.\ params->ppn=NULL), we generate a post \({}^2\)-Newtonian waveform, not a post \({}^{5/2}\)-Newtonian waveform. This is done not only for computationally efficiency, but also because the accuracy and reliability of the post \({}^{5/2}\)-Newtonian waveform is actually worse. You can of course specify a post \({}^{5/2}\)-Newtonian waveform with an appropriate assignment of params->ppn, but you do so at your own risk!
This routine also performs no sanity checking on the requested sampling interval \(\Delta t=\)params->deltaT, because this depends very much on how one intends to use the generated waveform. If you plan to generate actual wave functions \(h_{+,\times}(t)\) at the same sample rate, then you will generally want a sampling interval \(\Delta t<1/2f_\mathrm{max}\); you can enforce this by specifying a suitable params->fStopIn.
However, many routines (such as those in SimulateCoherentGW.h) generate actual wave functions by linear interpolation of the amplitude and phase data, which then need only be sampled on timescales \(\sim\dot{f}^{-1/2}\) rather than \(\sim f^{-1}\). More precisely, we would like our interpolated phase to differ from the actual phase by no more than some specified amount, say \(\pi/2\) radians. The largest deviation from linear phase evolution will typically be on the order of \(\Delta\phi\approx(1/2)\ddot{\phi}(\Delta t/2)^2\approx(\pi/4)\Delta f\Delta t\), where \(\Delta f\) is the frequency shift over the timestep. Thus in general we would like to have
\[ \Delta f \Delta t \lesssim 2 \]
for our linear interpolation to be valid. This routine helps out by setting the output parameter field params->dfdt equal to the maximum value of \(\Delta f\Delta t\) encountered during the integration.
| stat | UNDOCUMENTED |
| output | UNDOCUMENTED |
| params | UNDOCUMENTED |
Definition at line 341 of file GeneratePPNInspiral.c.
| void LALGeneratePPNAmpCorInspiral | ( | LALStatus * | stat, |
| CoherentGW * | output, | ||
| PPNParamStruc * | params | ||
| ) |
Computes a parametrized post-Newtonian inspiral waveform with ampltidude corrections.
Phase computed to 3.5PN. Amplitude computed to 2.5PN.
The ampitude corrrected gravitaional wave polarizations \(h_+\) and \(h_x\) are stored in output.h.
Warning! output.a is used to store the first three harmonics in alternate values, (i.e. [a1(0),a2(0),a3(0),a1(dt),a2(dt),a3(dt)...]) as this will be used for filtering with higher harmonic waveforms.
Although \(h_{+,\times} are computed, \)output.phi is also required for filtering.
| stat | UNDOCUMENTED |
| output | UNDOCUMENTED |
| params | UNDOCUMENTED |
Definition at line 182 of file GeneratePPNAmpCorInspiral.c.
| void LALGetInspiralParams | ( | LALStatus * | stat, |
| PPNParamStruc * | output, | ||
| GalacticInspiralParamStruc * | input, | ||
| RandomParams * | params | ||
| ) |
Computes the input parameters for a PPN inspiral.
This function takes a Galactic location and pair of masses from *input and uses them to set the PPNParamStruc fields output->position, output->mTot, output->eta, and output->d. The fields output->psi, output->inc, and output->phi are set randomly to reflect a uniform distribution in solid angle (that is, cosine of inclination is uniform between \(-1\) and 1, other angles are uniform between 0 and \(2\pi\)). The routine uses the random sequence specified by *params when given, but if *params=NULL a new sequence is started internally using the current execution time as a seed. The field input->geocentEndTime is ignored by this routine.
The other PPNParamStruc input fields are not touched by this routine, and must be specified externally before generating a waveform with this structure.
Galactocentric Galactic axial coordinates \(\rho\), \(z\), and \(l_G\) are transformed to geocentric Galactic Cartesian coordinates:
\begin{eqnarray} x_e & = & R_e + \rho\cos l_G \;,\\ y_e & = & \rho\sin l_G \;,\\ z_e & = & z \;, \end{eqnarray}
where
\[ R_e \approx 8.5\,\mathrm{kpc} \]
is the distance to the Galactic core (this constant will probably migrate into LALConstants.h eventually). These are converted to geocentric Galactic spherical coordinates:
\begin{eqnarray} d & = & \sqrt{x_e^2 + y_e^2 + z_e^2} \;,\\ b & = & \arcsin\left(\frac{z_e}{d_e}\right) \;,\\ l & = & \arctan\!2(y_e,x_e) \;. \end{eqnarray}
In the calculation of \(d\) we factor out the leading order term from the square root to avoid inadvertent overflow, and check for underflow in case the location lies on top of the Earth. The angular coordinates are then transformed to equatorial celestial coordinates \(\alpha\) and \(\delta\) using the routines in SkyCoordinates.h.
Definition at line 80 of file GetInspiralParams.c.
| #define GENERATEPPNINSPIRALH_ENUL 1 |
Unexpected null pointer in arguments.
Definition at line 168 of file GeneratePPNInspiral.h.
| #define GENERATEPPNINSPIRALH_EOUT 2 |
output field a, f, phi, or shift already exists
Definition at line 169 of file GeneratePPNInspiral.h.
| #define GENERATEPPNINSPIRALH_ETBAD 3 |
Bad sampling interval.
Definition at line 170 of file GeneratePPNInspiral.h.
| #define GENERATEPPNINSPIRALH_EFBAD 4 |
Bad starting frequency; could not get valid start time.
Definition at line 171 of file GeneratePPNInspiral.h.
| #define GENERATEPPNINSPIRALH_EPBAD 5 |
Bad post-Newtonian parameters.
Definition at line 172 of file GeneratePPNInspiral.h.
| #define GENERATEPPNINSPIRALH_EMBAD 6 |
Bad masses.
Definition at line 173 of file GeneratePPNInspiral.h.
| #define GENERATEPPNINSPIRALH_EDBAD 7 |
Bad distance.
Definition at line 174 of file GeneratePPNInspiral.h.
| #define GENERATEPPNINSPIRALH_EMEM 8 |
Out of memory.
Definition at line 175 of file GeneratePPNInspiral.h.
| #define GENERATEPPNINSPIRALH_EFSTOP 0 |
Reached requested termination frequency.
Definition at line 202 of file GeneratePPNInspiral.h.
| #define GENERATEPPNINSPIRALH_ELENGTH 1 |
Reached maximum length, or end of provided time series vector.
Definition at line 203 of file GeneratePPNInspiral.h.
| #define GENERATEPPNINSPIRALH_EFNOTMON 2 |
Frequency no longer increasing monotonically.
Definition at line 204 of file GeneratePPNInspiral.h.
| #define GENERATEPPNINSPIRALH_EPNFAIL 3 |
Evolution dominated by higher-order PN terms.
Definition at line 205 of file GeneratePPNInspiral.h.
| #define GENERATEPPNINSPIRALH_ERTOOSMALL 4 |
Orbital radius too small for PN approximation.
Definition at line 206 of file GeneratePPNInspiral.h.
1.9.4