LEGWORK: A python package for computing the evolution and detectability of stellar-origin gravitational-wave sources with space-based detectors

We present LEGWORK (LISA Evolution and Gravitational Wave Orbit Kit), an open-source Python package for making predictions about stellar-origin gravitational wave sources and their detectability in LISA or other space-based gravitational wave detectors. LEGWORK can be used to evolve the orbits of sources due to gravitational wave emission, calculate gravitational wave strains (using post-Newtonian approximations), compute signal-to-noise ratios and visualise the results. It can be applied to a variety of potential sources, including binaries consisting of white dwarfs, neutron stars and black holes. Although we focus on double compact objects, in principle LEGWORK can be used for any system with a user-specified orbital evolution, such as those affected by a third object or gas drag. We optimised the package to make it efficient for use in population studies which can contain tens-of-millions of sources. This paper describes the package and presents several potential use cases. We explain in detail the derivations of the expressions behind the package as well as identify and clarify some discrepancies currently present in the literature. We hope that LEGWORK will enable and accelerate future studies triggered by the rapidly growing interest in gravitational wave sources.


INTRODUCTION
The planned space-based gravitational wave detector LISA (Laser Interferometer Space Antenna) will present an entirely new view of gravitational waves by focusing on lower frequencies (10 −5 < f / Hz < 10 −1 ) than ground-based detectors. This will enable the study of many new source classes including mergers of supermassive black holes (e.g. Begelman et al. 1980;Klein et al. 2016;Bellovary et al. 2019), extreme mass ratio inspirals (e.g. Berti et al. 2006; Barack & Cutler 2007;Babak et al. 2017;Moore et al. 2017), and cosmological GW backgrounds (e.g. Bartolo et al. 2016;Caprini et al. 2016;Caldwell et al. 2019). However, this frequency regime is also of interest for the detection of local stellar-mass binaries during their inspiral phase. LISA is expected to detect Galactic stellar-origin binaries containing combinations of white dwarfs, neutron stars, and black holes, ranging from the numerous double white dwarf population, to the rare but loud double black hole population.
Each of these studies require making estimates of the signal-to-noise ratio of individual binary systems and possibly the slow gravitational wave inspiral that lead to the present-day parameters. So far, most studies made use of custom made codes which have not been made publicly available.
We believe that the large renewed interest in LISA and the stellar-origin sources it may detect will lead to many more studies in the near future that would need similar computations. This leads to a significant amount of redundancy which, at best results in extra work for each individual and at worst leads to an increased chance of introducing mistakes and inconsistencies when translating the necessary expressions to software.
LEGWORK is an open-source Python package designed to streamline the process of making predictions of LISA detection rates for stellar-origin binaries such that it is as fast, reliable and simple as possible. With LEGWORK one can evolve the orbits of a binary or a collection of binaries and calculate their strain amplitudes for any range of frequency harmonics. One can compute the sensitivity curve for LISA or other future gravitational wave detectors (e.g. TianQin's curve, or that of a custom instrument) and use it to compute the signal-to-noise ratio of a collection of sources. Furthermore, LEGWORK provides tools to visualise all of the results with easy-to-use plotting functions. Finally, LEGWORK is fully tested to check for consistency in the derivations described below.
Specifically, we implement the post-Newtonian expressions by Peters & Mathews (1963) and Peters (1964) for the evolution of binary orbits due to the emission of gravitational waves, equations for the strain amplitudes and signal-to-noise ratios of binaries from various papers (e.g. Flanagan & Hughes 1998;Finn & Thorne 2000;Cornish & Larson 2003;Barack & Cutler 2004;Moore et al. 2015) and approximations for the LISA and TianQin sensitivity curves given in Robson et al. (2019) and Huang et al. (2020) respectively. The post-Newtonian expressions are approximately of order 0.5 as they account for the orbital evolution from energy loss due to the emission of GWs, but don't include other effects such as spin-orbit coupling. We find this is an excellent approximation for stellar-origin sources in the LISA band as any higher order terms are either zero or negligible.
The open-source nature of the project means that new users as well as seasoned experts in the field can work together in a collaborative setting to consider new fea-tures and enhancements to the package as well as check the implementation. At the same time, with our thorough online documentation, derivations and tutorials, we hope LEGWORK can make this functionality more accessible to the broader scientific community.
We note that LEGWORK is not the only initiative of this kind. We highlight the "Gravitational Wave Universe Toolbox" presented by Yi et al. (2022) which was developed to simulate observations on the GW universe with different detectors covering the full gravitationalwave spectrum and source classes. We also highlight "GWPlotter" by (Moore et al. 2015), which provides an interactive plotting tool to compare the sensitivity of different gravitational wave detectors. LEGWORK differs from the tools listed here as we have tried to provide tools that are optimised to rapidly make predictions for large populations of stellar-origin sources and we have focused on space-based detectors.
LEGWORK has been developed with stellar-origin binary population studies in mind. We highlight two recent papers that use our package. Wagg et al. (2021) investigates several populations of potential LISA sources (double black holes, black hole neutron stars and double neutron stars). They use LEGWORK to evolve each synthesised source over the age of the Milky Way, make predictions about the LISA detectable population and explore how it varies with different binary physics assumptions. Additionally, Thiele et al. (2021) examines the implications of assuming a metallicity-dependent binary fraction for the formation of close double white dwarfs (WDWDs) in the LISA frequency band. They use LEGWORK to calculate the SNR of WDWDs and, in particular, apply a custom noise model that combines the LISA sensitivity curve from Robson et al. (2019) with a new fit for the Galactic confusion noise based on their WDWD population.
LEGWORK can be installed with pip or obtained from the GitHub repository 1 . All examples shown in this paper and code to reproduce the figures are available in the repository. Instructions for installation and basic usage are provided in the online documentation 2 which contain the most up-to-date instructions.
The paper is organised as follows. In Section 2 we give an overview of the capabilities of the LEGWORK package and its various modules. We detail the derivations of the equations relevant for LEGWORK in Section 3. In Section 4 we outline some example use cases of LEG-1 https://github.com/TeamLEGWORK/LEGWORK 2 https://legwork.readthedocs.io/en/latest/ WORK to demonstrate its use. Finally, in Section 5 we conclude and summarise our work.

PACKAGE OVERVIEW
The LEGWORK package is composed of seven modules that each focus on a particular aspect of calculations useful for gravitational wave sources that are detectable by space-based detectors. In Figure 1, we illustrate the general structure of the package with each of its modules. The source module is the central module of the package and provides an simple interface to the functions in the rest of the modules. For more complex analyses, users may want to interact directly with individual modules, particularly those in the top row of Figure 1 as they comprise the core functionality of LEGWORK. Below we explain the capabilities of each module in detail.
evol handles the orbital evolution of a binary due to the emission of gravitational waves. It includes functions for computing the merger times of both circular and eccentric binaries. In addition, you can use this module to evolve binary orbit parameters forward in time with any number or arrangement of timesteps. We discuss the relevant equations in Section 3.2.
strain contains two functions that compute a binary's gravitational wave strain and characteristic strain amplitude respectively. Each of these functions is capable of computing the strain for an array of binaries at any number of timesteps and evaluated at any number of frequency harmonics. We discuss the relevant equations in Section 3.3.
psd is used for evaluating the effective noise power spectral density of a detector at different frequencies.
The module currently contains the LISA and TianQin sensitivity curves that can be tweaked by adjusting parameters such as the observation time, response function and even the arm length. For the Galactic confusion noise, a user can choose one of three models (Robson et al. 2019;Huang et al. 2020;Thiele et al. 2021), use their own custom model, or turn it off entirely. Additionally, this module allows the user to specify a custom detector sensitivity curve. We discuss the relevant equations in Section 3.4.
snr uses the functions in evol, strain and psd to compute the signal-to-noise ratio of sources. It contains four functions that cover the permutations of whether a source is circular or eccentric and stationary in frequency space on the timescale of the mission or evolving. We discuss the relevant equations in Section 3.5.
visualisation contains several wrappers for plotting 1-and 2-dimensional distributions with histograms, scatter plots, and kernel density estimator (KDE) plots in order to quickly analyse a collection of sources. In addition, it provides functions for plotting sources directly onto a sensitivity curve.
source provides a direct and simple interface to the functions in other modules through the Source Class. You can instantiate this Class with an array of sources and use it compute their strains or signal-to-noise ratios directly. Moreover, depending on the user's choice of allowed gravitational wave luminosity error, the Class dynamically decides on the number of frequency harmonics needed to capture the full signal of each binary and at what eccentricity to no longer consider a binary circular. This Class also provides a quick means of evolving the sources, visualising the parameters of each source and allows you to plot the binaries on the sensitivity curve.
utils is a collection of miscellaneous utility functions mainly consisting of conversions between variables as well as constants and expressions from Peters (1964). We discuss the relevant equations in Section 3.1.

Units and automated testing
To ensure stability with physical units, all quantities included in LEGWORK use the astropy.units module (Astropy Collaboration et al. 2013. This means that all inputs to LEGWORK can be given in the units of the user's choice and will be automatically converted.
Furthermore, all of the source code in LEGWORK is fully tested with continuous integration in the LEG-WORK GitHub repository. We employ several unit tests to ensure consistency between each of the use cases described below. For example, we require that the SNR calculation for circular and stationary binaries produces consistent output whether LEGWORK uses the stationary and circular approximations or not. Similarly, we verify that the antenna patterns described below produce the expected values when averaged over source positions, inclinations, and polarisations.

Optimistations
We developed LEGWORK with an emphasis on increasing the efficiency of these computations in order to make simulations of large populations of systems tractable. We ensured that the entirety of LEGWORK is vectorised and thus scales well with larger populations. In addition, we made a several specific optimisations to further increase the speed of calculations.
Firstly, we find that the runtime of calculating strains and SNRs for large populations of sources is mainly limited by the computation of (1) the relative gravitationalwave power in each harmonic for eccentric systems (see Eq. 5) (2) the sensitivity curve of the given detector. Therefore, in order to significantly reduce the runtime of strain and SNR computations, by default LEGWORK automatically interpolates these functions upon instantiation of any Source class with a large number of sources. Thereafter all functions use the tabulated values instead of calculating them exactly.
In certain cases one can apply approximations in place of the general SNR calculations. Although it is possible to use each of these approximations directly through the snr module, the source module will automatically apply the most appropriate function for each individual source. LEGWORK dynamically classifies each source as one of four types, which cover the permutations of whether a binary is effectively circular or eccentric and whether or not it is stationary in frequency space on the timescale of the LISA mission. Then when computing the SNR, it applies a different function to each type of the source. This avoids computing unnecessary, time intensive integrals.
Moreover, for increasingly eccentric sources, gravitational waves are emitted in an increasing number of higher frequency harmonics. Although the total SNR is formally calculated as a the sum of the SNR over an infinite number of harmonics, in practice it is sufficient to only consider a subset. In the source module, the user can provide gw lum tol, a maximum allowable tolerance for the accuracy of the gravitational wave luminosity. Given this tolerance, LEGWORK automatically calculates the required number of harmonics to satisfy this tolerance and thus minimise the computation time.

DERIVATIONS
In this section, we present a derivation of the equations used in LEGWORK. We emphasise that these are not new derivations, conversely, they are in fact given frequently in the literature. However, they are often incomplete, unclear and, in some cases, contain spurious constant factors that arise from invalid combinations of previous work. Here, we aim to present a clear, clean and concise explanation of the expressions we use in LEGWORK.
A symbol in an equation directly links to the relevant online LEGWORK documentation for the implementation of that equation, which additionally contains a link to exact code used to reproduce the equation. These derivations are also given in more detail in the LEGWORK documentation, where we show each of the intervening steps in more detail.

Conversions and definitions (utils)
We start these derivations by defining some useful conversions and definitions. The chirp mass of a binary is the mass quantity measured by LISA and is given by where m 1 and m 2 are the primary and secondary masses of the binary. It is often convenient to convert between orbital frequency, f orb , and the semi-major axis, a, of a binary and this can be accomplished with Kepler's third law.
where G is the gravitational constant. Inversely, For circular binaries, gravitational-wave emission occurs at twice the orbital frequency (f GW = 2f orb ). However, for eccentric binaries, we need to consider all frequency harmonics of gravitational-wave emission. These are defined such that the n th harmonic frequency is It will be important to know the relative gravitationalwave power radiated into the n th harmonic for a binary with eccentricity e for the strain and SNR calculations. This is given by Peters & Mathews (1963, Eq. 20) 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 equivalent circular binary. This enhancement factor is (Peters & Mathews 1963, Eq. 17) F (e) ≡ ∞ n=1 g(n, e) = 1 + 73 24 e 2 + 37 96 e 4 (1 − e 2 ) 7/2 .
A useful rule of thumb is that F (0.5) ≈ 5.0, or in words, a binary with eccentricity 0.5 loses energy to gravitational waves at 5 times the rate of an equivalent circular binary.

Circular binaries
For a circular binary, the orbital evolution due to gravitational-wave emission can be calculated analytically, as the rate at which the separation of the binary shrinks is simply a function of its mass and the current separation (Peters 1964, Eq. 5 where the constant β is defined as where c is the speed of light, m 1 is the primary mass and m 2 is the secondary mass. This gives the semimajor axis of a circular binary as a function of time, t, as (Peters 1964, Eq. 5.9) where a 0 is the initial semi-major axis. Moreover, we can solve for the merger time, the time until the binary will merge, by setting the final semi-major axis in Eq. 9 to zero.

Eccentric binaries
The orbital evolution is more complex for eccentric binaries since the semi-major axis and eccentricity both evolve simultaneously and depend on one another. The final expression cannot be solved analytically and require numerical integration. The semi major axis, a, and eccentricity, e, are related as (Peters 1964, Eq. 5.11) a(e) = c 0 e 12/19 (1 − e 2 ) 1 + 121 304 where c 0 satisfies the initial conditions such that a(e 0 ) = a 0 . The time derivative of the eccentricity, e, is which we can integrate to find e(t) and convert to a(t) using Eq. 11 ( ). Inverting this function and applying the fact that we know that e → 0 when the binary merges gives the merger time (Peters 1964, Eq. 5.14) For very small or very large eccentricities we approximate this integral using the following expressions (given in unlabelled equations after Peters 1964, Eq. 5.14) t merge, (1−e 2 ) 1 = 768 425 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 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 which gives t 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 LEG-WORK.

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 th harmonic of the binary will spend approximately f n /ḟ n seconds (or f 2 n /ḟ n cycles) in the vicinity of a frequency f n (Finn & Thorne 2000). 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 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 th harmonic, h 0,n , as (e.g Finn & Thorne 2000; Moore et al. 2015 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 th harmonic is given by (e.g. Barack & Cutler 2004, Eq. 56;Flanagan & Hughes 1998, Eq. 5 where D L is the luminosity distance to the source (note that for Milky Way sources, or any sources with redshift ∼ 0, this is simply the distance to the source),Ė n is the power radiated in the n th harmonic andḟ n is the rate of change of the n th harmonic frequency. The power radiated in the n th harmonic can be expressed as (Peters & Mathews 1963, Eq. 19) where g(n, e) is given in Eq. 5. By substituting a for f orb (using Eq. 3) and applying the definition of the chirp mass (Eq. 1) we obtain a more useful form for making gravitational wave predictionṡ The last term needed to define the characteristic strain in Eq. 18 is the rate of change of the n th harmonic frequency as a result of gravitational wave inspiral, which we can write asḟ 3 Note that this is factor of 2 different from Finn & Thorne (2000). This is because the factor of 2 is already included in the Robson et al. (2019) sensitivity curve and so is removed here.
We can find an expression for df n / da by substituting Eq. 3 into Eq. 4 and differentiating The rate at which the semi-major axis decreases is (Peters 1964, Eq. 5.6) Substituting Eq. 22 and Eq. 23 into Eq. 21 gives an expression forḟ ṅ f n = 48n 5π which, as above withĖ n , we can recast using Kepler's third law and the definition of the chirp masṡ With definitions of bothĖ n andḟ n , we are now in a position to find an expression for the characteristic strain by plugging Eq. 20 and Eq. 25 into Eq. 18:

Strain
In order to obtain an expression for the strain amplitude of gravitational waves in the n th harmonic, we can use Eq. 17 and plug in Eq. 25 and Eq. 26

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 unknown positions, inclinations, and polarisations, we use an average for the detector. However, for sources where these quantities are known we can consider the amplitude modulation of the strain due to the average motion of LISA's orbit.
We write that the position of the source on the sky is given by the ecliptic coordinates (θ, φ), the inclination of a source is ι and the polarisation of a source (determined by its orientation relative to the detector) is given by ψ. We follow the results of Cornish & Larson (2003) to define the amplitude modulation. However, we adapt their expression to remain in the frequency domain and follow the conventions of more recent papers (e.g. Babak et al. 2021, Eq. 67) to write the amplitude modulation as where F 2 + orb and F 2 × orb , the orbit-averaged detector responses, are defined as and D + D × orb = 243 512 cos θ sin 2φ 2 cos 2 φ − 1 D 2 × orb = 3 512 120 sin 2 θ + cos 2 θ + 162 sin 2 2φ cos 2 θ , D 2 + orb = 3 2048 487 + 158 cos 2 θ + 7 cos 4 θ − 162 sin 2 2φ 1 + cos 2 θ 2 .
The orbital motion of LISA smears the source frequency by roughly 10 −4 mHz due to the antenna pattern changing as the detector orbits, the Doppler shift from the motion, and the phase modulation from the + and × 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 the amplitude modulation is only implemented in LEGWORK for quasi-circular binaries to remain consistent with the calculation in Cornish & Larson (2003). Since the expected use case of LEGWORK is estimation of the detectability of large populations of mHz stellar-remnant binaries, for which predictions are uncertain by orders of magnitude in some cases, an extension of Cornish & Larson (2003) which includes eccentric binaries is out of the current scope of LEG-WORK.

LISA
For the LISA sensitivity curve we use the equations from Robson et al. (2019). The effective strain spectral density of the noise is defined as where P n (f ) is the power spectral density of the detector noise and R(f ) is the sky and polarisation averaged signal response function of the instrument. Alternatively if we expand out P n (f ), approximate R(f ) and simplify we find (Robson et al. 2019, Eq. 1) where L = 2.5 Gm is detector arm length, f * = c/2πL = 19.09 mHz is the transfer frequency, is the single-link optical metrology noise (Robson et al. 2019, Eq. 10), is the single test mass acceleration noise (Robson et al. 2019, Eq. 11) and is the galactic confusion noise (Robson et al. 2019, Eq. 14), where the amplitude A is fixed as 9 × 10 −45 and the various parameters change over time and are listed in Robson et al.

TianQin
For the TianQin sensitivity curve we use the power spectral density given in Huang et al. (2020, Eq. 13) where L = √ 3 × 10 5 km is the arm length, S a = 1 × 10 −30 m 2 s −4 Hz −1 is the acceleration noise, S x = 1 × 10 −24 m 2 Hz −1 is the displacement measurement noise and f * = c/2πL is the transfer frequency. Note that Eq. 39 includes an extra factor of 10/3 compared to Huang et al. (2020, Eq. 13). Huang et al. (2020) absorb this factor into the waveform rather than include it in the power spectral density. We include it to match the convention used by Robson et al. (2019) for the LISA sensitivity curve (see the factor of 10/3 in Eq. 35) so that the sensitivity curves can be compared fairly.

SNR for LISA (snr)
We note that this section draws heavily from Flanagan & Hughes (1998) Section II C.

Defining general SNR
In order to calculate the signal to noise ratio for a given source of gravitational waves (GWs) in a 6-link LISA detector, we need to consider the following parameters: • position of the source on the sky: (θ, φ) • direction from the source to the detector: (ι, β) • orientation of the source, which fixes the polarisation of the GW: ψ • 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, ρ, is given in the frequency domain as wheres(f ) is the Fourier transform of the signal, s(t), and P n (f ) is the one-sided power spectral density of the noise defined as as n(t) n(t) = ∞ 0 1 2 P n (f )df (c.f. Robson et al. 2019, Eq. 2). Here,s(f ) is implicitly also dependent on D L , θ, φ, ψ, ι, and β as where F +,× are the 'plus' and 'cross' antenna patterns of the LISA detector to the 'plus' and 'cross' strains, h +,× . Note throughout any parameters discussed with the subscript x +,× refers to both x + and x × .

Average over position and polarisation
Now, we can consider averaging over different quantities. In LISA's case, when averaged over all angles and polarisations, the antenna patterns are orthogonal and thus F + F × = 0. This means we can rewrite Eq. 43 as which can then be applied to Eq. 40 to give Then combining Eq. 45 and Eq. 46, we then find 4

Average over orientation
Now, we can average over the orientation of the source: (ι, β), noting that the averaging is independent of the distance D L . With this in mind, we can rewrite |h + | 2 + |h × | 2 in terms of two functions |H + | 2 and |H × | 2 , where h +,× =H +,× /D L . Given this, averaging over the source direction gives where we would like to expressH +,× (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. Press & Thorne 1972, Eq. 6) where the bar indicates an average over several cycles of the wave which is appropriate for LISA sources. We can transform Eq. 49 using Parseval's theorem, where we can write Note that the factor of frequency squared comes from the Fourier transform of the square of the time derivative in Eq. 50. Now, since A = D 2 L Ω and |h +,× | 2 = |H +,× | 2 /D 2 L , we find |h +,× | 2 dA = |H +,× | 2 dΩ ι,β , then we can write Eq. 50 in terms of |H +,× | 2 as Alternatively, by using Eq. 51 and performing a Fourier transform we can also write that From inspection of Eq. 52 and Eq. 53, we can write the spectral energy flux as

Fulled averaged SNR
We are now in a position to write an expression for the fully averaged SNR. Note for brevity we write ρ 2 when referring to ρ 2 (θ,φ,ψ),(ι,β) . The application of Eq. 54 to Eq. 48 yields This simplifies nicely to Finally, noting that dE/df = dE/dt × dt/df =Ė/ḟ , we can use the definition of the characteristic strain from Eq. 18 to finish up our position, direction, and orientation/polarisation averaged SNR as where we have used that the effective power spectral density of the noise is defined as S n (f ) = P n (f )/R(f ). 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. 57 to eccentric binaries with where h c,n is defined in Eq. 26 and S n in Eq. 35.

SNR Approximations
Although Eq. 58 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 ≈ f n,f and so the small interval allows one to approximate the integral. Note we refer to non-stationary binaries as 'evolving' here though many papers also use 'chirping'.
For an evolving and eccentric binary, no approximation can be made and the SNR is found using Eq. 58.
For a evolving and circular binary, the sum can be removed and so the SNR found as For a stationary and eccentric binary we can approximate the integral.
where we have applied Eq. 17 to convert between strains. This gives the following expression Finally, for a stationary and circular binary the signalto-noise ratio is For SNR calculations that take into account the amplitude modulation due to LISA's orbital motion (Section 3.3.3, we apply the calculations as described above but include the modulation in either the strain or characteristic strain as is appropriate.

USE CASES
In this section, we demonstrate LEGWORK's range of capabilities through a series of example use cases. The plots and results in each subsection are reproduced directly in online demos in the LEGWORK documentation, which are each based on individual Jupyter notebooks. These tutorials are linked at the start of each subsection with the icon.

Computing the SNR of a binary system -
The most fundamental use case of LEGWORK is to compute the SNR of an individual binary system. This can be accomplished in LEGWORK with only two lines of code -one line to set up the source and another to compute the SNR. As an example one could consider a binary with the parameters m 1 = m 2 = 10 M , d = 8 kpc, f orb = 10 −4 Hz, e = 0.2, where m 1 , m 2 are the primary and secondary masses, d is the distance to the source, f orb is the orbital frequency and e is the eccentricity. As we do not specify a position, polarisation or inclination, LEGWORK will calculate the SNR averaged over these quantities. If a user specifies the position, then if either of the polarisation or inclination are not specified then they will be randomly generated. LEGWORK does not allow users to specify polarisations or inclinations without positions as we presume that the user must know the position if they know the polarisation or inclination of the source. One can now instantiate a source in LEGWORK using these parameters (leaving the mission parameters as the default values, such that we compute the SNR for a 4-year LISA mission). As shown in the linked demo, LEGWORK quickly computes that the SNR of this binary is 4.49. In the background, LEGWORK decides which SNR approximation is most applicable given the eccentricity of the binary and whether it is stationary in frequency space.
This can be generalised to a population of many binary systems with ease. Instead of inputting single values for each parameter, one can input arrays of values where each entry corresponds to a different binary.
As an example, we can take the same parameters as above for 3 different binaries but vary the primary as m 1 = [5, 10, 15] M . Using LEGWORK we find that the SNR for each of these cases is ρ = [2. 47, 4.49, 7.85]. This also need not be limited to a 4-year LISA mission. With LEGWORK we can additionally specify various parameters for the detector. For example, using LEG-WORK we can find that, for a 5-year TianQin mission with no confusion noise, the SNR for each of these cases is ρ = [1.07, 1.95, 3.41].

Horizon distance -
A common question to consider with stellar mass sources in LISA is how far away a certain source could be detected. In other words, what is the horizon distance beyond which a source no longer has a SNR greater than some chosen threshold. We can explore this question using LEGWORK.
Let us compute the horizon distance for a grid of chirp masses and orbital frequencies. First, we can recall that the signal-to-noise ratio of a source is inversely proportional to its distance from LISA and so we can find the horizon distance, D hor as where ρ(D) is the signal-to-noise ratio as some distance D and ρ detect is the threshold above which we consider a source to be detectable. For the purpose of this example we will set ρ detect = 7. We can then, as a function of chirp mass and frequency, determine the maximum dis- Horizon Distance Figure 2. The horizon distance for circular, stellar-mass binaries in a 4-year LISA mission. The filled contours indicate the horizon distance for different orbital frequencies and chirp masses. We add white dotted contours at 8, 50, 800 kpc and 40 Mpc to highlight the distances to the centre of the Milky Way, the Magellanic Clouds, the Andromeda galaxy and the nearest groundbased gravitational wave detection (GW170817, Abbott et al. 2017) respectively. The diagonal black line shows the frequencies and chirp masses at which the merger time is equal to the observation time. This line emphasises the sharp contrast for the horizon distance for binaries that merge before the LISA mission finishes observing (to the right of the line). The horizontal black lines indicate the approximate location of some common double compact object types on this plot, with the assumed masses labelled below in solar masses.
tance for which the source of interest is detected based on this SNR threshold. This can be most efficiently accomplished using LEG-WORK's Source class, since LEGWORK can then compute the merger times and signal-to-noise ratio of each source with only two lines of code. We convert the SNRs to horizon distances using Eq. 66 and plot the result in Figure 2. The shape of the LISA sensitivity curve is clearly reflected in the horizon distance. This is because circular sources only emit at a single frequency and thus every feature in the sensitivity curve has a strong effect on the SNR (and therefore horizon distance) for circular sources. We also see that once the merger time of a source is shorter than the LISA mission length (shown by the black dotted line), the horizon distance sharply decreases. This is because if a source merges before the LISA mission concludes, it has less time to accumulate signal and thus has a lower SNR and horizon distance.
One can use Figure 2 to estimate the horizon distance for any circular source of interest. To illustrate this, we add solid black lines to indicate the typical chirp masses of some possible stellar mass gravitational wave sources as well as white dotted lines to show the distances to nearby galaxies as well as the nearest groundbased gravitational wave detection (GW170817, Abbott et al. 2017). For example, we see that a circular NSNS with an orbital frequency greater than a mHz is detectable in the Magellanic clouds.

The role of eccentricity -
The role of eccentricity is important to consider in the detection of gravitational waves with LISA and other space-based detectors, as sources can still have significant eccentricity during their inspiral phase. A high eccentricity has two major effects on the SNR of a gravitational wave source in LISA and we can investigate these effects using LEGWORK.
Let's consider three hypothetical systems that are identical apart from their eccentricities, e i . We see two effects in the signal-to-noise ratio here. First, increasing the eccentricity from essentially circular to e = 0.6 results in a higher signal-to-noise ratio (ρ = 31.7 → ρ = 50.2). This is because an eccentric binary has enhanced energy emission via gravitational waves (Peters & Mathews 1963). This means that an eccentric binary will not only inspiral faster than an otherwise identical circular binary, but also will always produce a stronger gravitational wave strain. We discuss the enhancement factor and its exact dependence on eccentricity in more detail in Section 3 (see specifically Eq. 6).
The second effect is more intriguing. We see that increasing the eccentricity from e = 0.6 to e = 0.9 results in a relative decrease in SNR (ρ = 50.2 → ρ = 38.8).
The reason for this is that eccentric binaries emit gravitational waves at many harmonic frequencies (unlike circular binaries, which emit predominantly twice the orbital frequency). This leads to the gravitational wave signal being diluted over many frequencies higher than the orbital frequency, where the higher the eccentricity, the more harmonics are required to capture all of the gravitational luminosity (see Figure 3 of Peters & Mathews 1963). Therefore, if the eccentricity is too high, the majority of the signal may be emitted at a frequency to which LISA is less sensitive.
We can illustrate this point with LEGWORK by calculating the SNR at each individual frequency harmonic Figure 3. An illustration of the effect of eccentricity on the detectability of a LISA source. The three sets of points are coloured by their eccentricity and each individual point corresponds to a harmonic frequency, where its height above the curve gives its SNR. We annotate each set of points with its total SNR and overlay the LISA sensitivity curve. The dotted vertical lines indicate the frequency at which the majority of the gravitational wave signal is concentrated.
using the snr module. We plot this distribution of signal over different frequency harmonics with the LISA sensitivity curve overlaid in Figure 3 using LEGWORK's visualisation module. A point is plotted for each harmonic of each source that has an SNR greater than unity and such that its height above the sensitivity curve corresponds to its SNR.
From Figure 3, we can better understand why a source with e = 0.9 has a lower SNR than the same source with e = 0.6. From the dotted lines, we can note that the signal from the e = 0.9 source is concentrated at a frequency of around 40 mHz. The LISA sensitivity at this point is much weaker than the 6 mHz at which the e = 0.6 source is concentrated. Therefore, although the strain from a more eccentric binary is stronger, the SNR is lower due to the increased noise in the LISA detector.
Overall, we can therefore conclude that for LISA sources of this nature, higher eccentricity will produce more detectable binaries only if the orbital frequency is not already at or above the minimum of the LISA sensitivity curve.
Another consideration for more massive binaries is whether the increased eccentricity will cause the binary to merge before the mission ends, which would cause a significant decrease in signal-to-noise ratio. We can also use LEGWORK to find how the merger time of a source varies with frequency and eccentricity over a grid of sources.
We plot the results of this calculation in Figure 4. This plot shows that, for most eccentricities, the merger time is largely determined by the orbital frequency. However, for high eccentricities (e > 0.8), the eccentricity leads to a significant reduction in the merger time. Additionally, we can see that any binary that is to the right of the dashed line on this plot at the start of a 4-year LISA mission, would merge before the mission ended. Therefore, if increasing a binary's eccentricity moved it to the right of this line, its SNR will decrease significantly.

Comparing gravitational wave detectors -
It may also be useful to consider how changing the specifications of the LISA detector, or using a different detector entirely, could affect the SNR of a particular source. LEGWORK is capable of adjusting the LISA mission specifications or using a different sensitivity curve and thus we can use it to explore these differences.
As a first step, we can use LEGWORK to plot a series of sensitivity curves in Figure 5. We show the LISA sensitivity curve for the default 4 year mission length but also illustrate how the curve changes for shorter mission lengths. At 0.5 and 2 years, we see a stronger noise level around 3 mHz as a result of the increased Galactic confusion noise. This noise decreases with increasing mission length since more individual foreground sources can be resolved and thus removed from the confusion noise. We also see that using an approximated response function smooths out the sensitivity curve at higher frequencies. Finally, the TianQin curve is higher than the 4 year LISA curve until around 5 mHz, beyond which it has a lower noise level.
Although comparing the sensitivity curves would suffice for a stationary and circular source (since it would remain at a single frequency), LEGWORK can also be Figure 5. The strain spectral density of the LISA detector with different specifications (Robson et al. 2019) and the TianQin detector (Huang et al. 2020). We show the LISA curve for three different mission lengths and once with an approximate response function.
used to see how the relative SNR between two detectors changes over a range of eccentricities and frequencies.
Using LEGWORK, we can compute the SNR of a grid of sources (spanning a range of frequencies and eccentricities) for both detectors. In Figure 6, we show the ratio of the SNR in LISA to the SNR in TianQin. This plot shows that for circular binaries, the SNR of the source in LISA is stronger up to an orbital frequency of approximately 2.5 mHz, beyond which the SNR of the source is stronger in TianQin. This transition frequency becomes lower with increasing eccentricity as one would expect since eccentric sources emit more at higher harmonics and thus higher frequencies.
4.5. Track SNR of a binary over time -As a binary inspirals its orbital frequency and eccentricity change and this in turn affects the SNR of the binary. For this use case we will demonstrate how LEG-WORK can be used to track the evolution of these parameters and pinpoint the moment at which a binary becomes detectable.
Let's consider a binary with the following initial parameters m 1 = m 2 = 15 M , d = 20 kpc, f orb,i = 3 × 10 −5 Hz, e i = 0.5, and use LEGWORK's evol module to evolve the system until 100 years before its merger with 1000 linearly spaced timesteps, recording the eccentricity and frequency at each timestep.
We plot this evolution of the eccentricity and frequency in the top two panels of Figure 7 as a function of the time before the merger. We see that the binary  Figure 6. The ratio of the SNR in LISA (for a 4-year mission) to the SNR in TianQin. The dashed line indicates the transition at which the SNR is equal in both detectors. We annotate the regions in which either detector has a higher SNR and also annotate the mass and distance of each source in the grid.
circularises and increases its orbital frequency as it inspirals as we would expect.
To take this a step further, we can consider the binary at each timestep to be a separate source with the current eccentricity and frequency. It is then trivial to use LEGWORK to calculate the SNR for each of these 'sources' and thus attain the SNR evolution, which we plot in the last panel of Figure 7.
We see that the SNR increases monotonically over time and sharply increases as the binary approaches its merger. Around 1 Myr before the merger, the SNR reaches the detection threshold and thus could then be seen by a 4-year LISA mission.
Note that LEGWORK could also be used in this way to find the SNR of any system in which the orbital evolution is known. Thus for a triple system or a binary experiencing gas drag (as long as the evolution of the eccentricity and orbital frequency is known) LEGWORK is entirely capable of calculating the SNR evolution. 4.6. LISA verification binaries - Kupfer et al. (2018) presents the LISA verification binaries, a collection of known binary systems that have gravitational-waves that are strong enough to be detected by LISA. In LEGWORK, we provide easy access to this data through the VerificationBinaries class, which is a subclass of Source. This means that the class works identically to Source, but it has the verification . The evolution of a binary system's eccentricity (top), orbital frequency (middle) and SNR (bottom). Each panel is annotated with the constant parameters of the system. The SNR is calculated for a 4-year LISA mission. We annotate in a line the bottom panel at SNR = 7 to highlight the moment at which the source becomes detectable.
binary data (such as their masses and orbital frequencies) pre-loaded into the variables. In addition to the base variables, this class also includes the designation of each binary and the SNR that Kupfer et al. (2018) computed. We note that this SNR differs from the SNR that LEGWORK gives for each binary. This is because Kupfer et al. (2018) run a full, detailed LISA simulation for these binaries, whilst we follow the orbit-averaged approach (see Section 3.3.3), since the former would be intractable for any large number of sources.
As an example of how you could use this data, we use the VerificationBinaries class to plot the sources on the LISA sensitivity curve with the SNR calculated by Kupfer et al. (2018) for a 4-year LISA mission and colour the binaries by their primary mass. We see that the verification binaries tend to be detected with frequencies between 1-10 mHz and have masses estimated to be between 0.1-1 M .

CONCLUSION & SUMMARY
We have presented LEGWORK, a package designed to aid in calculations for stellar-origin binary sources of mHz gravitational-wave observatories like LISA. We outlined the implementation of orbital evolution due to gravitational-wave emission, gravitational wave strain, SNR, and visualisation modules and provided a detailed derivation for each of the equations required for each module. Finally, we provided several use case examples for how LEGWORK can be used to better understand the detectability of compact-object binaries.
We are grateful to Stas Babak, Floor Broekgaarden, Tom Callister, Will Farr, Yi-Ming Hu, Stephen Justham, Valeria Korol, Mike Lau, Tyson Littenberg, Ilya Mandel, Alberto Sesana, Lieke van Son, the CCA GW group, the BinCosmos group and the COMPAS group for stimulating discussions that influenced and motivated us to complete this project. We thank the BinCosmos group for testing an early version of the package and providing useful feedback. In particular, we thank Lieke van Son for her innovation in inventing the name LEGWORK! TW thanks Floor Broekgaarden for first suggesting that he investigate LISA and the derivation of the SNR calculation. KB thanks Shane L. Larson and Kyle Kremer for several conversations about the SNR derivation. We additionally thank the anonymous reviewer for their useful comments, which helped to improve the quality and clarity of this paper.