Note

This tutorial was generated from a Jupyter notebook that can be downloaded here. If you’d like to reproduce the results in the notebook, or make changes to the code, we recommend downloading this notebook and running it with Jupyter as certain cells (mostly those that change plot styles) are excluded from the tutorials.

Derivations and Equation Reference

This guide explains the origin and derivation of the equations used in LEGWORK functions. This follows closely Section 3 of the LEGWORK paper but with some extra steps shown!

At the end of this document (here) is a table that relates each of the functions in LEGWORK to an equation in this document.

Conversions and Definitions (utils)

This section contains a miscellaneous collection of conversions and definitions that are useful in the later derivations. First, the chirp mass of a binary is defined as

\begin{equation} \mathcal{M}_{c} = \frac{(m_1 m_2)^{3/5}}{(m_1 + m_2)^{1/5}}, \label{eq:chirpmass} \end{equation}

where \(m_1\) and \(m_2\) are the primary and secondary mass of the binary. This term often shows up in many equations and hence is easier to measure in gravitational wave data analysis than the individual component masses.

Kepler’s third law allows one to convert between orbital frequency, \(f_{\rm orb}\), and the semi-major axis, \(a\), of a binary. For convenience we show it here

\begin{equation} a = \left(\frac{G(m_1 + m_2)}{(2 \pi f_{\rm orb})^{2}}\right)^{1/3},\qquad f_{\rm orb} = \frac{1}{2 \pi} \sqrt{\frac{G(m_1 + m_2)}{a^3}}. \label{eq:kepler3rd} \end{equation}

As we deal with eccentric binaries, the different harmonic frequencies of gravitational wave emission become important. We can write that the relative power radiated into the \(n^{\rm th}\) harmonic for a binary with eccentricity \(e\) is [PetersMathews63] (Eq. 20)

\begin{equation} \begin{aligned} g(n, e) = \frac{n^{4}}{32} & \left\{ \left[ J_{n-2}(n e)-2 e J_{n-1}(n e)+\frac{2}{n} J_{n}(n e)+2 e J_{n+1}(n e)-J_{n+2}(n e)\right]^{2}\right.\\ &\left.+\left(1-e^{2}\right)\left[J_{n-2}(n e)-2 J_{n}(n e)+J_{n+2}(n e)\right]^{2}+\frac{4}{3 n^{2}}\left[J_{n}(n e)\right]^{2}\right\}, \end{aligned} \label{eq:g(n,e)} \end{equation}

where \(J_{n}(v)\) is the Bessel function of the first kind. Thus, the sum of \(g(n, e)\) over all harmonics gives the factor by which the gravitational wave emission is stronger for a binary of eccentricity \(e\) over an otherwise identical circular binary. This enhancement factor is defined by Peters as [PetersMathews63] (Eq. 17)

\begin{equation} F(e) = \sum_{n = 1}^{\infty} g(n, e) = \frac{1 + (73 / 24) e^2 + (37 / 96) e^4}{(1 - e^2)^{7/2}}. \label{eq:eccentricity_enhancement_factor} \end{equation}

Note that \(F(0) = 1\) as one would expect. A useful number to remember is that \(F(0.5) \approx 5.0\), or in words, a binary with eccentricity \(0.5\) loses energy to gravitational waves at a rate about \(5\) times higher than a similar circular binary.

Binary Evolution (evol)

Circular binaries

For a circular binary, the evolution can be calculated analytically as the rate at which the binary shrinks can be readily integrated [Peters64] (Eq. 5.6)

\begin{equation} \frac{\mathrm{d}a}{\mathrm{d}t}_{e = 0} = - \frac{\beta}{a^3}, \end{equation}

where \(\beta\) is

\begin{equation} \beta(m_1, m_2) = \frac{64}{5} \frac{G^3}{c^5} m_1 m_2 (m_1 + m_2). \label{eq:beta_peters} \end{equation}

Integratig this gives the semi-major axis of a circular binary as a function of time as [Peters64] (Eq. 5.9)

\begin{equation} a(t, m_1, m_2) = [a_0^4 - 4 t \beta(m_1, m_2)]^{1/4}, \label{eq:a_over_time_circ} \end{equation}

where \(a_0\) is the initial semi-major axis. Moreover, we can set the final semi-major axis in Eq. \eqref{eq:a_over_time_circ} equal to zero and solve for the inspiral time ([Peters64] Eq. 5.10)

\begin{equation} t_{\rm merge, circ} = \frac{a_0^4}{4 \beta} \label{eq:t_merge_circular} \end{equation}

Eccentric binaries

Eccentric binaries are more complicated because the semi-major axis and eccentricity both evolve simultaneously and depend on one another. These equations cannot be solved analytically and require numerical integration. Firstly, we can relate \(a\) and \(e\) with [Peters64] (Eq. 5.11)

\begin{equation} a(e) = c_0 \frac{e^{12/19}}{(1 - e^2)} \left(1 + \frac{121}{304} e^2\right)^{870/2299}, \label{eq:a_from_e} \end{equation}

where \(c_0\) is

\begin{equation} c_0(a_0, e_0) = a_0 \frac{(1 - e_0^2)}{e_0^{12/19}} \left(1 + \frac{121}{304} e_0^2\right)^{-870/2299} \label{eq:c0_peters} \end{equation}

where \(a_0\) and \(e_0\) are the initial semi-major axis and eccentricity respectively.

where \(c_0\) is defined in Eq. \eqref{eq:c0_peters} – such that the initial conditions are satisfied. Then we can numerically integrate [Peters64] (Eq. 5.13)

\begin{equation} \frac{\mathrm{d}e}{\mathrm{d}t} = -\frac{19}{12} \frac{\beta}{c_{0}^{4}} \frac{e^{-29 / 19}\left(1-e^{2}\right)^{3 / 2}}{\left[1+(121 / 304) e^{2}\right]^{1181 / 2299}}, \label{eq:dedt} \end{equation}

