# Derivations and Equation Reference
This guide explains the origin and derivation of the equations used in ``LEGWORK`` functions. Let's go through each of the modules and build up to an equation for the signal-to-noise ratio for a given LISA source.

At the end of this document ([here](#Equation-to-Function-Table)) 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 <cite data-cite="Peters1963"></cite> (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](https://mathworld.wolfram.com/BesselFunctionoftheFirstKind.html). 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 <cite data-cite="Peters1963"></cite> (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.

For binary evolution Peters and Mathews introduced two constants that are useful for calculations though without physical meaning. First, from <cite data-cite="Peters1964"></cite> (Eq. 5.9)

\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}

And additionally from <cite data-cite="Peters1964"></cite> (Eq. 5.11)
\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. 

## 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. This gives the semi-major axis of a circular binary as a function of time as <cite data-cite="Peters1964"></cite> (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 and $\beta$ is defined in Eq. \eqref{eq:beta_peters}. We can use this to also get the frequency evolution by using Kepler's third law Eq. \eqref{eq:kepler3rd}.

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 (<cite data-cite="Peters1964"></cite> 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 <cite data-cite="Peters1964"></cite> (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 defined in Eq. \eqref{eq:c0_peters} -- such that the initial conditions are satisfied. Then we can numerically integrate <cite data-cite="Peters1964"></cite> (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 <cite data-cite="Peters1964"></cite> (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 <cite data-cite="Peters1964"></cite> 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}

## Gravitational Wave Strains (`strain`)

### Characteristic Strain
The characteristic strain from a binary in the $n^{\rm th}$ harmonic is defined as follows (e.g. <cite data-cite="Barack&Cutler2004"></cite> Eq. 56; <cite data-cite="Flanagan+1998"></cite> 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 <cite data-cite="Peters1963"></cite> (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 <cite data-cite="Peters1964"></cite> (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}

### Strain
The strain can be found by dividing the characteristic strain by the square root of twice the number of cycles. This is explained in <cite data-cite="Finn&Thorne2000"></cite> and in further detail in <cite data-cite="Moore+2015"></cite>. The physical reasoning behind the idea is that the binary will spend a certain amount of time in the vicinity of some frequency $f$ and cause a similar gravitational wave strain. This leads to the signal 'accumulating' and resulting in a larger signal-to-noise ratio. Therefore, the characteristic strain represents the strain measured by the detector over the duration of the mission, whilst the strain is what is emitted by the binary at each instantaneous moment.

Following this logic, we can convert between characteristic strain and strain with (e.g. <cite data-cite="Finn&Thorne2000"></cite> see text before Eq. 2.2)

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

Therefore, using this, 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 closely follow the results of <cite data-cite="Cornish2003"></cite> to write down the amplitude modulation as

\begin{equation}
A_{\rm{mod}}^{2}=\frac{1}{2} \left[\left(1+\cos ^{2} \iota\right)^{2}\left\langle F_{+}^{2}\right\rangle+4 \cos ^{2} \iota\left\langle F_{\times}^{2}\right\rangle\right],
\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 <cite data-cite="Robson+2019"></cite>, which we list here for your convenience.

The *effective* LISA noise power spectral density is defined as (<cite data-cite="Robson+2019"></cite> 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 (<cite data-cite="Robson+2019"></cite> Eq. 1)

\begin{equation}
S_{\rm n}(f) = \frac{10}{3 L^2} \left(P_{\rm OMS}(f) + \frac{4 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 (<cite data-cite="Robson+2019"></cite> 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 (<cite data-cite="Robson+2019"></cite> 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 (<cite data-cite="Robson+2019"></cite> Eq. 14), where the amplitude $A$ is fixed as $9 \times 10^{-45}$ and the various parameters change over time:

<center>

| 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 |
</center>

### 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 <cite data-cite="Huang+2020"></cite> 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}

Note that this expression includes an extra factor of 10/3 compared Eq. 13 in <cite data-cite="Huang+2020"></cite>, since <cite data-cite="Huang+2020"></cite> absorbs the factor into the waveform but we instead follow the same convention as <cite data-cite="Robson+2019"></cite> for consistency and include it in this 'effective' PSD function instead.

## Signal-to-Noise Ratios for a 6-link (3-arm) LISA (`snr`)
Please note that this section draws heavily from <cite data-cite="Flanagan+1998"></cite> 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. <cite data-cite="Robson+2019"></cite> 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 <cite data-cite="Robson+2019"></cite>, 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 <cite data-cite="Flanagan+1998"></cite> 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. <cite data-cite="Press&Thorne1972"></cite> Eq. 6)

\begin{equation}
\label{eq:energy_flx}
\frac{dE}{dAdt} = \frac{1}{16\pi} \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} \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} \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} (2\pi f)^2 \left(|\tilde{h}_{+}|^2 + |\tilde{h}_{\times}|^2 \right) \\
& = \int_{-\infty}^{+\infty}df \int dA \frac{\pi f^2}{4} \left(|\tilde{h}_{+}|^2 + |\tilde{h}_{\times}|^2 \right) \\
& = \int_{0}^{\infty}df \int dA \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{\pi f^2}{2} \left(|\tilde{h}_{+}|^2 + |\tilde{h}_{\times}|^2 \right) \\
& = \int_{0}^{\infty}df \frac{\pi f^2}{2} \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} \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} \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}{(\pi D_L)^2}\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} (and use $c=G=1$),

\begin{equation}
h_{c}^2 = \frac{1}{(\pi D_L)^2} \left(\frac{2\dot{E}}{\dot{f}}\right),
\end{equation}

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}
\rho^2_{\rm circ, evol} = \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}
\rho^2_{\rm ecc, stat} &= \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}, \\
&= \sum_{n=1}^{\infty} \frac{\Delta f_n \cdot h_{c, n}^2}{f_n^2 S_{\rm n}(f_n)}, \\
&= \sum_{n=1}^{\infty} \frac{\dot{f}_n \Delta T \cdot h_{c, n}^2}{f_n^2 S_{\rm n}(f_n)}, \\
&= \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)}, \\
\rho^2_{\rm ecc, stat} &= \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}
\rho^2_{\rm circ, stat} = \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](https://github.com/TeamLEGWORK/LEGWORK/issues/new) 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()](../api/legwork.utils.chirp_mass.rst) |
|$a$ | \ref{eq:kepler3rd} | [legwork.utils.get_a_from_forb()](../api/legwork.utils.get_a_from_f_orb.rst) |
|$f_{\rm orb}$ | \ref{eq:kepler3rd} | [legwork.utils.get_forb_from_a()](../api/legwork.utils.get_f_orb_from_a.rst) |
|$g(n, e)$ | \ref{eq:g(n,e)} | [legwork.utils.peters_g()](../api/legwork.utils.peters_g.rst) |
|$F(e)$ | \ref{eq:eccentricity_enhancement_factor} | [legwork.utils.peters_f()](../api/legwork.utils.peters_f.rst) |
|$\beta$ | \ref{eq:beta_peters} | [legwork.utils.beta()](../api/legwork.utils.beta.rst) |
|$a_{\rm circ}(t), f_{\rm orb, circ}(t)$ | \ref{eq:a_over_time_circ} | [legwork.evol.evol_circ()](../api/legwork.evol.evol_circ.rst) |
|$t_{\rm merge, circ}$ | \ref{eq:t_merge_circular} | [legwork.evol.get_t_merge_circ()](../api/legwork.evol.get_t_merge_circ.rst) |
|$e(t), a(t), f_{\rm orb}(t)$ | \ref{eq:dedt} | [legwork.evol.evol_ecc()](../api/legwork.evol.evol_ecc.rst) |
|$t_{\rm merge}$ | \ref{eq:t_merge_eccentric} | [legwork.evol.get_t_merge_ecc()](../api/legwork.evol.get_t_merge_ecc.rst) |
|$h_{c,n}$ | \ref{eq:char_strain} | [legwork.strain.h_c_n()](../api/legwork.strain.h_c_n.rst) |
|$h_n$ | \ref{eq:strain} | [legwork.strain.h_0_n()](../api/legwork.strain.h_0_n.rst) |
|$S_{\rm n}(f)$ | \ref{eq:LISA_Sn} | [legwork.psd.power_spectral_density()](../api/legwork.psd.power_spectral_density.rst) |
|$\rho$ | \ref{eq:snr_general} | [legwork.source.Source.get_snr()](../api/legwork.source.Source.rst#legwork.source.Source.get_snr) |
|$\rho_{\rm e, e}$ | \ref{eq:snr_general}  | [legwork.snr.snr_ecc_evolving()](../api/legwork.snr.snr_ecc_evolving.rst) |
|$\rho_{\rm c, e}$ | \ref{eq:snr_chirp_circ}  | [legwork.snr.snr_circ_evolving()](../api/legwork.snr.snr_circ_evolving.rst) |
|$\rho_{\rm e, s}$ | \ref{eq:snr_stat_ecc} | [legwork.snr.snr_ecc_stationary()](../api/legwork.snr.snr_ecc_stationary.rst) |
|$\rho_{\rm c, s}$ | \ref{eq:snr_stat_circ} | [legwork.snr.snr_circ_stationary()](../api/legwork.snr.snr_circ_stationary.rst) |

## References

07f122d86d44e0733ba5cb153a4d1e39f8a6ee4a
ipynb
Jupyter Notebook
docs/notebooks/Derivations.ipynb
arfon/LEGWORK
91ca299d00ed6892acdf5980f33826421fa348ef
MIT
2021-09-28T21:53:24.000Z
2022-02-05T14:29:44.000Z
docs/notebooks/Derivations.ipynb
arfon/LEGWORK
91ca299d00ed6892acdf5980f33826421fa348ef
MIT
2021-10-31T15:04:26.000Z
2022-03-15T19:01:40.000Z
docs/notebooks/Derivations.ipynb
katiebreivik/LEGWORK
07c3938697ca622fc39d9617d74f28262ac2b1aa
MIT
2021-11-18T09:20:53.000Z
2022-03-16T11:30:44.000Z
Qwen/Qwen-72B
1. YES
2. YES