The construction and use of LISA sensitivity curves

We describe the construction and use of LISA sensitivity curves, the computation of signal-to-noise ratios, and how to plot signal strengths against the sensitivity curve. Figure 1 shows an example of a sensitivity/source plot taken from the LISA L3 mission proposal [1]. The idea, in plotting signal and noise curves in this manner, is that the height a signal is above the sensitivity curve indicates how loud it will be.

Zoom In Zoom Out Reset image size

The literature on this topic can be very confusing, with a profusion of conventions and notation. Unfortunately, some of the choices that have now become standard are misleading, but it is probably too late to change the conventions now. Sensitivity curves are useful for making a quick assessment of what signals may be detectable. While not used for actual data analysis, the sensitivity curve, and signal representations that are shown with them, are designed to represent the quantities that are used in the data analysis. See [2] for a review of gravitational wave sensitivity curves. While useful, the sky-averaged sensitivity can be misleading as there is often significant variation in the sensitivity with sky location. To this end, we also provide expressions and tools for computing the signal-to-noise ratio as a function of sky location, averaged over inclination and polarization angles.

The goal of this paper is to provide a convenient reference that explains how LISA signal-to-noise ratios are calculated using the latest reference design parameters. Python code and a Jupyter notebook for generating the results shown in this paper can be downloaded from GitHub [3]. The Jupyter notebook can be edited and executed directly in an internet browser without the need to install any software using the MyBinder version [4].

The LISA sensitivity curve can be well approximated by the equation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Sn} S_n(\,f) = \frac{10}{3 L^2} \left(P_{\rm OMS}(\,f) + \frac{4P_{\rm acc}(\,f)}{(2\pi f)^4} \right) \left(1 + \frac{6}{10}\left(\frac{f}{f_*}\right)^2 \right) + S_c(\,f), \nonumber \end{align} \tag{ 1 }$$

where L  =  2.5 Gm, $f_*= c/(2\pi L)= 19.09$ mHz, and expressions for $P_{\rm OMS}(\,f)$, $P_{\rm acc}(\,f)$ and S_c(f) are given in equations (10), (11) and (14) below. Here we explain how this curve is computed and how it can be used (and sometimes mis-used).

The simplest type of sensitivity curve, and the one used by the ground-based detector community, is the power spectral density of the detector noise P_n(f), or the amplitude spectral density $\sqrt{P_n(\,f)}$. The mean-squared noise in the frequency band $[\,f_1,f_2]$ is just the integral of P_n(f) over that band. But for a detector like LISA, where signals may have wavelengths that are shorter than the arms of the detector, it is conventional to include the ensuing arm-length penalty in the sensitivity curve [5]. The strain spectral sensitivity is then defined in terms of the square root of the effective noise power spectral density

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S_n(\,f) = \frac{P_n(\,f)}{{{\mathcal R}}(\,f)}, \nonumber \end{align} \tag{ 2 }$$

where ${{\mathcal R}}(\,f)$ is the sky and polarization averaged signal response function of the instrument. The signal response function ${{\mathcal R}}(\,f)$ relates the power spectral density of the incident gravitational wave signals to the power spectral density of the signal recorded in the detector. As such, it might have been more logical to have this multiplying the signal power, but early on [5–8] it was decided to apply the inverse of this factor to the noise power to define a sensitivity curve—c'est la vie.

The response function can be computed by working in the frequency domain, where the gravitational wave amplitude in the detector, $\tilde{h}(\,f)$, is related to the plus and cross gravitational wave amplitudes via

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{h}(\,f) = F^+(\,f) \tilde{h}_+(\,f) + F^\times(\,f) \tilde{h}_\times(\,f) \, , \nonumber \end{align} \tag{ 3 }$$

where $F^+(\theta,\phi,\psi,f)$ and $F^\times(\theta, \phi,\psi,f)$ are the (complex) frequency dependent detector response functions, which depend on the sky location $(\theta,\phi)$ and polarization angle $\psi$ of the source. The sky/polarization averaged spectral power of the signal in the detector, $\langle \tilde{h}(\,f) \tilde{h}^*(\,f)\rangle$ is related to the raw spectral signal power $|\tilde{h}_+(\,f)|^2 +| \tilde{h}_\times(\,f)|^2$ by the response function:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Rdef} & \langle \tilde{h}(\,f) \tilde{h}^*(\,f)\rangle = \langle F^+(\,f) F^{+*}(\,f)\rangle |\tilde{h}_+(\,f)|^2 + \langle F^\times(\,f)F^{\times*}(\,f) \rangle |\tilde{h}_\times(\,f)|^2 \nonumber \\ & \quad = {{\mathcal R}}(\,f)\left(|\tilde{h}_+(\,f)|^2 +| \tilde{h}_\times(\,f)|^2\right) \nonumber \end{align} \tag{ 4 }$$