to find \(e(t)\) and use this in conjunction with Eq. \eqref{eq:a_from_e} to solve for \(a(t)\), which can in turn be converted to \(f_{\rm orb}(t)\).

Furthermore, we can invert this to find the inspiral time by using that \(e \to 0\) when the binary merges which gives [Peters64] (Eq. 5.14)

\begin{equation} t_{\rm merge} = \frac{12}{19} \frac{c_{0}^{4}}{\beta} \int_0^{e_0} \frac{\left[1+(121 / 304) e^{2}\right]^{1181 / 2299}}{e^{-29 / 19}\left(1-e^{2}\right)^{3 / 2}} \mathrm{d}e \label{eq:t_merge_eccentric} \end{equation}

For very small or very large eccentricities we can approximate this integral using the following expressions (given in unlabelled equations after [Peters64] Eq. 5.14)

\begin{equation} t_{\rm merge,\, e^2 \ll 1} = \frac{c_0^4}{4 \beta} \cdot e_0^{48 / 19} \end{equation}

\begin{equation} t_{\rm merge,\, (1 - e^2) \ll 1} = \frac{768}{425} \frac{a_0^4}{4 \beta} (1 - e_0^2)^{7/2} \end{equation}

The standard threshold employed by LEGWORK for small eccentricities is e = 0.15 and for large eccentricities is e = 0.9999 (as this approximates \(t_{\rm merge}\) with an error below roughly 2%), though we note that this can be customised by the user if desired.

In addition, we implement the fit to Eq. 13 from Mandel (2021) that approximates the merger time as

\begin{equation} t_{\rm merge} \approx t_{\rm merge, circ} (1 - e_0^2)^{7/2} (1 + 0.27 e_0^{10} + 0.33 e_0^{20} + 0.2 e_0^{1000}) \end{equation}

which gives \(t_{\rm merge}\) with an error below 3% for eccentricities below 0.9999. We additionally add a rudimentary polynomial fit to further reduce this error to below 0.5%. The user may specify whether to use this fit or perform the full integral when calculating merger times in LEGWORK.

Gravitational Wave Strains (strain)

Characteristic Strain

The strength of a gravitational wave in a detector at any one moment is determined by the strain amplitude, \(h_0\). However, for stellar-origin sources at mHz frequencies, the signal can be present in the detector for many years. This means that, the \(n^{\rm th}\) harmonic of the binary will spend approximately \(f_n / \dot{f}_n\) seconds (or \(f_n^2 / \dot{f}_n\) cycles) in the vicinity of a frequency \(f_n\) [FinnThorne00]. This leads to the signal ‘accumulating’ at the frequency \(f_n\).

Therefore, to account for the integration of the signal over the mission, we instead use the ‘characteristic’ strain amplitude of the \(n^{\rm th}\) harmonic, \(h_{c, n}\), which is the term present in the general signal-to-noise ratio equation. This can be related the to the strain amplitude in the \(n^{\rm th}\) harmonic, \(h_{0, n}\), as (e.g. [FinnThorne00] [MooreColeBerry15])

\begin{equation} h_{c, n}^2 = \left(\frac{f_n^2}{\dot{f}_n} \right) h_n^2. \label{eq:strain-charstrain} \end{equation}

Note

Note that this is factor of 2 different from [FinnThorne00]. This is because the factor of 2 is already included in the [RobsonCornishLiu19] sensitivity curve and so we remove it here.

The characteristic strain represents the strain measured by the detector over the duration of the mission (approximated as a single broad-band burst), whilst the strain amplitude is the strength of the GW emission at each instantaneous moment. For a stellar mass binary, the characteristic strain in the \(n^{\rm th}\) harmonic is given by (e.g. [BarackCutler04] Eq. 56; [FlanaganHughes98] Eq. 5.1)

\begin{equation} h_{c,n}^2 = \frac{1}{(\pi D_L)^2} \left( \frac{2 G}{c^3} \frac{\dot{E_n}}{\dot{f_n}} \right), \label{eq:char_strain_dedf} \end{equation}

where \(D_L\) is the luminosity distance to the binary, \(\dot{E}_n\) is the power radiated in the \(n^{\rm th}\) harmonic and \(\dot{f}_n\) is the rate of change of the \(n^{\rm th}\) harmonic frequency. The power radiated in the \(n^{\rm th}\) harmonic is given by [PetersMathews63] (Eq. 19)

\begin{equation} \dot{E}_n = \frac{32}{5} \frac{G^{4}}{c^5} \frac{m_{1}^{2} m_{2}^{2}\left(m_{1}+m_{2}\right)}{a^{5}} g(n, e), \label{eq:edot_peters} \end{equation}

where \(m_1\) is the primary mass, \(m_2\) is the secondary mass, \(a\) is the semi-major axis of the binary and \(e\) is the eccentricity. Using Eq. \eqref{eq:chirpmass} and Eq. \eqref{eq:kepler3rd}, we can recast Eq. \eqref{eq:edot_peters} in a form more applicable for gravitational wave detections that is a function of only the chirp mass, orbital frequency and eccentricity.

\begin{align} \dot{E}_n &= \frac{32}{5} \frac{G^{4}}{c^5} \left(m_{1}^{2} m_{2}^{2}\left(m_{1}+m_{2}\right)\right) g(n, e) \cdot \left(\frac{(2 \pi f_{\rm orb})^{2}}{G(m_1 + m_2)}\right)^{5/3} \\ \dot{E}_n &= \frac{32}{5} \frac{G^{7/3}}{c^5} \frac{m_{1}^{2} m_{2}^{2}}{\left(m_{1}+m_{2}\right)^{2/3}} (2 \pi f_{\rm orb})^{10/3} g(n, e) \\ \dot{E}_n(\mathcal{M}_c, f_{\rm orb}, e) &= \frac{32}{5} \frac{G^{7 / 3}}{c^{5}}\left(2 \pi f_{\mathrm{orb}} \mathcal{M}_{c}\right)^{10 / 3} g(n, e) \label{eq:edot} \end{align}

