The initial detections of gravitational waves (GWs) emitted during the coalescences of black holes and neutron stars [1, 2] have confirmed a fundamental prediction of general relativity and opened a new observational window into the universe in the same breath. Further detections by the LIGO and Virgo observatories have pushed the number of observed GW signals into the decades [3] , and hence set off the era of GW astronomy. This success has ignited the desire for detections that allow for better estimation of the information buried in the signals that the GW observatories record, thus allowing us to reveal new phenomena and novel physics as well as shed light into greater depths of the cosmos.
Such desires pose challenges on many essential components in the figurative GW detection machinery, with the development of new waveform (WF) models and the planning of future detector designs being two important tasks to tackle. The former is imperative to lower the systematic biases in what kind of signals the observatories will be able to detect, to increase the information that can be extracted from the detector output, and possibly to lower the computational cost of WF model evaluation. The inclusion of higher order spherical harmonic modes beyond the quadrupole [4, 5] is one such example illustrating the impact of improved WF models as it enhances the detectability of binaries whose orbits are highly inclined with respect to the line of sight and allows for an improved parameter estimation by lifting the degeneracy between luminosity distance and inclination angle.
The detector challenge is even more fundamental: new, ground-breaking analyses require 'richer' data to forge through. The current, so called second-generation (2G) of GW detectors, LIGO Hanford and Livingston [6] as well as Virgo [7] , will be joined by two more detectors, KAGRA [8] and LIGO-India [9] , and further sensitivity improvements are planned to go in effect in the years to follow [10] . Nevertheless the limits of these facilities to address future, scientific goals are already * sub284@psu.edu apparent [11, 12] . Thus the GW community has been pursuing the conceptualization of the next generation of detectors. The proposals range from the Voyager efforts, envisioning a multitude of upgrades to the LIGO facilities [13] , to completely new, third-generation (3G) observatories, namely the Einstein Telescope [12] and the Cosmic Explorer effort [14] . In addition, efforts are under way to search for GW from space, with the Laser Interferometer Space Antenna (LISA) [11] .
Independently of what WF models shall be developed and adopted, or which facilities continue or start operating, the GW community will need tools to benchmark the science output of any combination of these two components in a quick and efficient way. One such tool is the Fisher information formalism (FIF) [15, 16] which has proven to be a viable method for GW benchmarking. Its major caveat is its dependence on high signal-to-noise ratios to give reliable results, thus making it mostly suitable to exceptional events in current detectors. Nevertheless, the formalism is very useful and widely applicable in studies entailing future generations of detectors. For an extended review on the caveats of the FIF we refer to [17] .
In this work, we present a new Python package, gwbench, that implements the FIF in an easy-to-use manner and further provides means to compute and access a host of quantities that are necessary for benchmarking. For a selection of GW detectors, a so-called network, and given a WF model, gwbench can compute: plus and cross polarizations of the WF, detector power spectral densities, antenna patterns, location phase factors, detector responses, detector and network signal to noise ratios, and measurement errors in the WF model parameters from the FIF. Two particularly substantial features are the inclusion of WFs from the LSC Algorithm Library (LAL) [18, 19] and the capability to include the effects of Earth's rotation in the detector antenna patterns.
The paper is structured as follows: Section II introduces the Cosmic Explorer trade study which stimulated the development of this package. In Section III we briefly review the FIF and aforementioned core concepts utilized by gwbench. Next, we show a straightforward example application of gwbench to showcase its functionality in Section IV and proceed to present a few validation tests of the core methods of gwbench in Section V. In Section VI, we describe the application of the code in two previous papers and conclude this work in Section VII.
The Cosmic Explorer Project 1 aims to develop among other things the science case for a next-generation GW detector, the Cosmic Explorer (CE). CE is a promising US proposal for a ground-based detector design which is targeting to answer a wide variety of scientific questions [20] [21] [22] [23] [24] as part of a global network of GW observatories. These questions are grouped under three overarching science goals: (i) dynamics of dense matter and extreme environments, (ii) black holes and neutron stars through cosmic time, and (iii) extreme gravity and fundamental physics. The questions encompass a broad range of fields relevant to GW astronomy including but not restricted to neutron star physics, cosmology, tests of general relativity, and the feasibility of multi-messenger observations. gwbench was developed to perform an exhaustive trade study to gauge CE's performance in reaching the set science goals. The trade study examines some 454 global GW detector networks combining 50 CE configurations with plausible background detectors such as upgraded 2G facilities or the Einstein Telescope. Besides the number of networks, the computational requirements are also driven by the list of science questions to be answered.
Hence, we decided to implement the FIF in gwbench, as it allows us to perform the required tasks effectively for many events and different waveform models in each network of interest. The measurement error estimates computed from the FIF are processed to obtain pre-defined performance metrics which in return allow us to evaluate the science potential of the various network configurations. Such metrics are for example the signal-to-noise ratios, 90%-credible sky areas, expected event rates, and measurement error estimates on WF parameters, but also on inferred quantities such as neutron star radii.
As such, gwbench has an important role to play for the development of CE's science case. The broad range of science goals that can be examined with the code showcases that gwbench is a powerful tool that has a wide variety of applications in gravitational-wave astronomy. We present two of our recent studies that were enabled by gwbench in Section VI.
gwbench implements a variety of standard concepts and aims to make them easily accessible to the user. We briefly review them [25] in this section.
Waveform polarizations: Gravitational waves are emitted by a host of astrophysical and cosmological phenomena, with the coalescence of compact binaries as the most prominent example. The incoming GWs can be decomposed into two polarizations, plus h + (t; Θ GW ) and cross h × (t; Θ GW ), and carry an abundance of information about their source systems, encoded in the parameter vector Θ GW . This vector depends on the WF model. In case of the coalescence of two black holes the standard set of such parameters, Θ GW = {m 1 , m 2 , χ 1 , χ 2 , D L , ι, t c , φ c } contains the companion masses m 1 , m 2 and dimensionless spin vectors χ 1 , χ 2 , the binary's luminosity distance D L , and the inclination angle ι of its orbital plane with respect to the line of sight, as well as two integration parameters: the time t c and phase φ c of coalescence.
In the remainder of this manuscript, we will use the Fouriertransforms of the WF and other quantities in the frequencydomain:
(3.1)
Antenna patterns and location phase factor: The quadrupolar nature of GWs results in sky-position dependent responses from detectors to incoming waves. The responses are governed by the so-called antenna pattern functions F + and F × which are functions of the source's sky position-right ascension α and declination δ-and polarization angle ψ. Besides, the relative motion between the source and detector results in a frequency-dependence in the Fourier-domain; albeit being negligible for signals that are short compared to the time scale of the motion [26, 27] .
Another group of important quantities that encode the spatial distribution of detectors in networks are the location phase factors F lp which shift the phase of the detector responses appropriately for each detector:
wherer and d are the unit vector pointing to the source and the position vector of the detector, respectively, from the center of Earth. α d and β d are the detector longitude and latitude, respectively. Detector response: The aforementioned detector response H is the GW strain in the frequency-domain measured by a detector, combining the WF polarizations with the detector's antenna patterns and location phase factor. It differs from the actually recorded detector output s(t) = h(t) + n(t) in the time-domain or S ( f ) = H( f ) + N( f ) in the frequency-domain, which also contains contributions from the detector noise n(t) or N( f ), respectively. The detector response is given by
where Θ = {Θ GW , α, δ, ψ}. Power spectral density: The sensitivity of a detector determines its ability to measure the strain H. It is governed by systematics due to detector noise n(t) and characterized by the noise auto-correlation κ = n(t 1 ) n(t 2 ) 2 . In the following, we assume the noise to be stationary and Gaussian, with mean n = 0. The former assumption ensures that κ only depends on the time-difference t * = t 1 − t 2 . Thus, the one-sided power spectral density (PSD) S n -the Fourier-transform of κ for positive frequencies and thus its frequency-domain counterpart-is given by
S n ( f ) indicates the sensitivity of the detector at different GW frequencies.
Noise-weighted scalar product and signal-to-noise ratio: The existence of noise in the detector output affects the detection and measurement of the GW signal strain, obscuring the difference between various detector responses to some level. One way to assess the overlap of these responses in a manner that takes the detector's sensitivity into account is by defining a noise-weighted scalar product [25] ,
where H, G are two detector responses, e.g. for two different WF models or the same model evaluated at different parameters. Given this scalar product, we can calculate the signal-tonoise ratio (SNR) ρ of a given detector response H for a given detector PSD S n as
indicating the relative loudness of the detector response compared to the noise in the detector output [25] . Fisher information formalism: The FIF or Fisher analysis is well-established in GW data analysis and discussed in detail in [15, 16] . Here, we give a brief summary of how it allows us to estimate parameter error bounds.
Since we assumed that the noise n = s − h behaves like Gaussian with zero mean, the same holds for its Fourier transform N = S − H and thus the probability of the noise can be written as
where Θ and p 0 are the parameter vector of the detector response H and the prior on these parameters, respectively. Then the value of Θ at peak probability is a good estimate of the true value Θ * for the given GW signal S , assuming a large SNR.
The probability is the greatest when the exponential is the largest 3 . Around the maximal value, the exponent E = S − H, S − H can be expanded as
can be simplified for large SNR where the first-order derivatives of H dominate the second-order ones [15] to give the Fisher information matrix (FIM) Γ 
Thus the assumption of Gaussianity motivates the association of the FIM to be the inverse of the covariance matrix Σ ≡ Γ −1 . The diagonal and off-diagonal elements of Σ denote the variances and covariances of the parameters, respectively, due to the uncertainty introduced by the detector noise and give 1σ-error estimates via σ Θ i = √ Σ ii . Detector network and associated quantities: The above definitions of the SNR ρ and FIM Γ depend on quantities like the detector response and PSD and, ergo, are inherently detector-specific. In a network, detectors work in tandem to improve detectability and parameter estimation which can be captured by defining network versions of the SNR and FIM 13) where N d is the number of detectors in the network and ρ n and Γ n are the SNR and FIM of the n-th detector, respectively. Consequently, the network covariance matrix and error estimates follow from Γ net as defined before. Condition number of a matrix: If the covariance matrix is obtained via numerical inversion, the FIM needs to be wellconditioned for the task. This can be checked with the condition number c Γ = e max /e min , where e max , e min are the largest and smallest eigenvalues of Γ. If the condition number exceeds a critical value-a conservative choice for double-precision is c crit = 10 15 [28] -the inversion result cannot be trusted generically and should thus be disregarded.
The installation instructions for gwbench and the source code are publicly available on GitLab: https://gitlab. com/sborhanian/gwbench. The package consists of 33 modules which implement the concepts discussed in Section III as well as several WFs models. In this section, we will demonstrate with the aid of a few examples how to use the main module of gwbench, network, to perform GW benchmarking. These examples are available as Python scripts at https://gitlab.com/sborhanian/gwbench/-/tree/ master/example_scripts. network contains the Network class, see Appendix A, which is designed to handle the concepts introduced in Section III-detector initialization (antenna patterns, location phase factors, and PSDs), WF evaluation, and both SNR and error benchmarking-by conducting and facilitating the tasks of the other modules in the package. We will refer to the appendices as needed to discuss these modules and aspects of the package.
We begin with a step-by-step showcase of the basic usage of the Network module following the example script quick start.py. In this example we estimate the SNR, measurement errors, and 90%-credible sky area for a a signal with the same masses as the first detected GW event GW150914 [1] , if it was measured by three detectors at the LIGO sites in Hanford, WA and Livingston, LA and the Virgo site in Cascina, Italy. Each detector is set to the aLIGO detector sensitivities [6] and we perform the Fisher analysis with the WF model TaylorF2, see Appendix B.
We begin by importing the numpy and gwbench.network modules and specify the network spec which tells the Network which Detector instances to initialize: ['aLIGO H','aLIGO L','aLIGO V'] will result in three detectors with aLIGO technology, located and oriented like the Hanford, Livingston, and Virgo detectors (refer to the Appendices C and D for other sites and sensitivities available in gwbench). import numpy as np from gwbench import network
While the Network class is the main interaction point for the user, single-detector tasks are relayed to the Detector class. Appendix A compares these two classes to show their similarities and how they synergize to conduct the GW benchmarking together.
Next, we choose the WF model of interest, TaylorF2, and initialize the Waveform object inside our Network. Waveform is another important class that handles WF selection and evaluation. Appendix B contains a brief description of this class and the full list of available WF models in gwbench, at the time of writing this paper.
Finally, we need to set a few Network variables like the frequency array and injection parameters used for WF evaluation. Further, the code needs to know with respect to which of these parameters to take derivatives for the FIM. Besides, we can set the Network to convert the derivatives ∂ p and errors σ p for certain parameters p into the cosine (∂ cos(p) and σ cos(p) ) or logarithmic (∂ log(p) and σ log(p) ) versions. The final option tells the code whether to take Earth's rotation into account for the computation of the antenna pattern functions. The high frequency cutoff of 67.6 Hz corresponds to the GW frequency at the innermost stable circular orbit of a compact binary coalescence with total mass M = M/η 3/5 = 65 M for the example injections with chirp mass M = 28.1 M and symmetric mass ratio η = 0.247. The Network is prepared for benchmarking and we calculate the WF polarizations and their derivatives with respect to the chosen parameters with the following function calls.
net.calc_wf_polarizations() net.calc_wf_polarizations_derivs_num() The first line only computes the plus and cross polarizations in the frequency-domain, net.hfp and net.hfc, while the second obtains a dictionary net.del hfpc containing the derivatives of the two WF polarizations, in addition to that. The structure of these derivative dictionaries is detailed in Appendix E.
Next, we setup antenna patterns, location phase factors, and PSDs for each Detector in the Network . 
We note that the Detector automatically truncates the frequency array, if necessary, to the range dictated by the given detector's PSD. This truncation translates to all consecutive computations which depend on the frequency, e.g. in the detector response or the location phase factor, and is necessary to avoid array-length mismatches within a Detector instance. On the other hand, it can result in a mismatch when comparing frequency-dependent quantities of two Detector instances or between the Network and one of its Detectors.
With antenna patterns and location phase factors prepared, we can calculate the detector responses and their derivatives, analogously to the case for the WF polarizations. With the previously computed PSDs, we obtain readily the detector and network SNRs.
The function calculates both ρ and ρ 2 for the total network and each detector in the network. The values are saved in the respective Network and Detector instances.
Finally, we compute the error estimates on the parameters specified in deriv symbs string and the estimates of the 90%-credible sky area.
net.calc_errors() net.calc_sky_area_90() net.print_detectors() net.print_network() calc errors contains several steps: It starts by calculating the FIMs for all the Detector instances and then the Network one from those. Next, it computes and saves the condition number c Γ for each FIM Γ, and checks if the matrices are wellconditioned (c Γ < 10 15 ). If this criterion is fulfilled, the code will proceed to invert each FIM to obtain all the covariance matrices Σ and consequently the error estimates. Finally, it calculates the inversion accuracy = ||Γ·Σ− I|| max , where I and || · || max are the identity and maximum matrix norm, respectively.
Analogously to the SNR computations, the Network and Detector instances save their respective quantities. The final two function calls simply print the contents of all the Detector instances inside the Network as well its own contents.
Further remarks: For the sake of simplicity, we omitted non-default function arguments in the example above. The Network class is designed to perform 'default'-benchmarking after setting all the necessary variables via net.set wf vars and net.set net vars without the need to pass further variables during the consecutive function calls. Nevertheless, certain functions give the user the option to tune the functionality of the code. In the following, we will present these options for the main GW benchmarking functions showcased in the example.
net.setup ant pat lpf psds: This function takes the following three arguments (default values) F lo (-np.inf), F hi (np.inf), and psd file dict (None). The first two are user-settable frequency cutoffs which truncate the frequency array, if they lie within the PSD-intrinsic frequency range. The third arguments allows the user to specify PSD files that are not part of gwbench. It has to be passed as
where the det key are the exact keys that each Detector contains, the file are the relative paths to the PSD files, and the booleans bool tell the code whether the file contains PSDs S n (bool=0) or amplitude spectral densities √ S n (bool=1). A caveat of using external PSD files is that the Network and Detector labels will not reflect this choice and will need manual changes. Fortunately, these can be done quite easily, if needed, in run time.
gwbench can naturally be used to just compute the antenna patterns, location phase factors, and PSDs of detectors as shown in the example script compute ant pat lpf psd.py.
net.calc det responses derivs num, net.calc wf polarizations derivs num: These functions calculate the derivatives of the detector responses and WF polarizations numerically via the numdifftools package [29] . The numeric derivatives require four special arguments: the step size step, the differentiation scheme method, the differentiation order order, and the derivative order n. The respective default values in gwbench are 1e-7, 'central', 2, and 1. We refer to the documentation of numdifftools for further information. gwbench further allows the use of symbolic derivatives for certain waveforms via the sympy package. Their evaluation is faster, but requires a few extra steps compared to the numeric ones. Their usage is showcased in Section IV C, while Section V A compares the output of the numeric and symbolic differentiation methods.
net.calc snrs det responses, net.calc errors, and net.calc sky area 90: These three functions have a common argument only net which tells the code whether to calculate the respective quantities just for the Network or also for each Detector within. The default is 0, computing both network and detector quantities.
Besides, net.calc errors takes two more arguments (default values) cond sup (1e15) and by element (0). The first determines the condition number supremum c * that the code uses to assess whether a computed FIM Γ is wellconditioned (c Γ < c * ). Passing None sets it to infinity, effectively making every FIM well-conditioned. The second argument tells the code whether to obtain a single inversion error for the entire FIM via the maximum norm (default) or n errors i . n is the dimension of the FIM while each i is the maximum of the elements of the i-th row and column combined. The single inversion error is equal to the maximum of the n errors.
Compute dictionaries of lambdified sympy derivatives: While numeric derivatives are computed and evaluated in the same, computationally expensive step, symbolic derivatives via sympy are first computed generically and then evaluated for a given set of arguments. The latter is faster than the numeric approach as shown in Section IV B, but requires the generation of the so-called lambdified derivative functions ahead of time. In the following, we will showcase an example script that can handle this step: generate lambdified functions.py.
The user has to specify the WF, the derivative variables, the locations of interest, and whether to take Earth's rotation into account. By default, the script is setup to (i) use the WF model TaylorF2, (ii) take derivatives with respect to chirp mass M, symmetric mass ratio η, luminosity distance D L , time and phase of coalescence t c , φ c , inclination angle ι, right ascension α, declination δ, and polarization angle ψ, (iii) use all locations available in gwbench, and (iv) include Earth's rotation.
Next, the code initializes a Waveform object and checks whether the derivative variables actually form a subset of the WF and antenna pattern arguments. The actual computation is then handled by dr.generate_det_responses_sym( wf,deriv_symbs_string,locs=locs,use_rot=use_rot)
which takes the Waveform, the derivative variables, the locations of interest, and the rotation setting. The function generates two sets of dictionaries containing different generic derivative expressions stored as lambdified sympy functions. The first type contains the derivatives of the WF polarizations while the second stores the derivatives of the detector responses for all detector locations implemented in gwbench. These dictionaries are then saved in a folder named lambdified functions using the dill package. The output files contain one dictionary each and the file names are structured as 'par_deriv_WFM_{wf_model_name}_ VAR_{variables}_DET_{location}.dat'
where {wf model name} and {variables} are set according to the values passed to generate det responses sym. {location} is replaced by pl cr in case of the WF polarizations. Symbolic derivative evaluation for Fisher analysis: sym gw benchmarking.py showcases how to load and evaluate the symbolic derivatives generated with generate lambdified functions.py and how to use the injections module to generate random samples of input parameters for the detector response evaluation (see Appendix F). The script is mimicked by num gw benchmarking.py which differs by computing the derivatives numerically; it can be seen as a more generic expansion of the example shown in Section IV B.
Most of the script simply prepares the benchmarking and hence we focus on the usage of symbolic derivatives, in order to avoid repetition. As before, the Network is the central hub and requires the same setup as shown in Section IV B to prepare the Network and Detector instances. In contrast to the earlier numeric example, the function call net.calc_det_responses_derivs_num() is replaced by net.load_det_responses_derivs_sym() net.calc_det_responses_derivs_sym() which first load the generated lambdified derivative functions and then evaluate them. This is all that is needed besides ensuring that the variables wf model name, wf other var dic, deriv symbs string, and use rot are set to the same values as during the generation of the lambdified derivatives.
In the following, we present benchmarking results obtained with gwbench. We compare numeric and symbolic derivatives and show their effects on the error estimates. We further present error estimates from various WF models for particularly interesting real events, namely GW170817 [2] , GW190521 [30] , and GW190814 [31] .
A. Comparison of numeric and symbolic derivative methods gwbench contains both numeric and symbolic differentiation methods. The former can be used with all WF models, while the latter only works with some of the implemented WFs, as specified in Appendix B. In order to compare the two approaches for the purposes of GW benchmarking, we will compute the partial derivatives of the detector responses H and the error estimates for two WF models, TaylorF2 (tf2) and TaylorF2 with tidal effects (tf2 tidal).
Both WFs share the same set of input parameters, (M, η, χ 1,z , χ 2,z , D L , t c , φ c , ι, α, δ, ψ), with the tidal WF further depending two tidal deformability parameter combinations (Λ, δΛ) [32] FIG. 1. The plots show the relative residuals r, see Equation (5.3), between the numeric and symbolic partial derivatives of the detector responses for 11 (12) input parameters and WF model tf2 (tf2 tidal). The ∼ 20 very poor residuals in the chirpmass derivatives correspond to ill-conditioned Fisher matrices and are therefore disregarded for further analysis.
where the component tidal deformabilities Λ 1 and Λ 2 are defined as Λ i = 2 3 k 2 C 5 i with Love number k 2 , compactness C i = m i R i , mass m i , and radius R i . The partial derivatives are taken with respect to all available parameters except δΛ for the tidal WF. The detector responses are evaluated for three detectors, each set to a 40km compactbinary optimized CE1 PSD, see Appendix D.
We performed this comparison for a total of 1000 random sets of input parameters for each WF. In both cases we chose the injections to be uniformly distributed in co-moving volume up to a redshift of z = 0.5 as well as over sky positions and orientation angles. The masses and spins are also sampled uniformly for both WFs albeit the ranges are chosen to resemble binary black hole or binary neutron star systems, respectively: 5 M ≤ m 1 , m 2 ≤ 100 M and −0.75 ≤ χ 1,z , χ 2,z ≤ 0.75 for tf2 and 1 M ≤ m 1 , m 2 ≤ 2 M and −0.05 ≤ χ 1,z , χ 2,z ≤ 0.05 for tf2 tidal. Figure 1 shows the relative residuals of all partial derivatives of the detector responses for both waveforms and all three detectors. The residuals are defined via
where i runs over the input frequency array and x ∈ {ln(M), η, χ 1,z , χ 2,z , ln(D L ), t c , φ c , cos(ι), α, cos(δ), ψ, ln(Λ)}. The inclusion of tidal effects results in two trends for the relative residuals of the partial derivatives in Figure 1 : The errors of the derivatives with respect to the so-called intrinsic parameters, (M, η, χ 1,z , χ 2,z , t c , φ c ,Λ, δΛ), worsen but are also more tightly constrained, while the derivatives with respect to the extrinsic parameters (D L , ι, α, δ, ψ) appear to not be affected. The only exceptions are the chirpmass and symmetric mass ratio derivatives, with the former showing improved residuals with tidal effects and the latter showing more scatter. The reason for the discrepancy in the behavior of the intrinsic and extrinsic derivatives is found in the correlation of the intrinsic parameters to each other, thus worsening the numeric derivatives when adding another intrinsic parameter, whereas the extrinsic parameters stay unaffected.
Overall, the partial derivatives appear to be well-behaved with relative residuals well below 10 −4 . The tidal parameter derivatives perform the worst, but stay below an error of 10% for most of the samples. Finally, the ∼ 20 very poor residuals in chirpmass derivatives actually correspond to ill-conditioned Fisher matrices and are therefore disregarded for further analysis. Figure 2 shows the relative residuals
The plots show the relative residualsr, see Equation (5.4) , between the 1σ-error bounds, estimated via the FIF from numeric or symbolic derivatives, for 11 (12) input parameters and WF model tf2 (tf2 tidal).
of the 1σ-errors bounds estimated from the FIF using the numeric and symbolic derivatives of the detector responses. While the relative residuals on the parameter error estimates from Figure 2 do not show the actual discrepancy between the two derivative methods, they do show their effect on the final product of the Fisher analysis. It is important to note that the number of injections shown has decreased due to some parameter combinations yielding ill-conditioned Fisher matrices for which the further analysis is aborted. In fact, it appears that the inclusion of tidal effects has a stabilizing effect for the FIF, resulting in many more well-conditioned Fisher matrices. We remark that the mentioned differences in the behavior of the residuals between tf2 and tf2 tidal could be simply a symptom of the different input parameter sets which we leave for a future study.
It appears that the correlations between intrinsic or extrinsic parameters smudges the differences in the numeric and symbolic derivatives within each of these parameter groups and the inclusion of tidal effects has little influence on the residuals. The only exception is presented in the case of the luminosity distance error bounds where the tidal WF resulted in improved errors. Overall, the relative residuals between the parameter error bounds, estimated with numeric and symbolic derivative methods, stay under 1% and mostly under 0.1%. Thus, we can use either approach confidently to perform the FIF.
B. Application to real events GW170817 and GW190814 are among the loudest events observed by the LIGO and Virgo detectors so far with SNRs ρ GW170817 = 32 and ρ GW190814 = 25, respectively. These are exceptional candidates to check the FIF and GW benchmarking against. GW190521 is another interesting event representing signals of average loudness with ρ GW190521 = 14 and very short duration as well as heavy binaries with merger-and ringdown-dominated signals. Neither of these events have well determined sky positions and binary orientations. We applied gwbench to 1000 realizations of the various angles while we fixed the masses, spins, and luminosity distances of the systems to the reported median values [2, 3, 30] . All three events were detected by three detectors located at the LIGO Hanford, LIGO Livingston, and Virgo sites. Since the actual sensitivities of these detectors during the detections are determined by the facilities, we used the respective PSDs for this analysis [33] and ran the FIF for 4 different WF models: tf2 tidal for GW170817, IMRPhenomD and IMRPhenomHM for GW190521, and tf2, IMRPhenomD, and IMRPhenomHM for G190814.
We compare the reported values for the SNR ρ, the fractional error estimates of the chirpmass ∆M/M and luminosity distance ∆D L /D L , as well as the 90%-credible sky area Ω to our benchmarking results. While this comparison contains only four quantities, we still perform the FIF with respect to 11 (GW190521 and GW190814) or 12 (GW170817) parameters, in order to not underestimate the error bounds. These 11 and 12 parameters are the same sets as in Section V A above for tf2 and tf2 tidal, respectively. Besides, in order to allow for an easier comparison with the measured values we converted the fractional errors from the standard 1σ-output of the FIF to the respective 90% credible bounds.
The benchmarking posteriors for the 1000 realizations are presented in Figure 3 indicating that they capture the measured errors well for at least one of the WF models per event. The figure shows that tf2 tidal performs really well for an inspiral dominated, binary neutron star signal like GW170817. This does not hold for signals from the heavier binaries. For such systems, the error estimates from tf2 do not capture the measured values due to missed post-inspiral contributions (GW190814), while the FIF completely failed for a strongly merger-and ringdown-dominated signal (GW190521). In fact, for both of these systems it was necessary to use IMRPhenomHM, a WF with higher modes, and not just IMRPhenomD in order to fully capture the measured errors. These findings are in line with the WFs used for the actual parameter estimation analyses of these events and they show that gwbench provides a consistent tool for the purposes of GW benchmarking.
gwbench has already been applied in two works: [34, 35] . In [34] we examined the possibility to use dark sirens-GW sources without an electromagnetic counterpart-to resolve the tension in the measurement of the Hubble constant. In [35] we studied how multiband observations of GW signals with the space-based LISA and ground-based 3G observatories will open up new possibilities to conduct multiparameter tests of general relativity in the line of [36] .
A. "Dark sirens to resolve the Hubble-Lemaître tension"
We examined the prospect to use dark sirens to perform precision measurements of the Hubble constant and thus provide a GW-assisted method to break the Hubble-Lemaître tension. The proposed method makes use of the potential of future ob-servatories to measure sky positions of close-by dark sirens well enough to pin-point the host galaxy of the GW event. This provides us independent measurements, the distance from the GW signal and the galaxy's redshift from electromagnetic follow-up observations, which allow us to estimate the Hubble constant.
We simulated three populations of dark sirens: one resembled binary black holes from the GWTC-1 catalog [3] , the second focused on the heavy sub-population of such BBHs with component masses m 1 , m 2 > 25 M , while the third was GW190814-like [31] . The GWTC-1 catalog summarizes the detections from the first two observing runs of the LIGO and VIRGO detectors and GW190814 was the most asymmetric binary observed to date.
For each population we performed the Fisher analysis in four networks containing three possible detector technologies of the future ('Advanced+', 'Voyager', and 3G). In order to leverage the sensitivities of such future networks we used the lalsimulation WF IMRPhenomHM [4] which also models subdominant, higher-order spherical harmonic modes beyond the quadrupole. These higher modes improve the distance measurement by breaking the degeneracy between the luminosity distance and the inclination angle of the binary orbit to the line of sight, thus ultimately yielding better measurements of the Hubble constant.
Using gwbench we were able to perform the FIF using numeric derivatives on 10 4 binaries per population. We estimated the number of detections we can expect from the studied networks that yield precise measurements of the dark siren's sky position. Finally, we propagated the error bounds on the luminosity distance measurement to estimate the precision of a Hubble constant measurement from such binaries. We showed that, if redshift errors are well controlled, these dark siren observations will yield ∼ 2%-level measurement precisions already in the Advanced+ era of the current detector generation.
B. "Multiparameter tests of general relativity using multiband gravitational-wave observations"
Parameterized post-Newtonian tests [36] present a crucial resource to perform theory-agonistic tests of general relativity which can probe generically potential deviations of GW signals from Einstein's theory. This class of tests is especially important when verifiable GW predictions from alternate theories of gravity are sparse. The tests are based on modifications to the expansion coefficients in the post-Newtonian approximation to the WFs emitted by compact binary coalescences, such that these deviations vanish if general relativity is correct. Unfortunately, the current tests are handicapped by the strong correlations between the introduced deviation parameters. Therefore, parameterized post-Newtonian tests are limited to single-coefficient modifications for now. We showed that the combination of the information present in the low-frequency LISA and high-frequency 3G bands breaks this degeneracy. Thus such multiband observations allow for true multiparameter tests of general relativity.
In our study we examined stellar mass binary black hole coalescences. Their inspiral signals will last several years in the LISA band (∼ 0.1 − 100 mHz) before they reach the 3G band (∼ Hz-kHz) where they finally merge. Unfortunately, these signals, while extremely loud in 3G networks, are going to be faint to non-visible in LISA, making the 'joint' observation a difficult task. The issue lies in the detection methodology for such quiet signals. The current state-of-the-art are blind, matched filtering searches [37, 38] without a-priori knowledge of the GW signal. Such blind searches require GW template banks which cover a large volume in the parameter space. Thus, the sizes of the banks become computationally infeasible for the signals of interest in the LISA band [39] .
Fortunately, we can leverage the high-fidelity detections with 3G detectors to dig out these quiet signals from LISA's archival data. In contrast to a blind search the 3G detection confirms the existence of the signal we are searching for and further provides some information about it, allowing for the generation of signal-specific template banks that decrease the computational cost immensely. This is enabled in two-ways: A 3G network will estimate all WF parameters except chirpmass and symmetric mass ratio to better precision than LISA could, thus decreasing the parameter dimension of the template bank to two and fixing the other parameters to the 3G measurement. The search volume can be further narrowed down to within the error bounds estimated from the 3G detection for both chirpmass and symmetric mass ratio.
Using gwbench, we tackled the question of feasibility to search the archival LISA data following a 3G detection, strengthening our proposal to perform multiband multi-parameter tests of general relativity. For this purpose we generated a set of 500,000 binary black holes up to a redshift of z = 10. Roughly 200 of these binaries should be visible in LISA with an SNR ≥ 4 which we set as the threshold for visibility in the archival search. gwbench allowed us to estimate the SNR and 1σ error bounds on the WF parameters for a network consisting of one Einstein Telescope and two 40 km, CE2 detectors in compact binary optimization.
Given the LISA SNR ρ LISA and 3G errors on chirpmass σ M c and symmetric mass ratio σ η , we used gwbench to calculate the detector-dependent volume V via the metric g [40, 41] 
The inner product in the definition of the metric is taken with respect to the LISA PSD. Finally, the number of templates N is given by the ratio of the total volume to the volume of the error ball around each template for a given minimal match mm
Using mm = 0.95, we found that such targeted searches of LISA's archival data following the detection of a binary black hole merger in a network of 3G detectors should allow for template bank sizes of the order 10 3 -10 4 . This is substantially less than the what was predicted before [39] . Ultimately, gwbench enabled us to check the feasibility of multibanding to perform multi-parameter tests of general relativity by providing the tools to calculate estimates of both the 3G errors bounds and the LISA template numbers.
In this work we presented gwbench, a Python package for fast and straightforward applications of the Fisher information formalism for the purposes of GW benchmarking. The package is written with the usage for detector networks in mind. Thus, it is structured around the Network class which acts as a hub for user interaction and facilitates all the possible computations. The package gives easy access to a variety of GW benchmarking quantities such as the frequency-domain WF polarizations, antenna pattern functions and location phase factors as well as noise PSDs for various detector locations and technologies, frequency-domain detector responses, SNRs, Fisher and covariance matrices, and the product of the FIF, the 1σ error bounds for input parameters of choice.
There are two types of WFs available, both of which support numeric derivatives via the numdifftools package: LAL WFs as well as those implemented gwbench. The latter further support symbolic derivatives via sympy which are faster and more accurate. The inclusion of LAL WF support is very important as it gives access to a host of WF models that are well tested and commonly used in the literature. gwbench contains nine existing, planned, or fiducial detector locations and 23 different detector technologies ranging from aLIGO to various CE sensitivity curves. The addition of user-defined locations and technologies is also possible. Another specialty of this package is the inclusion of the effects of Earth's rotation in the antenna patterns and detector responses, and thus the FIF. This is particularly important for long-lived GW signals such as from binary neutron stars.
gwbench is and has been used for several research projects of which we presented two that are in public already. In both cases, gwbench allowed us to generate and calculate parameter errors and even template numbers for a matched filtered search for large numbers of sources and detector networks. gwbench represents a unique code that implements the FIF and GW benchmarking in an easy-to-use manner for a wide range of problems, WFs, and network options and allows for straightforward application in large scale simulations for improved statistics of the scientific claims. ACKNOWLEDGMENTS I thank Bangalore Sathyaprakash, Anuradha Gupta, Patrick Godwin, and everyone in the Cosmic Explorer Project for helpful comments and discussions during the development of gwbench. I further thank Kevin Kuns, Varun Srivastava, Evan Hall, Matthew Evans, and Stefan Ballmer for providing the sensitivity curves used in gwbench. I also thank all front line workers combating the CoVID-19 pandemic and those maintaining ordinary life without whose support this work would not have been possible. I acknowledge the support by NSF grant PHY-1836779. Computing resources for the development and testing of this code were provided by the Pennsylvania State University.
Appendix A: Network and Detector classes GW detector networks consist, as per name, of a number of detectors that work in tandem to observe GWs to improve jointly signal strength, sky coverage, and parameter estimation. Leaning on this understanding, the Network class is designed as a GW benchmarking hub for a chosen set of detectors. The class structure allows for straightforward handling of the relevant quantities and network methods and thus provides the user with tools to load, calculate, and manipulate quantities like waveforms, SNRs, or errors estimates.
The Detector class is written with similar syntax and methodology as the Network class. This is very natural because network quantities like SNR or error estimates are handled in the same manner as the respective single detector ones. This choice allows for clearer codes and faster development, but also better understanding of both classes, which helps when using Network objects. Thus, • what frequency array and injection parameters to use for function evaluation, and with respect to which parameters to take derivatives for error estimation via FIF,
• what Waveform to initialize,
• which parameters x to convert to cos(x) or ln(x) to estimate errors on the latter, and whether to take Earth's rotation into account.
The bottom section shows the variables which can be loaded and computed with instance methods based on the user-set variables:
• Detector-specific quantities like PSDs, antenna patterns, and location phase factors are loaded and stored in each Detector instance.
• The WF polarizations and the respective derivative dictionaries are detector-agnostic and saved in the Network.
• Each Detector handles the respective detector response and derivative dictionaries.
• The detector and network SNRs (and their squares) are stored in the respective instances, where each Detector also saves the integrand of the square SNR.
• Each detector and the network contain their respective copy of the FIM, its condition number, a wellconditioned flag, the covariance matrix, the inversion error, and an error dictionary containing the 1σ errors extracted from the covariance matrix.
Appendix B: Implemented waveform models and the Waveform class gwbench handles WF selection and evaluation within the Waveform class. The available WFs are listed in Table II together with their labels, their input parameters, and whether they allow for numeric or symbolic differentiation. One of the strong suites of gwbench is the implementation of WF models from LAL in the Waveform class and the capability to use them for the FIF with numeric derivatives. The class contains five instance variables • .wf model name: internal name of the WF model,
• .wf other var dic: extra input parameters, • .wf symbs string: string containing the basic parameters (space-separated),
• .hfpc np: numpy function,
• .hfpc sp: sympy function, and few methods of which two are of particular interest to the user:
• .get sp expr(): load the sympy function to obtain a symbolic expression which can be manipulated,
• .eval np func(f,inj params): evaluates the numpy function for the given frequency array f and injection parameters inj params.
(The injection parameters inj params can be either provided as a dictionary with keys equal to the parameters defined in .wf symbs string or as a list with the same ordering given by .wf symbs string.)
Only the wf model name and wf other var dic are needed to fully initialize a Waveform instance: from gwbench import wf_class wf = wf_class.Waveform( wf_model_name = wf_model_name, wf_other_var_dic = wf_other_var_dic)
The LAL WF models do not provide any of the aforementioned symbolic functionality, enabled by sympy, as this requires a direct implementaion of the WF model to ensure proper sympy function calls.
Appendix C: Implemented detector locations -antenna patterns and location phase factor two modules antenna pattern np.py and antenna pattern sp.py which differ in their intended usage: the former is used to calculate the respective quantities directly whereas the latter is required for the generation of the symbolic derivatives via sympy.
The core concepts for a single detector are demonstrated in [25] , whereas [42] outlines the more intricate methods required for the analysis of a network of detectors where calculations are carried out in a universal frame centered at Earth's center. In either case, F p , F c , and F lp depend on a few signal-specific quantities: right ascension α, declination δ, and polarization angle ψ. In a network, these functions further depend on the arrival time t E of the signal at Earth's center as well as three detector-and location-specific angles: the longitude α d and latitude β d of the detector site as well as the angle γ d of the Figure 5 , as functions of the frequency between 1 and 4096 Hz. The bottom presents the observable redshifts for equal-mass binaries at SNR 100 (solid) and 12 (dotted) as functions of the total binary mass between 1 and 1000 M (in the source frame). The systems are assumed to be optimally orientated for an L-shaped detector with the respective ASD. The SNRs were obtained for the LAL WF IMRPhenomHM. The shaded regions show the respective ranges for events with SNRs between 12 and 100.
detector's y-arm with respect to due East 5 .
All these considerations omit the effects of Earth's rotation on the antenna patterns and location phases factors. These effects are negligible for short-lived, transient signals lasting only a few minutes or less. They cannot be neglected for 3G detectors which will see signals that can be up to a few days 5 Since the Einstein Telescope consists of three V-shaped detectors ET1, ET2, and ET3, that form a equilateral triangle, γ d is fiducially set to correspond to ET1. The fixed geometry allows the calculation of the detector quantities for all three detectors.
long. Hence, gwbench contains the functionality to take Earth's rotation into account which follows the calculations detailed in [26, 27] .
Appendix D: Implemented detector technologies -detector PSDs Figures 5 and 6 show the 23 available PSDs [43] in gwbench. The former is a graphical summary presenting the hierarchy of these functions in terms of detector generation (2G and 3G) and design concepts (Advanced+ vs Voyager vs ET vs CE and more) while the latter shows the amplitude spectral densities √ S n (h) as functions of frequencies. Figure 6 further presents the observable redshift ranges for equal-mass binaries with total mass between 1 and 1000 M (source frame), if observed with SNRs between 12 and 100. Table IV presents  the respective label and key options IV. The available detector PSDs are presented together with the respective labels and internal keys as used by gwbench. CBO stands for compact-binary optimized and PMO for post-merger optimized. The former PSD will be tuned towards increased sensitivity at low frequencies, while the latter is geared towards high frequencies. There are two examples illustrating the nomenclature of the CE configurations.
Appendix E: Derivative dictionaries gwbench obtains the partial derivatives of the WF polarizations and detector responses either numerically or symbolically using the numdifftools or sympy packages, respectively. The former calculates the evaluated derivatives directly which is computationally expensive but applicable to most WFs. sympy on the other hand works with symbolic expressions that require modifications to the WF source code. In return they enable the analytic computation of generic derivative expressions which can be stored in the form of lambdified functions for fast evaluation in run time.
The evaluated partial derivatives are stored in dictionaries which are structured as follows • dim': specifies spin dimension (values: 1 and 3; onedimensional spin sampling disregards 'geom' and samples using Cartesian geometry; further 1 returns zeroes for the x-and y-components to sample aligned injections).
• 'chi lo': specifies the minimum spin value (smallest spin value in each dimension for 'cartesian' geometry or smallest spin magnitude for 'spherical' geometry).
• 'chi hi': specifies the maximum spin value (largest spin value in each dimension for 'cartesian' geometry or largest spin magnitude for 'spherical' geometry).
Angle sampler: The sky location and orientation angles sampling is handled by angle sampler: iotas, ras, decs, psis = injections.angle_sampler( num_injs,seed)
Redshift and luminosity distance sampler: The redshift and luminosity distance sampling is handled by redshift lum distance sampler: zs, DLs = injections.redshift_lum_distance_sampler( cosmo_dict,num_injs,seed)
The cosmo dict contains three definite keys
• 'sampler': specifies the redshift sampling method (values: 'inversion' and 'rejection'; the inversion method is faster, while the rejection method is more accurate, especially for very jagged redshift distributions).
• 'zmin': specifies the minimum redshift.
• 'zmax': specifies the maximum redshift.
By default the code samples redshifts from the redshift distribution astropy.cosmology.Planck15 provided in the astropy package. The user can change this behavior by adding three further keys to the cosmo dict whose values are then passed to astropy.cosmology.LambdaCDM:
• 'H0': specifies Hubble constant.
• 'Om0': specifies the mass density Ω m,0 at the current time.
• 'Ode0': specifies dark energy density Ω Λ,0 at the current time.
Injection parameter sampler: injections further contains a function that combines the four other samplers into one function call: injections CBC params redshift. Hence it depends on the same input as the previous four methods: data = injections.injections_CBC_params_redshift( cosmo_dict,mass_dict,spin_dict,redshifted, num_injs,seed,file_path)
The two additional inputs, redshifted and file path, control whether the sampled masses are already redshifted (1) accordingly to the sampled redshifts or not (0) as well as if the sampled data should be stored in a file at the specified path (default is None), respectively.