where ${{\mathcal R}}(\,f) = \langle F^+(\,f) F^{+*}(\,f) \rangle = \langle F^\times(\,f) F^{\times*}(\,f) \rangle$, and the angle brackets indicate the sky/polarization average

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle X \rangle \equiv \frac{1}{4\pi^2} \int_0^\pi {\rm d}\psi \int_0^{2 \pi} {\rm d}\phi \int_0^\pi X \, \sin\theta \, {\rm d}\theta. \nonumber \end{align} \tag{ 5 }$$

For a right-angle interferometer operating in the long wavelength limit, such as LIGO/Virgo, the antenna patterns are real and independent of frequency, and are given by [9]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & F^+ = \frac{1}{2}(1+\cos^2\theta) \cos(2\phi) \cos(2\psi) -\cos\theta \sin2\phi \sin2\psi \nonumber \\ & F^\times = \frac{1}{2}(1+\cos^2\theta) \cos(2\phi) \sin(2\psi) + \cos\theta \sin2\phi \cos2\psi. \nonumber \end{align} \tag{ 6 }$$

For LIGO we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\mathcal R}} = \langle {F^+}^2\rangle = \langle {F^\times}^2\rangle = \frac{1}{32} \int_{-1}^{1} (1+6 x^2 + x^4) {\rm d}x = \frac{1}{5}. \nonumber \end{align} \tag{ 7 }$$

In the LIGO literature this factor is applied to the signals [9, 10] (see for example, equation (2.39) of [10]), leaving the sensitivity curve to be just the power spectral density of the noise. The full expressions for F⁺ (f) and $F^\times(\,f)$ for the Michelson-style interferometry signals for LISA are much more complicated than those for LIGO (they are given in equations (5), (6), (16) and (17) of [11].) For a 3 arm LISA, there are two independent channels for f  <  f_* and three for f  >  f_*, where $f_* = c/(2\pi L)$ is the transfer frequency [12]. The standard convention is to define ${{\mathcal R}}(\,f)$ as being summed over the plus and cross channels (often referred to as the A and E channels following [12]). For sources that have frequency components f  >  f_*, it is more accurate to consider the three-channel expressions given in [12]. The full expression for ${{\mathcal R}}(\,f)$ is not known in closed form, but to leading order is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\mathcal R}}(\,f) = \frac{3}{10} -\frac{507}{5040}\left(\frac{f}{f_*}\right) + \dots . \nonumber \end{align} \tag{ 8 }$$

The first term, 3/10, is a factor of $2 \times \sin^2(60^\circ)= 3/2$ larger than the corresponding LIGO result due to the LISA having two low-frequency channels, and arms that make an angle of $60^\circ$, as opposed to the $90^\circ$ angle for LIGO. The full expression for ${{\mathcal R}}(\,f)$ has to be computed numerically [5], and has the form shown in figure 2. The transfer function can be well-fit by the curve

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Rfit} {{\mathcal R}}(\,f) = \frac{3}{10} \frac{1} {\left(1 + 0.6 (\,f/f_*)^2 \right)}. \nonumber \end{align} \tag{ 9 }$$

Zoom In Zoom Out Reset image size

Note that many publications quote the number 3/20 for the low frequency limit of $R(\,f)$. The factor of two larger value quoted in equation (9) comes from summing over the two independent low-frequency data channels.

The current mission performance requirement (MPR) model for the power spectral density of the LISA noise P_n(f) is taken from equation (3) of the LISA Science Requirements Document V1 [13]. The single-link optical metrology noise is quoted as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{oms} P_{\rm OMS} = (1.5 \times 10^{-11}~{\rm m})^2~{\rm Hz}^{-1}, \nonumber \end{align} \tag{ 10 }$$

and the single test mass acceleration noise is quoted as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{acc} P_{\rm acc} = (3 \times 10^{-15}~{\rm m}~{\rm s}^{-2})^2 \left(1+ \left(\frac{0.4~{\rm mHz}}{f}\right)^2 \right)~{\rm Hz}^{-1}. \nonumber \end{align} \tag{ 11 }$$