The last term needed to define the characteristic strain in Eq. \eqref{eq:char_strain_dedf} is the rate of change of the \(n^{\rm th}\) harmonic frequency. We can first apply the chain rule and note that

\begin{equation} \dot{f}_{n} = \frac{\mathrm{d}f_{n}}{\mathrm{d} a} \frac{\mathrm{d} a}{\mathrm{d} t}. \label{eq:fdot_chainrule} \end{equation}

The frequency of the \(n^{\rm th}\) harmonic is simply defined as \(f_n = n \cdot f_{\rm orb}\) and therefore we can find an expression for \(\mathrm{d} {f_{n}} / \mathrm{d} {a}\) by rearranging and differentiating Eq. \eqref{eq:kepler3rd}

\begin{align} f_{n} &= \frac{n}{2 \pi} \sqrt{\frac{G(m_1 + m_2)}{a^3}}, \\ \frac{\mathrm{d}f_{n}}{\mathrm{d} a} &= -\frac{3 n}{4 \pi} \frac{\sqrt{G(m_1 + m_2)}}{a^{5/2}}. \label{eq:dfda} \end{align}

The rate at which the semi-major axis decreases is [Peters64] (Eq. 5.6)

\begin{equation} \frac{\mathrm{d} a}{\mathrm{d} t} = -\frac{64}{5} \frac{G^{3} m_{1} m_{2}\left(m_{1}+m_{2}\right)}{c^{5} a^{3}} F(e). \label{eq:dadt} \end{equation}

Substituting Eq. \eqref{eq:dfda} and Eq. \eqref{eq:dadt} into Eq. \eqref{eq:fdot_chainrule} gives an expression for \(\dot{f}_{n}\)

\begin{align} \dot{f}_n &= -\frac{3 n}{4 \pi} \frac{\sqrt{G(m_1 + m_2)}}{a^{5/2}} \cdot -\frac{64}{5} \frac{G^{3} m_{1} m_{2}\left(m_{1}+m_{2}\right)}{c^{5} a^{3}} F(e), \\ \dot{f}_n &= \frac{48 n}{5 \pi} \frac{G^{7/2}}{c^5} \left(m_1 m_2 (m_1 + m_2)^{3/2}\right) \frac{F(e)}{a^{11/2}}, \end{align}

which, as above with \(\dot{E}_n\), we can recast using Kepler’s third law and the definition of the chirp mass

\begin{align} \dot{f}_n &= \frac{48 n}{5 \pi} \frac{G^{7/2}}{c^5} \left(m_1 m_2 (m_1 + m_2)^{3/2}\right) F(e) \cdot \left(\frac{(2 \pi f_{\rm orb})^{2}}{G(m_1 + m_2)}\right)^{11/6}, \\ &= \frac{48 n}{5 \pi} \frac{G^{5/3}}{c^5} \frac{m_1 m_2}{(m_1 + m_2)^{1/3}} \cdot (2 \pi f_{\rm orb})^{11/3} \cdot F(e), \\ \dot{f}_n(\mathcal{M}_c, f_{\rm orb}, e) &= \frac{48 n}{5 \pi} \frac{\left(G \mathcal{M}_c \right)^{5/3}}{c^5} (2 \pi f_{\rm orb})^{11/3} F(e) \label{eq:fdot} \end{align}

With definitions of both \(\dot{E}_n\) and \(\dot{f}_n\), we are now in a position to find an expression for the characteristic strain by plugging Eq. \eqref{eq:edot} and Eq. \eqref{eq:fdot} into Eq. \eqref{eq:char_strain}:

\begin{align} h^2_{c,n} &= \frac{1}{(\pi D_L)^2} \left( \frac{2 G}{c^3} \frac{\frac{32}{5} \frac{G^{7 / 3}}{c^{5}}\left(2 \pi f_{\mathrm{orb}} \mathcal{M}_{c}\right)^{10 / 3} g(n, e)}{\frac{48 n}{5 \pi} \frac{\left(G \mathcal{M}_c\right)^{5/3}}{c^5} (2 \pi f_{\rm orb})^{11/3} F(e)} \right) \\ &= \frac{1}{(\pi D_L)^2} \left( \frac{2^{5/3} \pi^{2/3}}{3} \frac{(G \mathcal{M}_c)^{5/3}}{c^3} \frac{1}{f_{\rm orb}^{1/3}} \frac{g(n, e)}{n F(e)} \right) \end{align}

This gives a final simplified expression for the characteristic strain amplitude of a GW source.

\begin{equation} h_{c,n}^2(\mathcal{M}_c, D_L, f_{\rm orb}, e) = \frac{2^{5/3}}{3 \pi^{4/3}} \frac{(G \mathcal{M}_c)^{5/3}}{c^3 D_L^2} \frac{1}{f_{\rm orb}^{1/3}} \frac{g(n, e)}{n F(e)} \label{eq:char_strain} \end{equation}

Using the conversion between characteristic strain and regular strain, in addition to Eq. \eqref{eq:fdot} and Eq. \eqref{eq:char_strain}, we can write an expression for the strain amplitude of gravitational waves in the \(n^{\rm th}\) harmonic

\begin{align} h_n^2 &= \left(\frac{\dot{f}_n}{f_n^2} \right) h_{c, n}^2, \\ h_n^2 &= \left(\frac{48 n}{5 \pi} \frac{\left(G \mathcal{M}_c \right)^{5/3}}{c^5} F(e) \cdot \frac{(2 \pi f_{\rm orb})^{11/3}}{n^2 f_{\rm orb}^2} \right) \left(\frac{2^{5/3}}{3 \pi^{4/3}} \frac{(G \mathcal{M}_c)^{5/3}}{c^3 D_L^2} \frac{1}{n f_{\rm orb}^{1/3}} \frac{g(n, e)}{F(e)} \right), \end{align}

This gives a final simplified expression for the strain amplitude of a GW source.

\begin{equation} h_n^2(\mathcal{M}_c, f_{\rm orb}, D_L, e) = \frac{2^{28/3}}{5} \frac{(G \mathcal{M}_c)^{10/3}}{c^8 D_L^2} \frac{g(n, e)}{n^2} \left(\pi f_{\rm orb} \right)^{4/3} \label{eq:strain} \end{equation}

Amplitude modulation for orbit averaged sources

Because the LISA detectors are not stationary and instead follow an Earth-trailing orbit, the antenna pattern of LISA is not isotropically distributed or stationary. For sources that have a known position, inclination, and polarisation, we can consider the amplitude modulation of the strain due to the average motion of LISA’s orbit. We follow the results of [CornishLarson03] to define the amplitude modulation. However, we follow the conventions of more recent papers (e.g. [BabakHewitsonPetiteau21], Eq.67) to write the amplitude modulation as

\begin{equation} A_{\rm{mod}}^{2}= \frac{1}{4} \left(1+\cos ^{2} \iota\right)^{2}\left\langle F_{+}^{2}\right\rangle+ \cos ^{2} \iota\left\langle F_{\times}^{2}\right\rangle, \label{eq:amp_mod} \end{equation}

where \(\left\langle F_{+}^{2}\right\rangle\) and \(\left\langle F_{\times}^{2}\right\rangle\), the orbit-averaged detector responses, are defined as

\begin{equation} \left \langle F_{+}^{2} \right\rangle = \frac{1}{4}\big(\cos ^{2} 2 \psi\left\langle D_{+}^{2}\right\rangle -\sin 4 \psi\left\langle D_{+} D_{\times}\right\rangle +\sin ^{2} 2 \psi\left\langle D_{\times}^{2}\right\rangle\big), \label{eq:response_fplus} \end{equation}

\begin{equation} \left\langle F_{\times}^{2} \right\rangle = \frac{1}{4}\big(\cos^{2} 2 \psi \left \langle D_{\times}^{2} \right\rangle +\sin 4 \psi \left \langle D_{+} D_{\times} \right\rangle +\sin ^{2} 2 \psi \left \langle D_{+}^{2} \right\rangle \big), \label{eq:response_fcross} \end{equation}

and

\begin{equation} \left\langle D_{+} D_{\times} \right\rangle = \frac{243}{512} \cos \theta \sin 2 \phi \left(2 \cos ^{2} \phi-1\right) \left(1+\cos ^{2} \theta\right), \label{eq:d_plus_cross} \end{equation}\begin{equation} \left\langle D_{\times}^{2} \right\rangle = \frac{3}{512}\big(120 \sin ^{2} \theta +\cos ^{2} \theta + 162 \sin ^{2} 2 \phi \cos ^{2} \theta\big), \label{eq:d_cross} \end{equation}\begin{equation} \left\langle D_{+}^{2} \right\rangle = \frac{3}{2048}\big[487+158 \cos ^{2} \theta+7 \cos ^{4} \theta -162 \sin ^{2} 2 \phi\left(1+\cos ^{2} \theta\right)^{2}\big]. \label{eq:d_plus} \end{equation}

In the equations above, the inclination is given by \(\iota\), the right ascension and declination are given by \(\phi\) and \(\theta\) respsectively, and the polarisation is given by \(\psi\).

The orbital motion of LISA smears the source frequency by roughly \(10^{-4}\,\rm{mHz}\) due to the antenna pattern changing as the detector orbits, the Doppler shift from the motion, and the phase modulation from the \(+\) and \(\times\) polarisations in the antenna pattern. Generally, the modulation reduces the strain amplitude because the smearing in frequency reduces the amount of signal build up at the true source frequency.

We note that since the orbit averaging is carried out in Fourier space, this requires the frequency to be monochromatic and thus is only implemented in LEGWORK for quasi-circular binaries. We also note that since the majority of the calculations in LEGWORK are carried out for the full position, polarisation, and inclination averages, we place a pre-factor of \(5/4\) on the amplitude modulation in the software implementation to undo the factor of \(4/5\) which arises from the averaging of Equations \eqref{eq:response_fplus} and \eqref{eq:response_fcross}.

Sensitivity Curves (psd)

LISA

For the LISA sensitivity curve, we follow the equations from [RobsonCornishLiu19], which we list here for your convenience.

The effective LISA noise power spectral density is defined as ([RobsonCornishLiu19] Eq. 2)

\begin{equation} S_{\rm n}(f) = \frac{P_n(f)}{\mathcal{R}(f)} + S_c(f), \end{equation}

where \(P_{\rm n}(f)\) is the power spectral density of the detector noise and \(\mathcal{R}(f)\) is the sky and polarisation averaged signal response function of the instrument. Alternatively if we expand out \(P_n(f)\), approximate \(\mathcal{R}(f)\) and simplify we find ([RobsonCornishLiu19] Eq. 1)

\begin{equation} S_{\rm n}(f) = \frac{10}{3 L^2} \left(P_{\rm OMS}(f) + 2 \left( 1 + \cos^2 \left( \frac{f}{f_*} \right) \right)\frac{P_{\rm acc}(f)}{(2 \pi f)^4} \right) \left(1 + \frac{6}{10} \left(\frac{f}{f_*} \right)^2 \right) + S_c(f) \label{eq:LISA_Sn} \end{equation}

where \(L = 2.5\,\mathrm{Gm}\) is detector arm length, \(f^* = 19.09 \, \mathrm{mHz}\) is the response frequency,

\begin{equation} P_{\rm OMS}(f) = \left(1.5 \times 10^{-11} \mathrm{m}\right)^{2}\left(1+\left(\frac{2 \mathrm{mHz}}{f}\right)^{4}\right) \mathrm{Hz}^{-1} \end{equation}

is the single-link optical metrology noise ([RobsonCornishLiu19] Eq. 10),

\begin{equation} P_{\rm acc}(f) = \left(3 \times 10^{-15} \mathrm{ms}^{-2}\right)^{2}\left(1+\left(\frac{0.4 \mathrm{mHz}}{f}\right)^{2}\right)\left(1+\left(\frac{f}{8 \mathrm{mHz}}\right)^{4}\right) \mathrm{Hz}^{-1} \end{equation}