The total noise in a Michelson-style LISA data channel is then [11]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Pn} P_n(\,f) = \frac{P_{\rm OMS}}{L^2} + 2(1+\cos^2(\,f/f_*)) \frac{P_{\rm acc}}{(2\pi f)^4 L^2} . \nonumber \end{align} \tag{ 12 }$$

Note that the Michelson-style response has four contributions from the optical metrology noise and 16 from the test mass acceleration noise. We convert from displacement to strain by dividing by the round-trip light travel distance 2L, so, for example, the factor of $4 P_{\rm OMS}$ gets divided by (2L)², leading to the expression seen in (12). The same factor of $1/(2L)$ is also applied to the path-length change caused by the gravitational wave, so it cancels out in the likelihood function and the SNR, and the choice to divide by 2L is an unimportant convention. A good analytic model for the sensitivity curve that is sufficient for most purposes is given by combining (9) and (12):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Sn2} S_n(\,f) = \frac{10}{3 L^2} \left(P_{\rm OMS} + 2(1+\cos^2(\,f/f_*)) \frac{P_{\rm acc}}{(2\pi f)^4} \right) \left(1 + \frac{6}{10}\left(\frac{f}{f_*}\right)^2 \right). \nonumber \end{align} \tag{ 13 }$$

The $\cos^2(\,f/f_*)$ term multiplying the acceleration noise approaches unity in the frequency range where the acceleration noise dominates, allowing us to write the instrumental contribution to the sensitivity curve in the form shown in equation (1). Figure 3 compares the amplitude spectral density of the instrument noise, $\sqrt{P_n(\,f)}$, to the numerically computed sensitivity curve $\sqrt{S_n(\,f)}$, and the analytic approximation to the sensitivity curve from equation (1). The analytic expression provides a very good fit to the numerically derived sensitivity curve.

Zoom In Zoom Out Reset image size

In addition to the instrument noise, unresolved galactic binaries will act as an effective noise source (though one that is not stationary). The galactic confusion noise goes down as the mission progresses and more foreground sources are removed. Estimates for the confusion noise using the new LISA design are given in [14], and are well fit by the function

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{gf} S_c(\,f) = A \, f^{- 7/3}\, {\rm e}^{-f^{\alpha} + \beta f \sin(\kappa f)} \left[1+{\rm tanh}(\gamma(\,f_k-f))\right]~{\rm Hz}^{-1} \nonumber \end{align} \tag{ 14 }$$

with fit parameters given in table 1. Note that the amplitude quoted here is half the value quoted in [14] since here we are using two-channel sensitivity curves. The full sensitivity curve is found by adding S_c(f) to S_n(f).

Table 1. Parameters of the analytic fit the Galactic confusion noise as described by equation (14). The amplitude A has been fixed to $9\times 10^{-45}$. Note that the knee frequency f_k decreases with observation time and $\gamma$ increase with observation time, leading to a steeper drop off in confusion noise.

	6 mo	1 yr	2 yr	4 yr
$\alpha$	0.133	0.171	0.165	0.138
$\beta$	243	292	299	−221
$\kappa$	482	1020	611	521
$\gamma$	917	1680	1340	1680
fk	0.002 58	0.002 15	0.001 73	0.001 13

Figure 4 shows the contribution of the galactic confusion noise assuming a four year mission, along with the updated sensitivity curve that includes the confusion noise. Note that the confusion noise shown here is an average value—in practice it will vary over a year as the LISA antenna pattern sweeps across the galaxy.

Zoom In Zoom Out Reset image size

The majority of LISA sources will be binaries of various masses and mass ratios. For simplicity we will focus here on quasi-circular, non-spinning comparable mass binaries, and only consider the dominant quadrupole harmonic. Extreme mass ratio binaries, which may be highly eccentric, require a more involved treatment. While we ignore spin, the signal-to-noise ratios we compute should be good to within a factor of two or so for spinning systems, with the main difference coming from the location of the inner-most stable circular orbit [15, 16] and orbital 'hang-up' [17] (note that spin precession causes the inclination and polarization angles to vary and yields SNRs for generic systems that are closer to the angle averaged expression). Our model for the the dominant harmonic of a non-precessing quasi-circular binary is then

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle & \tilde{h}_+(\,f)= A(\,f) \frac{(1+\cos^2\iota)}{2} \, {\rm e}^{{\rm i}\Psi(\,f)}\nonumber \\ & \tilde{h}_\times(\,f) = {\rm i}A(\,f) \cos\iota \, {\rm e}^{{\rm i}\Psi(\,f)} \, , \nonumber \end{align} \tag{ 15 }$$