is the single test mass acceleration noise ([RobsonCornishLiu19] Eq. 11) and

\begin{equation} S_{c}(f)=A f^{-7 / 3} e^{-f^{\alpha}+\beta f \sin (\kappa f)}\left[1+\tanh \left(\gamma\left(f_{k}-f\right)\right)\right] \mathrm{Hz}^{-1} \end{equation}

is the galactic confusion noise ([RobsonCornishLiu19] Eq. 14), where the amplitude \(A\) is fixed as \(9 \times 10^{-45}\) and the various parameters change over time:

parameter	6 months	1 year	2 years	4 years
\(\alpha\)	0.133	0.171	0.165	0.138
\(\beta\)	243	292	299	-221
\(\kappa\)	482	1020	611	521
\(\gamma\)	917	1680	1340	1680
\(f_{k}\)	0.00258	0.00215	0.00173	0.00113

TianQin

We additionally allow other instruments than LISA. We have the TianQin sensitivity curve built in where we use the power spectral density given in [HuangHuKorol+20] Eq. 13.

\begin{equation} \begin{split} S_{N}(f) &= \frac{10}{3 L^{2}}\left[\frac{4 S_{a}}{(2 \pi f)^{4}}\left(1+\frac{10^{-4} H z}{f}\right)+S_{x}\right] \\ & \times\left[1+0.6\left(\frac{f}{f_{*}}\right)^{2}\right] \end{split} \label{eq:tianqin} \end{equation}

where \(L = \sqrt{3} \times 10^5\) km} is the arm length, \(S_a = 1 \times 10^{-30}\) m\(^2\) s\(^{-4}\) Hz\(^{-1}\) is the acceleration noise, \(S_x = 1 \times 10^{-24}\) m\(^2\) Hz\(^{-1}\) is the displacement measurement noise and \(f_* = c / 2 \pi L\) is the transfer frequency. Note that Eq. \ref{eq:tianqin} includes an extra factor of \(10 / 3\) compared to Eq.13 of [HuangHuKorol+20]. [HuangHuKorol+20] absorb this factor into the waveform rather than include it in the power spectral density. We include it to match the convention used by [RobsonCornishLiu19] for the LISA sensitivity curve (see the factor of \(10/3\) in Eq. \ref{eq:LISA_Sn}) so that the sensitivity curves can be compared fairly.

Signal-to-Noise Ratios for a 6-link (3-arm) LISA (snr)

Please note that this section draws heavily from [FlanaganHughes98] Section II C. We go through the same derivations here in more detail than in a paper and hopefully help clarify all of the different stages.

Defining the general SNR

In order to calculate the signal to noise ratio for a given source of gravitational waves (GWs) in the LISA detector, we need to consider the following parameters:

• position of the source on the sky: (\(\theta\), \(\phi\))

• direction from the source to the detector: (\(\iota\), \(\beta\))

• orientation of the source, which fixes the polarisation of the GW: \(\psi\)

• the distance from the source to the detector: \(D_L\)

Then, assuming a matched filter analysis of the GW signal \(s(t) + n(t)\) (where \(s(t)\) is the signal and \(n(t)\) is the noise), which relies on knowing the shape of the signal, the signal to noise ratio, \(\rho\), is given in the frequency domain as

\begin{equation} \rho^2(D_L, \theta, \phi, \psi, \iota, \beta) = \frac{\langle s(t)^{\star}s(t)\rangle}{\langle n(t)^{\star}n(t)\rangle} = 2 \int_{-\infty}^{+\infty} \frac{|\tilde{s}(f)|^2}{P_{\rm n}(f)} df = 4 \int_0^{\infty} \frac{|\tilde{s}(f)|^2}{P_{\rm n}(f)} df, \label{eq:snr_general_start} \end{equation}

where \(\tilde{s}(f)\) is the Fourier transform of the signal, \(s(t)\), and \(P_{\rm n}(f)\) is the one sided power spectral density of the noise defined as as \(\langle n(t)^{\star}n(t)\rangle = \int_0^{\infty} \frac{1}{2}P_{\rm n}(f) df\) (c.f. [RobsonCornishLiu19] Eq. 2). Here, \(\tilde{s}(f)\) is implicitly also dependent on \(D_L, \theta, \phi, \psi, \iota,\) and \(\beta\) as

\begin{equation} |\tilde{s}(f)|^2 = |F_+(\theta, \phi, \psi)\tilde{h}_+(t, D_L, \iota, \beta) + F_{\times}(\theta, \phi, \psi)\tilde{h}_{\times}(t, D_L, \iota, \beta)|^2, \label{eq:signal} \end{equation}

where \(F_{+,\times}\) are the ‘plus’ and ‘cross’ antenna patterns of the LISA detector to the ‘plus’ and ‘cross’ strains, \(h_{+,\times}\). Note throughout any parameters discussed with the subscript \(x_{+,\times}\) refers to both \(x_{+}\) and \(x_{\times}\).

In LISA’s case, when averaged over all angles and polarisations, the antenna patterns are orthogonal thus \(\langle F_+ F_{\times}\rangle = 0\). This means we can rewrite Eq. \ref{eq:signal} as

\begin{equation} |\tilde{s}(f)|^2 = |F_+(\theta, \phi, \psi)\tilde{h}_+(t, D_L, \iota, \beta)|^2 + |F_{\times}(\theta, \phi, \psi)\tilde{h}_{\times}(t, D_L, \iota, \beta)|^2, \end{equation}

which can then be applied to Eq. \eqref{eq:snr_general_start} as

\begin{equation} \rho^2(D_L, \theta, \phi, \psi, \iota, \beta) = 4 \int_0^{\infty} \frac{|F_+\tilde{h}_+|^2 + |F_{\times}\tilde{h}_{\times}|^2}{P_{\rm n}(f)} df. \label{eq:snr_general_simpler} \end{equation}

Average over position and polarisation

Now, we can consider averaging over different quantities. In particular, we can average over the sky position and polarisation as

\begin{equation} \label{eq:position_orientation_ave} \langle \rho \rangle^2_{\theta,\phi,\psi} = 4 \int_0^{\infty} df \int \frac{d\Omega_{\theta,\phi}}{4\pi} \int \frac{d\psi}{\pi} \frac{|F_+(\theta,\phi,\psi)\,\tilde{h}_+(\iota,\beta)|^2 + |F_{\times}(\theta,\phi,\psi)\,\tilde{h}_{\times}(\iota,\beta)|^2}{P_{\rm n}(f)}. \end{equation}

From [RobsonCornishLiu19], we can write the position and polarisation average of the signal response function of the instrument, \(\mathcal{R}\), as

\begin{equation} \label{eq:response} \mathcal{R} = \langle F_+F^{\star}_+ \rangle = \langle F_{\times}F^{\star}_{\times} \rangle, \,\,\rm{where}\,\, \langle F_{+,\times}F^{\star}_{+,\times} \rangle = \int \frac{d\Omega_{\theta,\phi}}{4 \pi} \int \frac{d\psi}{\pi} |F_{+,\times}|^2. \end{equation}

Then combining Eq. \eqref{eq:position_orientation_ave} and Eq. \eqref{eq:response}, we then find

\begin{equation} \label{eq:averaged_antenna_simp} \langle \rho \rangle^2_{\theta,\phi,\psi} = 4 \int_0^{\infty} df \mathcal{R}(f)\, \left(\frac{|\tilde{h}_+|^2 + |\tilde{h}_{\times}|^2}{P_{\rm n}(f)}\right) \end{equation}

Note that this is written in [FlanaganHughes98] for the LIGO response function which is \(\mathcal{R} = \langle F_{+,\times} \rangle ^2 = 1/5\).

Average over orientation

Now, we can average over the orientation of the source: \((\iota, \beta)\), noting that the averaging is independent of the distance \(D_L\). Then we can rewrite \(|\tilde{h}_+|^2 + |\tilde{h}_{\times}|^2\) in terms of two functions \(|\tilde{H}_+|^2\) and \(|\tilde{H}_{\times}|^2\), where \(\tilde{h}_{+,\times} = \tilde{H}_{+,\times}/D_L\). Then, averaging over the source direction gives

\begin{equation} \label{eq:averaged_all} \langle \rho \rangle^2_{(\theta,\phi,\psi),(\iota,\beta)} = \frac{4}{D_L^2} \int_0^{\infty} df \mathcal{R}(f)\,\int \frac{d\Omega_{\iota,\beta}}{4 \pi} \frac{|\tilde{H}_+|^2 + |\tilde{H}_{\times}|^2}{P_{\rm n}(f)}, \end{equation}

where we would like to express \(\tilde{H}_{+,\times}(f)^2\) in terms of the energy spectrum of the GW. To do this, we note that the local energy flux of GWs at the detector is given by (e.g. [PressThorne72] Eq. 6)

\begin{equation} \label{eq:energy_flx} \frac{dE}{dAdt} = \frac{1}{16\pi} \frac{c^3}{G} \overline{\left[\left(\frac{dh_{+}}{dt}\right)^2 + \left(\frac{dh_{\times}}{dt}\right)^2\right]}, \end{equation}

where the bar indicates an average over several cycles of the wave which is appropriate for LISA sources. We can transform Eq. \eqref{eq:energy_flx} using Parseval’s theorem, where we can write

\begin{align} \int_{-\infty}^{+\infty}dt\int dA \frac{dE}{dAdt} & = \int_{-\infty}^{+\infty}dt\int dA \frac{1}{16\pi} \frac{c^3}{G}\overline{\left[\left(\frac{dh_{+}}{dt}\right)^2 + \left(\frac{dh_{\times}}{dt}\right)^2\right]} \\ & = \int_{-\infty}^{+\infty}df \int dA \frac{1}{16\pi} \frac{c^3}{G} \Big[\big((-2\pi if)|\tilde{h}_{+}|\big)^2 + \big((-2\pi if)|\tilde{h}_{\times}|\big)^2\Big] \\ & = \int_{-\infty}^{+\infty}df \int dA \frac{1}{16\pi} \frac{c^3}{G} (2\pi f)^2 \left(|\tilde{h}_{+}|^2 + |\tilde{h}_{\times}|^2 \right) \\ & = \int_{-\infty}^{+\infty}df \int dA \frac{c^3}{G} \frac{\pi f^2}{4} \left(|\tilde{h}_{+}|^2 + |\tilde{h}_{\times}|^2 \right) \\ & = \int_{0}^{\infty}df \int dA \frac{c^3}{G} \frac{\pi f^2}{2} \left(|\tilde{h}_{+}|^2 + |\tilde{h}_{\times}|^2 \right). \label{eq:Parseval} \end{align}

Note that we perform a Fourier transform of the square of the time derivatives in the second line. Now, since \(A = D_L^2 \Omega\) and \(|\tilde{h}_{+,\times}|^2 = |\tilde{H}_{+,\times}|^2 / D_L^2\), we know

\begin{equation} \label{eq:little_h_to_big} |\tilde{h}_{+,\times}|^2 dA = |\tilde{H}_{+,\times}|^2 d\Omega_{\iota,\beta}, \end{equation}

then we can write Eq. \eqref{eq:Parseval} in terms of \(|H_{+,\times}|^2\) as

\begin{align} \int_{-\infty}^{+\infty}dt\int dA \frac{dE}{dAdt} & = \int_{0}^{\infty}df \int dA \frac{c^3}{G} \frac{\pi f^2}{2} \left(|\tilde{h}_{+}|^2 + |\tilde{h}_{\times}|^2 \right) \\ & = \int_{0}^{\infty}df \frac{\pi f^2}{2} \frac{c^3}{G} \int d\Omega \left(|\tilde{H}_{+}|^2 + |\tilde{H}_{\times}|^2 \right). \label{eq:Parseval_2} \end{align}

We can note that by using Eq. \eqref{eq:little_h_to_big} and performing a Fourier transform we also have that