where $\iota$ describes the inclination of the orbit relative to the line of sight, and $A(\,f)$ and $\Psi(\,f)$ are the amplitude and phase of the wave. To compute the sky/polarization averaged SNR we only need to know $A(\,f)$, and in some cases, how the frequency evolves with time, $f(t)$.

Earlier we related the sky and polarization averaged power spectral density of the signal to the power spectral density seen in the detector via equation (4). For binary systems it is natural to extend the angle averaging to include the inclination angle:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \langle \tilde{h}(\,f) \tilde{h}^*(\,f)\rangle = {{\mathcal R}}(\,f) A^2(\,f) \frac{1}{2} \int_{-1}^{1} \left(\frac{\left(1+x^2\right)^2}{4} + x^2\right) {\rm d}x = \frac{4}{5} {{\mathcal R}}(\,f) A^2(\,f). \nonumber \end{align} \tag{ 16 }$$

Note that for LIGO we recover the well-known pre-factor $\sqrt{(4/5) {{\mathcal R}}(\,f)}= 2/5$ that is applied to the GW amplitude to account for averaging over the source location and orientation [18].

The amplitude signal-to-noise ratio $\rho$ for a deterministic signal $\tilde{h}(\,f)$ is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho^2 = 4 \int \frac{\vert \tilde{h}(\,f) \vert^2}{P_n(\,f)} \, {\rm d}f = 4 \int_{f=0}^\infty \frac{f \vert \tilde{h}(\,f) \vert^2}{P_n(\,f)} \, d(\ln f) \, . \nonumber \end{align} \tag{ 17 }$$

Averaging over sky location, inclination and polarization we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \overline{\rho^2} = \frac{16}{5} \int \frac{f {{\mathcal R}}(\,f) A^2(\,f)} {P_n(\,f)} \,d(\ln f) = \frac{16}{5} \int \frac{(2 f T) S_h(\,f)}{S_n(\,f)} \, d(\ln f). \nonumber \end{align} \tag{ 18 }$$

Where T is the observation time and S_h(f) is the one-sided, angle averaged, power spectral density of the signal,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S_h(\,f) = \frac{A^{2}(\,f)}{2 T} \, . \nonumber \end{align} \tag{ 19 }$$

If you took a Fourier transform of the data, d  =  h  +  n, then ignoring any correlations between the signal and the noise, the power spectral density of the data would equal to $S_d(\,f) = S_h(\,f) + P_n(\,f)$. In other words, S_h(f) is the power spectral density of the signal. The factor of $(2 f T)$ that appears in the expression for the optimal signal-to-noise shows that the signal is effectively boosted relative to the noise by using templates to coherently extract the signal. Rather than plotting the signal power directly (which often lies below the sensitivity curve), the convention is to plot $h_{\rm eff}^2=16 f (2 f T) S_h(\,f)/5$, to account for the boost we get from the coherent signal extraction.