\begin{equation} \int_{-\infty}^{+\infty}dt\int dA \frac{dE}{dAdt} = \int_{0}^{\infty}df \int d\Omega \frac{dE}{d\Omega df}. \label{eq:deriv_relation} \end{equation}

From inspection of Eq. \eqref{eq:Parseval_2} and Eq. \eqref{eq:deriv_relation}, we can write the spectral energy flux as

\begin{equation} \int d\Omega \frac{dE}{d\Omega df} = \frac{\pi f^2}{2} \frac{c^3}{G} \int d\Omega \left(|\tilde{H}_{+}|^2 + |\tilde{H}_{\times}|^2 \right) . \label{eq:se_flux} \end{equation}

Fully averaged SNR equation

We are now in a position to write an expression for the fully averaged SNR. Let’s take Eq. \eqref{eq:se_flux} and apply it to Eq. \eqref{eq:averaged_all}

\begin{equation} \langle \rho \rangle^2_{(\theta,\phi,\psi),(\iota,\beta)} = \frac{4}{D_L^2} \frac{G}{c^3} \int_0^{\infty} df \frac{1}{P_{\rm n}(f)/\mathcal{R}(f)} \int \frac{d\Omega}{4\pi} \frac{dE}{d\Omega df} \frac{2}{\pi f^2}. \end{equation}

This simplifies nicely to

\begin{equation} \langle \rho \rangle^2 = \frac{2 G}{(\pi D_L)^2 c^3}\int_0^{\infty}df \frac{dE}{df}\frac{1}{f^2 P_{\rm n}(f)/\mathcal{R}(f)}. \end{equation}

Finally, noting that \(dE/df = dE/dt \times dt/df = \dot{E}/\dot{f}\), we can use the definition of the characteristic strain from Eq. \eqref{eq:char_strain_dedf} to finish up our position, direction, and orientation/polarisation averaged SNR as

\begin{equation} \langle \rho \rangle^2_{(\theta,\phi,\psi),(\iota,\beta)} = \int_0^{\infty}df \frac{h_{c}^2}{f^2P_{\rm n}(f)/\mathcal{R}(f)} = \int_0^{\infty}df \frac{h_{c}^2}{f^2 S_{\rm n}(f)}, \label{eq:snr_finished_circ} \end{equation}

where we have used that the effective power spectral density of the noise is defined as \(S_{\rm n}(f) = P_{\rm n}(f) / \mathcal{R}(f)\). Note that this defintion is the sensitivity for a 6-link (3-arm) LISA detector in the long wavelength limit, which is appropriate for stellar mass binary LISA sources.

It is also important to note that this is only the SNR for a circular binary for which we need only consider the \(n = 2\) harmonic. In the general case, a binary could be eccentric and requires a sum over all harmonics. Thus we can generalise Eq. \eqref{eq:snr_finished_circ} to eccentric binaries with

\begin{equation} \langle \rho \rangle^2_{(\theta,\phi,\psi),(\iota,\beta)} = \sum_{n = 1}^{\infty} \langle \rho_n \rangle^2_{(\theta,\phi,\psi),(\iota,\beta)} = \sum_{n = 1}^{\infty} \int_0^{\infty} d f_n \frac{h_{c, n}^2}{f_n^2 S_{\rm n}(f_n)}, \label{eq:snr_general} \end{equation}

where \(f_n = n \cdot f_{\rm orb}\) (with \(n\) being the harmonic and \(f_{\rm orb}\) the orbital frequency), \(h_{c, n}\) is defined in Eq. \eqref{eq:char_strain} and \(S_{\rm n}\) in Eq. \eqref{eq:LISA_Sn}.

Different SNR approximations

Although Eq. \eqref{eq:snr_general} can be used for every binary, it can be useful to consider different cases in which we can avoid unnecessary sums and integrals. There are four possible cases for binaries in which we can use increasingly simple expressions for the signal-to-noise ratio. Binaries can be circular and stationary in frequency space.

• Circular binaries emit only in the \(n=2\) harmonic and so the sum over harmonics can be removed

• Stationary binaries have \(f_{n, i} \approx f_{n, f}\) and so the small interval allows one to approximate the integral

We refer to non-stationary binaries as ‘evolving’ here though many also use ‘chirping’.

For an evolving and eccentric binary, no approximation can be made and the SNR is found using Eq. \eqref{eq:snr_general}.

For an evolving and circular binary, the sum can be removed and so the SNR found as

\begin{equation} \langle \rho \rangle^2_{\rm c, e} = \frac{1}{4} \int_{f_{2, i}}^{f_{2, f}} \frac{h_{c, 2}^2}{f_2^2 S_{\rm n}(f_2)} \mathrm{d}{f} \label{eq:snr_chirp_circ} \end{equation}

For a stationary and eccentric binary we can approximate the integral.

\begin{align} \langle \rho \rangle^2_{\rm e, s} &= \frac{1}{4} \sum_{n=1}^{\infty} \lim_{\Delta f \to 0} \int_{f_{n}}^{f_{n} + \Delta f_n} \frac{h_{c, n}^2}{f_n^2 S_{\rm n}(f_n)} \mathrm{d}{f_n}, \\ &= \frac{1}{4} \sum_{n=1}^{\infty} \frac{\Delta f_n \cdot h_{c, n}^2}{f_n^2 S_{\rm n}(f_n)}, \\ &= \frac{1}{4} \sum_{n=1}^{\infty} \frac{\dot{f}_n \Delta T \cdot h_{c, n}^2}{f_n^2 S_{\rm n}(f_n)}, \\ &= \frac{1}{4} \sum_{n=1}^{\infty} \left(\frac{\dot{f}_n}{f_n^2} h_{c, n}^2 \right) \frac{T_{\rm obs}}{S_{\rm n}(f_n)}, \\ \langle \rho \rangle^2_{\rm e, s} &= \frac{1}{4} \sum_{n=1}^{\infty} \frac{h_{n}^2 T_{\rm obs}}{S_{\rm n}(f_n)}, \label{eq:snr_stat_ecc} \end{align}