For the waveform model we use the original phenomenological inspiral-merger-ringdown (IMR) model, known as PhenomA [19]. While more accurate models now exist, such as the latest PhenomP model [20, 21], which includes spin-precession, PhenomA is sufficient for making graphs and estimating SNRs. The PhenomA amplitude is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{PhA} A(\,f) &\equiv \sqrt{\frac{5}{24}}\frac{(G{{\mathcal M}}/c^3)^{5/6}f_0^{-7/6}}{\pi^{2/3}(D_L/c)}\, \left\{\begin{array}{@{}ll@{}} \left(\frac{f}{f_0}\right)^{-7/6} & {\rm if}~f < f_0 \nonumber \\ \nonumber \\ \left(\frac{f}{f_0}\right)^{-2/3} & {\rm if}~f_0 \leqslant f < f_1 \nonumber \\ \nonumber \\ w \, {{\mathcal L}}\left(\,f,\,f_1,\,f_2\right) &{\rm if}~f_1 \leqslant f < f_3, \end{array} \right.\nonumber \\ \nonumber \end{align} \tag{ 20 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f_k \equiv \frac{a_k \eta^2 + b_k \eta + c_k}{\pi (G M/c^3)}\, , \nonumber \end{align} \tag{ 21 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\mathcal L}}(\,f,f_1,f_2) \equiv \left(\frac{1}{2 \pi}\right) \frac{f_2}{(\,f-f_1)^2+f_2^2/4} , \label{eq:Lorenzian} \nonumber \end{align} \tag{ 22 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle w \equiv \frac{\pi f_2}{2} \left(\frac{f_0}{f_1}\right)^{2/3}. \nonumber \end{align} \tag{ 23 }$$

Here $M=m_1+m_2$ is the total mass, $\newcommand{\e}{{\rm e}} \eta = m_1 m_2/M^2$ is the symmetric mass ratio and ${{\mathcal M}}= (m_1 m_2){}^{3/5}/M^{1/5}$ is the chirp mass. The G's and c's have been included for those that are not used to working in natural units. Note that the combinations GM/c³ and D_L/c both have units of time. A useful number to remember is that the mass of the Sun, $G M_\odot /c^3$, is approximately 5 $\mu$s in natural units. The coefficients for the transition frequencies f_k are given in table 2. Roughly speaking, f₀ is the merger frequency, f₁ is the ringdown frequency, f₂ is decay-width of the ringdown and f₃ is the cut-off frequency.

Table 2. Polynomial coefficients of the transition frequencies.

	ak	bk	ck
f0	2.9740 $\times 10^{-1}$	4.4810 $\times 10^{-2}$	9.5560 $\times 10^{-2}$
f1	5.9411 $\times 10^{-1}$	8.9794 $\times 10^{-2}$	1.9111 $\times 10^{-1}$
f2	5.0801 $\times 10^{-1}$	7.7515 $\times 10^{-2}$	2.2369 $\times 10^{-2}$
f3	8.4845 $\times 10^{-1}$	1.2848 $\times 10^{-1}$	2.7299 $\times 10^{-1}$

The final ingredient we need for computing the SNR is the frequency range covered by the signal. For comparable mass black holes, with $M > 10^4 M_{\odot}$, the signal will sweep across the LISA band and merge in less than the mission lifetime. However, for lower mass systems, such as stellar origin black holes that will merger in the LIGO band a decade or so later, or for white dwarf binaries, which may be millions of years from merger, we need to specify the start and end frequencies for the SNR integration. To leading post-Newtonian order, the frequency as a function of time is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{ft} f(t) = \frac{1}{8 \pi (G{{\mathcal M}}/c^3)} \left(\frac{5 (G {{\mathcal M}}/c^3)}{t-t_c} \right)^{3/8}, \nonumber \end{align} \tag{ 24 }$$

where t_c is the time of coalescence. For example, an equal mass binary at z  =  3 with a total source-frame mass of $M = 10^6 \, M_{\odot}$ will have a GW frequency of $2.93\times10^{-5}$ Hz one year prior to merger. Note that it is the detector frame mass, M_z  =  M(1  +  z), that should be used in equations (20) and (24). For these high mass systems it makes sense to plot $h_{\rm eff}$ across the entire LISA band, and not worry about setting limits in the SNR integration. Tracks in $\sqrt{S_h}$ and $h_{\rm eff}$ for the aforementioned source are shown in figure 5. In contrast, as shown in figure 6, a source similar to GW150914 that is five years from merger when LISA turns on will sweep from f  =  16 mHz to f  =  29 mHz over the nominal four year mission lifetime.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

For galactic binaries the time to merger is typically very large compared to the mission lifetime, and the frequencies will evolve very little over the course of the mission. Taylor expanding (24) we find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle f(t) = f_{\rm in} + \frac{96}{5} \pi^{8/3} (G {{\mathcal M}}/c^3)^{5/3} f_{\rm in}^{11/3} (t-t_{\rm in}) + \dots \, , \nonumber \end{align} \tag{ 25 }$$

where $f_{\rm in}$ is the GW frequency at the start of the observation, at time $t_{\rm in}$. For typical galactic binaries the change in frequency $\Delta f$ during the mission lifetime is so small that it no longer makes sense to the plot the signals as tracks. Rather, the signals are plotted as points with an amplitude $h_{\rm GB}$ that follows from evaluating the SNR integral:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \overline{\rho^2} = \frac{16}{5} \int_{f_{\rm in}}^{\,f_{\rm in}+\Delta f} \frac{A^2(\,f)} {S_n(\,f)} \,{\rm d}f \approx \frac{16}{5} \frac{\Delta f A^2(\,f_{\rm in})} {S_n(\,f_{\rm in})} \equiv \frac{h^2_{\rm GB}(\,f_{\rm in})}{S_n(\,f_{\rm in})} \nonumber \end{align} \tag{ 26 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{hgb} h_{\rm GB} = \frac{8\, T^{1/2} (G {{\mathcal M}}/c^3)^{5/3}\pi^{2/3} f_{\rm in}^{2/3}}{5^{1/2} (D_L/c)} \, . \nonumber \end{align} \tag{ 27 }$$

For example, as shown in figure 6, SDSS J0651  +  2844 which has $D_L \sim 1$ Kpc, $m_1 \sim 0.5\, M_\odot$, $m_2 \sim 0.25\, M_\odot$, and $f_{\rm in} = 2.6$ mHz [22], will produce a strain spectral density of $h_{\rm GB} = 2.8 \times 10^{-18}~{\rm Hz}^{-1/2}$ and have an SNR of 140 assuming a four year mission lifetime. Of course, it not strictly correct to compute angle averaged SNRs for a source with a known sky location and orientation, nor does it make much sense to plot its amplitude against against an all-sky averaged sensitivity curve, but doing so allows us to put all LISA sources on a single graph.

The differing conventions between how slowly evolving and rapidly evolving signals are plotted can be problematic for stellar origin black hole binaries (SOBHBs). For example, if a GW150914 type system was 20 years from merger when LISA started observations, it would be emitting at a gravitational wave frequency of 9.5 mHz, and four years later it would be emitting at 10.4 mHz, producing a track that runs for just $\Delta \ln f = 0.09$. Since the frequency range is so short, the questions becomes do we treat the system as evolving, and plot a track as we do for massive black holes, or do we treat the system as non-evolving, and plot a point as we do for galactic binaries? Figure 7 shows that the two choices paint an inconsistent picture. If the track is longer, such as for a system that is five years from merger, the two representations look more consistent. To arrive at consistent representations, where sources appear at almost the same height when shown as evolving tracks or non-evolving points, we recommend switching from tracks to points when $\Delta \ln f < 0.5$.

Zoom In Zoom Out Reset image size

2.1. EMRIs and other complicated signals

Some sources produce signals that can not be accurately described by simple frequency domain models. Extreme mass ratio inspirals and rapidly precessing spinning black hole binaries fall into that category. The inclination and polarization averaging we used for quasi-circular binaries is not applicable these systems, but for simplicity we still plot these signals against the standard sky and polarization averaged sensitivity curve. Writing the sky and polarization averaged signal-to-noise squared as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \overline{\rho}_{\Omega,\psi}^2 = 4 \int \frac{\vert \tilde{h}(\,f) \vert^2}{S_n(\,f)} \, {\rm d}f = \int_{f=0}^\infty \frac{4 f^2 \vert \tilde{h}(\,f) \vert^2}{(\,f S_n(\,f))} \, d(\ln f) \, , \nonumber \end{align} \tag{ 28 }$$

indicates that an appropriate quantity to plot against the characteristic sensitivity $(\,f S_n(\,f)){}^{1/2}$ is the dimensionless characteristic strain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle h_c(\,f) = 2 f \left(\vert \tilde{h}_+(\,f) \vert^2+ \vert \tilde{h}_\times(\,f) \vert^2 \right)^{1/2}. \nonumber \end{align} \tag{ 29 }$$

We generate the Barycenter signals h₊ (t) and $h_\times(t)$, Fourier transform, and form h_c(f). To beautify the plots we smooth the numerically generated h_c(f) using a running average over  ∼100 frequency bins. As an example, we generated augmented analytic kludge (AAK) EMRI waveforms [23] using the code provided at GitHub [24] for a $20\, M_\odot$ stellar remnant black hole falling into a $10^6 M_\odot, \chi = 0.5$ spinning supermassive black hole at distance of 4 Gpc, starting four years before merger with an eccentricity of 0.5. This system has ${\rm SNR} = 52$, and the effective strain shown in figure 8.

Zoom In Zoom Out Reset image size

2.2. Sky dependent estimates

The sky averaged signal-to-noise ratios are useful for a first look at what systems might be detectable, but the signal-to-noise ratio can vary significantly across the sky, especially for short duration signals [25]. To incorporate the sky location dependence we must revert to using signals which have not been averaged over the sky location. The sky-location dependent SNRs are particularly useful for sources with known locations, such as the galactic verification binaries. We will continue to average over inclination angle and polarization angle as these are usually not well constrained.

The Michelson-type signal with spacecraft 1 at the vertex is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \quad s_{1}(t) = \frac{\delta \ell_{12}(t-2 L/c) + \delta \ell_{21}(t-L/c)}{2 L} - \frac{\delta \ell_{13}(t-2 L/c) + \delta \ell_{31}(t-L/c)}{2 L}, \nonumber \end{align} \tag{ 30 }$$

where the GW induced variation in LISA arm lengths between spacecraft i and j are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\delta \ell_{ij}(t)}{L} = \frac{1}{2}(1+\cos^{2}\iota)d^{+}_{ij}(t) h^{+}(\xi_{i}) + \cos \iota \, d^{\times}_{ij}(t)h^{\times}(\xi_{i}). \nonumber \end{align} \tag{ 31 }$$

The detector terms $d^{+,\times}_{ij} = d^{+,\times}_{ij}\left(t; {\bf x}_{i}, {\bf x}_{j}, \theta, \phi, \iota, \psi\right)$ [26] describes LISA's geometry through their dependence on the spacecraft position ${\bf x}_{i}$. The variable $\xi_{i}$ define surfaces of constant gravitational wave phase at spacecraft i. The frequency domain representation of this signal can be found using the stationary phase approximation [27]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\widetilde{\delta \ell}_{ij}(\,f)}{L} = \left[\frac{1}{2}(1+\cos^{2}\iota)d^{+}_{ij}\left(t_{*}\right) + {\rm i}\cos \iota \, d^{\times}_{ij}\left(t_{*}\right) \right] A(\,f) {\rm e}^{{\rm i}\left(\Psi(\,f) + \delta \Psi_{i}(\,f) \right)} \nonumber \end{align} \tag{ 32 }$$

where $A(\,f)$ and $\Psi(\,f)$ are the amplitude and phase. The stationary time is given by the relation $t_{*}(\,f) = \Psi'(\,f)/2\pi$. This is used to map $d^{+}_{ij}(t)$ to $d^{+}_{ij}(\,f)$ etc. The motion of the LISA detector also imparts a phase shift $\delta \Psi_{i}(\,f) = 2\pi f \hat{{\bf k}}\cdot{\bf x}_{i}(t_*(\,f))/c$, where $\hat{{\bf k}}$ defines the line of sight vector to the source. The PhenomA amplitude $A(\,f)$ is given in equation (20). The phase is written in terms of a power series expansion that is motivated the post Newtonian expansion, with additional terms that are found by fitting to numerical relativity simulations:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \quad \Psi(\,f) = 2 \pi f t_0 +\phi_0 + \psi_0 f^{-5/3}+\psi_2 f^{-1} + \psi_3 f^{-2/3} + \psi_4 f^{-1/3}+\psi_6 f^{1/3}. \nonumber \end{align} \tag{ 33 }$$

The expansion coefficients $\psi_i$ depend on the masses and are given in terms of a numerical look-up table [19].

The average over the source inclination and polarization angles, $\left< |\tilde{s}_{1}(\,f)|^{2}\right>_{\iota, \psi}$, where the subscripts to the angle brackets denote which variables are averaged over, can be shown to be equivalent to computing the signal at two fixed values of inclination and polarization:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left< |\tilde{s}_{1}(\,f)|^{2}\right>_{\iota, \psi} = \frac{8}{5}\left(|\tilde{s}_{1}(\,f)|^{2}|_{\iota = \frac{\pi}{2}, \psi = 0}+ |\tilde{s}_{1}(\,f)|^{2}|_{\iota = \frac{\pi}{2}, \psi = \frac{\pi}{4}} \right) \, . \nonumber \end{align} \tag{ 34 }$$

This allows us to compute the orientation averaged signal using just two calls to the waveform generator. The sky location dependent signal-to-noise ratio is then

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{SNRtime} \bar{\rho}_{\iota, \psi}^{2}(\theta,\phi) = 4 \int_{f=0}^{\infty} \frac{f^{2} \left< |\tilde{s}_{1}(\,f,\theta,\phi)|^{2}\right>_{\iota, \psi}}{f P_{n}(\,f)} {\rm d} \log f . \nonumber \end{align} \tag{ 35 }$$

To give a visual impression of the signal strength we plot the characteristic strain $\tilde{h}_{\mathrm{eff}}(\,f) = 2 f \left< |\tilde{s}_{1}(\,f)|\right>_{\iota, \psi}$ against the characteristic noise amplitude in the Michelson channel $\sqrt{f P_{n}(\,f)}$. For a galactic binary we take a similar approach to calculate SNRs as for the sky averaged case:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \bar{\rho}_{\iota, \psi}^{2}(\theta,\phi) = 4 \int_{f_{\mathrm{in}}}^{f_{\mathrm{in}}+\Delta f} \frac{\left< |\tilde{s}_{1}(\,f)|^{2}\right>_{\iota, \psi}}{P_{n}(\,f)} {\rm d} f \approx 4 \frac{\left< |\tilde{s}_{1}(\,f_{\mathrm{in}})|^{2}\right>_{\iota, \psi}}{P_{n}(\,f_{\mathrm{in}})} \Delta f . \nonumber \end{align} \tag{ 36 }$$

An example of a sky-location dependent sensitivity plot is shown in figure 9 for the same sources shown previously in figure 6, but now at two different sky locations ($\theta = 0.5,\phi = 2.3)$ and ($\theta = 1.1,\phi = 1.5)$. The signal-to-noise ratio for these sky-dependent sources are 3106 and 2017 for the super massive black hole binary, 207 and 151 for the galactic binary, and 4.39 and 4.51 for the stellar origin black hole binary.

Zoom In Zoom Out Reset image size

The oscillations in the tracks seen in figure 9 are due to the time dependent antenna pattern. For the supermassive black hole system most of the modulation is seen at low frequencies where the system spends many months. The evolution of the signal becomes far more rapid as it sweeps to higher frequencies, and the detector is effectively stationary on this timescale, so the amplitude no longer oscillates. The oscillations are very pronounced for the stellar origin black hole binary, which slowly evolves over the entire four year observation period.

The sky location of the source plays a very large role in its detectability and our ability to characterize the source parameters. For example, the standard deviation of the SNR across the sky for the equal mass binary black hole system shown in figure 6 is 76% of the sky averaged value. The galactic binary and LIGO binary had standard deviations of 77% and 39% of the sky averaged value respectively. In figure 10 we show the variation in the SNR across the sky for two binary black holes, one with detector frame total mass $4 \times 10^6 M_\odot$, and another with detector frame total mass $1.6 \times 10^4 M_\odot$. For the more massive system, most of the SNR is accumulated in the two days around merger, and since the detector is effectively stationary during this time, we recover the fixed quadrupolar antenna pattern. For the less massive system, where the merger lasts for a longer time, the sky map is more uniform.

Zoom In Zoom Out Reset image size

Figure 11 shows the signal-to-noise grows in time for the stellar origin black hole binary for various sky locations. We see that the rate of SNR growth depends on the sky location at any given time due to the changing orientation of the LISA antenna pattern. If the sources lies in a sensitive region of the antenna pattern then the SNR will grow quickly compared to when the source lies in an insensitive region.

Zoom In Zoom Out Reset image size

An additional factor that will impact the SNR as a function of sky location is the time variation of the galactic confusion noise. The confusion noise is loudest when the most the peak of the antenna pattern sweeps across the galaxy, and quietest when pointed away from the galaxy. The confusion noise will vary adiabatically such that S_c(f,t). Using the stationary phase approximation, the time dependence gets mapped to frequency dependence via t_*(f), which allows for the time variation of the confusion noise to be incorporated in the SNR integral (35). We defer an analysis of this effect to a future study.

We have explained the conventions and factors that enter into the calculation of angle-averaged signal-to-noise ratios for some of the key science targets for the future LISA gravitational wave observatory. Simple Python codes have also been provided to perform these calculations, including for cases where the sky location is known, as will be the case for galactic 'verification' binaries. In future we plan to extend the calculations to more complicated signals, and to include the time-variation in the confusion noise. We also plan to extend the calculations to include simple Fisher matrix based error estimates for the parameter estimation errors. Our hope is that these tools will help astronomers explore the science potential of the LISA mission.

We appreciate the input and feedback from Martin Hewitson, Emanuele Berti, Paul McNamara, and Davide Gerosa. TR and NJC appreciate the support of the NASA grant NNX16AB98G. CL acknowledges the support by the National Natural Science Foundation of China Grants No.11690021 and the UCAS Joint PhD Training Program.