where we have applied Eq. \eqref{eq:strain-charstrain} to convert between strains and labelled \(\Delta t = T_{\rm obs}\). Finally, for a stationary and circular binary the signal-to-noise ratio is simply

\begin{equation} \langle \rho \rangle^2_{\rm c, s} = \frac{1}{4} \frac{h_2^2 T_{\rm obs}}{S_{\rm n}(f_2)} \label{eq:snr_stat_circ} \end{equation}

That’s all for the derivations of the equations! If you are confused by something or think there is a mistake please feel free to open an issue on GitHub.

Continue reading for the function table and references!

Equation to Function Table

The following table gives a list of the functions in the modules and which equation numbers in this document that they come from.

Quantity	Equation	Function
\(\mathcal{M}_c\)	\ref{eq:chirpmass}	legwork.utils.chirp_mass()
\(a\)	\ref{eq:kepler3rd}	legwork.utils.get_a_from_forb()
\(f_{\rm orb}\)	\ref{eq:kepler3rd}	legwork.utils.get_forb_from_a()
\(g(n, e)\)	\ref{eq:g(n,e)}	legwork.utils.peters_g()
\(F(e)\)	\ref{eq:eccentricity_enhancement_factor}	legwork.utils.peters_f()
\(\beta\)	\ref{eq:beta_peters}	legwork.utils.beta()
\(a_{\rm circ}(t), f_{\rm orb, circ}(t)\)	\ref{eq:a_over_time_circ}	legwork.evol.evol_circ()
\(t_{\rm merge, circ}\)	\ref{eq:t_merge_circular}	legwork.evol.get_t_merge_circ()
\(e(t), a(t), f_{\rm orb}(t)\)	\ref{eq:dedt}	legwork.evol.evol_ecc()
\(t_{\rm merge}\)	\ref{eq:t_merge_eccentric}	legwork.evol.get_t_merge_ecc()
\(h_{c,n}\)	\ref{eq:char_strain}	legwork.strain.h_c_n()
\(h_n\)	\ref{eq:strain}	legwork.strain.h_0_n()
\(S_{\rm n}(f)\)	\ref{eq:LISA_Sn}	legwork.psd.power_spectral_density()
\(\rho\)	\ref{eq:snr_general}	legwork.source.Source.get_snr()
\(\rho_{\rm e, e}\)	\ref{eq:snr_general}	legwork.snr.snr_ecc_evolving()
\(\rho_{\rm c, e}\)	\ref{eq:snr_chirp_circ}	legwork.snr.snr_circ_evolving()
\(\rho_{\rm e, s}\)	\ref{eq:snr_stat_ecc}	legwork.snr.snr_ecc_stationary()
\(\rho_{\rm c, s}\)	\ref{eq:snr_stat_circ}	legwork.snr.snr_circ_stationary()

References

[BabakHewitsonPetiteau21]

Stanislav Babak, Martin Hewitson, and Antoine Petiteau. LISA Sensitivity and SNR Calculations. arXiv e-prints, pages arXiv:2108.01167, August 2021. arXiv:2108.01167.

[BarackCutler04]

Leor Barack and Curt Cutler. LISA capture sources: Approximate waveforms, signal-to-noise ratios, and parameter estimation accuracy. Physical Review D, 69(8):082005, April 2004. arXiv:gr-qc/0310125, doi:10.1103/PhysRevD.69.082005.

[CornishLarson03]

Neil J. Cornish and Shane L. Larson. LISA data analysis: Source identification and subtraction. \prd , 67(10):103001, May 2003. arXiv:astro-ph/0301548, doi:10.1103/PhysRevD.67.103001.

[FinnThorne00]

Lee Samuel Finn and Kip S. Thorne. Gravitational waves from a compact star in a circular, inspiral orbit, in the equatorial plane of a massive, spinning black hole, as observed by LISA. Physical Review D, 62(12):124021, December 2000. arXiv:gr-qc/0007074, doi:10.1103/PhysRevD.62.124021.

[FlanaganHughes98]

Éanna É. Flanagan and Scott A. Hughes. Measuring gravitational waves from binary black hole coalescences. I. Signal to noise for inspiral, merger, and ringdown. Physical Review D, 57(8):4535–4565, April 1998. arXiv:gr-qc/9701039, doi:10.1103/PhysRevD.57.4535.

[HuangHuKorol+20]

Shun-Jia Huang, Yi-Ming Hu, Valeriya Korol, Peng-Cheng Li, Zheng-Cheng Liang, Yang Lu, Hai-Tian Wang, Shenghua Yu, and Jianwei Mei. Science with the TianQin Observatory: Preliminary results on Galactic double white dwarf binaries. \prd , 102(6):063021, September 2020. arXiv:2005.07889, doi:10.1103/PhysRevD.102.063021.

[MooreColeBerry15]

C. J. Moore, R. H. Cole, and C. P. L. Berry. Gravitational-wave sensitivity curves. Classical and Quantum Gravity, 32(1):015014, January 2015. arXiv:1408.0740, doi:10.1088/0264-9381/32/1/015014.

[Peters64]

P. C. Peters. Gravitational Radiation and the Motion of Two Point Masses. Physical Review, 136(4B):1224–1232, November 1964. doi:10.1103/PhysRev.136.B1224.

[PetersMathews63]

P. C. Peters and J. Mathews. Gravitational Radiation from Point Masses in a Keplerian Orbit. Physical Review, 131(1):435–440, July 1963. doi:10.1103/PhysRev.131.435.

[PressThorne72]

William H. Press and Kip S. Thorne. Gravitational-Wave Astronomy. \araa , 10:335, January 1972. doi:10.1146/annurev.aa.10.090172.002003.

[RobsonCornishLiu19]

Travis Robson, Neil J. Cornish, and Chang Liu. The construction and use of LISA sensitivity curves. Classical and Quantum Gravity, 36(10):105011, May 2019. arXiv:1803.01944, doi:10.1088/1361-6382/ab1101.