GWFAST: A Fisher Information Matrix Python Code for Third-generation Gravitational-wave Detectors

In the most general case, restricting to quasi-circular orbits, the parameters needed to describe the signal emitted by a coalescing binary system are 15 for a binary black hole (BBH), 17 for a binary neutron star (BNS), and 16 for a neutron star black hole binary (NSBH), i.e., ${\boldsymbol{\theta }}=\{{{ \mathcal M }}_{c},\eta ,{d}_{L},\theta ,\phi ,\iota ,\psi ,{t}_{c},{{\rm{\Phi }}}_{c},{\chi }_{1,x},{\chi }_{2,x},{\chi }_{1,y},{\chi }_{2,y},{\chi }_{1,z},{\chi }_{2,z},{{\rm{\Lambda }}}_{1},{{\rm{\Lambda }}}_{2}\}$ (see, e.g., Maggiore 2007), where: ${{ \mathcal M }}_{c}$ denotes the detector-frame chirp mass, η is the symmetric mass ratio, d_L is the luminosity distance to the source, θ and ϕ are the sky position coordinates, defined as θ = π/2−δ and ϕ = α (with α and δ R.A. and decl., respectively), ι is the inclination angle of the binary orbital angular momentum with respect to the line of sight, ψ is the polarization angle, t_c is the time of coalescence, Φ_c is the phase at coalescence, χ_i,c is the dimensionless spin of the object i = {1, 2} along the axis c = {x, y, z}, and Λ_i is the dimensionless tidal deformability of the object i, which is zero for a BH. GWFAST expresses the chirp mass in units of solar masses, M_⊙, the luminosity distance in gigaparsecs, the time of coalescence as a GPS time, in seconds, and all of the angles in radians. Also including a small eccentricity in the orbit, there is one more parameter to consider, e₀, which denotes the eccentricity at a given reference frequency. ³

To use GWFAST, the parameters θ of the events in the catalog to analyze have to be stored in a dictionary, with the same keys as in the following example:

Code example 1. Input parameters for waveform models and FIM calculations in GWFAST.

It is also possible to pass some entries with different commonly used parameterizations, namely:

• (a)  
  the sky position coordinates can be given in terms of R.A., α, and decl., δ, always in radians, in place of θ and ϕ;
• (b)  
  the time of coalescence can be provided as a Greenwich Mean Sidereal Time (GMST) in days, under the entry name 'tcoal', which takes the place of 'tGPS';
• (c)  
  in the nonprecessing case, one can choose the spin parameters χ_s , χ_a instead of χ_1,z , χ_1,z , defined as

  $$\begin{eqnarray}&&{\chi }_{s}=\displaystyle \frac{1}{2}({\chi }_{1,z}+{\chi }_{2,z}),\qquad {\chi }_{a}=\displaystyle \frac{1}{2}({\chi }_{1,z}-{\chi }_{2,z});\end{eqnarray} \tag{ 1 }$$
• (d)  
  it is possible to use the combinations of the tidal deformabilities $\tilde{{\rm{\Lambda }}},\delta \tilde{{\rm{\Lambda }}}$ in place of Λ₁, Λ₂, whose definitions are (Wade et al. 2014)

  $$\begin{eqnarray}\begin{array}{rcl}\tilde{{\rm{\Lambda }}} & = & \displaystyle \frac{8}{13}\left[(1+7\eta -31{\eta }^{2})({{\rm{\Lambda }}}_{1}+{{\rm{\Lambda }}}_{2})\right.\\ & & \left.+\,\sqrt{1-4\eta }(1+9\eta -11{\eta }^{2})({{\rm{\Lambda }}}_{1}-{{\rm{\Lambda }}}_{2})\right],\end{array}\end{eqnarray} \tag{ 2a }$$

  $$\begin{eqnarray}\begin{array}{rcl}\delta \tilde{{\rm{\Lambda }}} & = & \frac{1}{2}\left[\sqrt{1-4\eta }\left(1-\frac{13,272}{1319}\eta +\frac{8944}{1319}{\eta }^{2}\right)({{\rm{\Lambda }}}_{1}+{{\rm{\Lambda }}}_{2})\right.\\ & & +\,\left(1-\frac{15,910}{1319}\eta +\frac{32,850}{1319}{\eta }^{2}+\frac{3380}{1319}{\eta }^{3}\right)\\ & & \left.\times ({{\rm{\Lambda }}}_{1}-{{\rm{\Lambda }}}_{2})\right];\end{array}\end{eqnarray} \tag{ 2b }$$
• (e)  
  if using a waveform model that includes the contribution of unaligned spin components (precessing spins) it is possible to substitute the entries ι, χ_1,x , χ_2,x , χ_1,y , χ_2,y , χ_1,z , χ_2,z with θ_JN , χ₁, χ₂, θ_s,1, θ_s,2, ϕ_JL , ϕ_1,2. These are, respectively, the angle between the total angular momentum and the line of sight, θ_JN , the magnitudes of the spin vectors, χ_i , the angles between the spin vectors and the orbital angular momentum, θ_s,i , the azimuthal angle of the orbital angular momentum and the total angular momentum, ϕ_JL , and the difference in azimuthal angle between the two spin vectors, ϕ_1,2.

A summary of the parameters, their physical symbol, and their name in GWFAST is provided in Table 1.

Table 1. Summary of the Parameters Used in GWFAST to Describe the GW Signal

Parameter Symbol	Parameter Description	Name in GWFAST	Units in GWFAST	Physical Range
${{ \mathcal M }}_{c}$	detector-frame chirp mass	'Mc'	M⊙	(0, +∞)
η	symmetric mass ratio	'eta'	⋯	(0, 0.25]
dL	luminosity distance	'dL'	Gpc	(0, +∞)
θ, ϕ/α, δ(a)	sky position	'theta', 'phi'/	rad	[0, π], [0, 2π]/
		'ra', 'dec'		[0, 2π], [−π/2, π/2]
ι(e)	inclination angle w.r.t. orbital	'iota'	rad	[0, π]
	angular momentum
ψ	polarization angle	'psi'	rad	[0, π]
tc (b)	time of coalescence GPS/GMST	'tGPS'/'tcoal'	s/day	[0, +∞)/[0, 1]
Φc	phase at coalescence	'Phicoal'	rad	[0, 2π]
χi,c (c),(e)	spin component of object i = {1, 2}	'chi1x', 'chi1y', 'chi1z'	⋯	[−1, 1], (S2)
	along axis c = {x, y, z}	'chi2x', 'chi2y', 'chi2z'
Λi (d)	adimensional tidal deformability	'Lambda1', 'Lambda2'	⋯	[0, +∞)
	of object i = {1, 2}
e0	orbital eccentricity	'ecc'	⋯	[0, 1)
χs , χa (c)	symmetric and asymmetric spin	'chiS', 'chiA'	⋯	[−1, 1], [−1, 1]
	components; see Equation (1)
$\tilde{{\rm{\Lambda }}}$, $\delta \tilde{{\rm{\Lambda }}}{\,}^{({\rm{d}})}$	adimensional tidal deformability	'LambdaTilde',	⋯	[0, +∞),
	combinations; see Equation (2a)	'deltaLambda'		(−∞ , +∞)
θJN (e)	inclination angle w.r.t. total	'thetaJN'	rad	[0, π]
	angular momentum
χi (e)	spin magnitude of object i = {1, 2}	'chi1', 'chi2'	⋯	[0, 1]
θs,i (e)	spin tilt of object i = {1, 2}	'tilt1', 'tilt2'	rad	[0, π]
ϕJL (e)	azimuthal angle between orbital	'phiJL'	rad	[0, 2π]
	and total angular momentum
ϕ1,2 (e)	difference in azimuthal angle	'phi12'	rad	[0, 2π]
	between the spin vectors

Note. The first column reports the symbol used to denote a parameter, the second denotes a brief description of its physical meaning, the third depicts its name in GWFAST, the fourth shows the physical units of the parameter adopted in GWFAST, and the last column denotes its physical range. Parameters describing the same physical quantities, which thus have to be provided alternatively, are followed by a superscript in the first column, matching the one reported in the list in Section 2.1. S² in the χ_i,c stresses that the three components of a spin vector are not independent, but defined on a sphere, i.e., ${\chi }_{i,x}^{2}+{\chi }_{i,y}^{2}+{\chi }_{i,z}^{2}\leqslant 1$.

Download table as:  ASCIITypeset image

We use Fourier domain waveform models. In particular, at the time of writing, the code implements some selected waveform models in Python, and an interface to the LIGO Algorithm Library, LAL, which allows us to use all of the waveforms available in that library.

In particular, the following waveform models are directly available in GWFASTin a Python implementation:

• (i) a restricted post-Newtonian (PN) waveform model (Buonanno et al. 2009; Ajith 2011; Mishra et al. 2016), also with its tidal (Wade et al. 2014) and moderate eccentric (Moore et al. 2016) extensions. This is an inspiral-only waveform, but can still be used to describe signals coming from BNS mergers, whose major contribution to the S/N comes from the inspiral. There is no limitation in the parameters range, except for the eccentricity, which cannot exceed e₀ ∼ 0.1 for comparable mass systems; ⁴
• (ii) a full inspiral-merger-ringdown waveform model (Husa et al. 2016; Khan et al. 2016), tuned with numerical relativity (NR) simulations, which can efficiently be used to simulate signals coming from BBH mergers, with nonprecessing spins up to ∣χ_z ∣ ∼ 0.85 and mass ratios up to q = m₁/m₂ ∼ 18;
• (iii) tidal extension of the previous model (Dietrich et al. 2019), which can be used to accurately describe signals coming from BNS mergers. It includes spin terms up to higher order and a filter to terminate the waveform after merger. The validity has been assessed for masses ranging from 1–3 M_⊙, spins up to ∣χ_z ∣ ∼ 0.6, and tidal deformabilities up to Λ_i ≃ 5000;
• (iv) full inspiral-merger-ringdown waveform model (London et al. 2018; Kalaghatgi et al. 2020), which takes into account not only the quadrupole of the signal, but also the subdominant multipoles (l, m) = (2, 1), (3, 2), (3, 3), (4, 3), and (4, 4), which can be particularly relevant to better describe the signal coming from BBH systems. The calibration range is the same of the IMRPhenomD model;
• (v) full inspiral-merger-ringdown waveform model (Pannarale et al. 2015; Dietrich et al. 2019), which can describe the signal coming from the merger of an NS and a BH, with mass ratios up to q ∼ 100, also taking into account tidal effects and the impact of the possible tidal disruption of the NS.

These waveform models have been translated from their C implementation in the LIGO Algorithm Library (LAL) into a pure Python version. We carefully checked that our Python implementation accurately reproduces the original LAL waveforms, as can be seen on some example events in Figures 2 and 3. ⁵

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The implementation in Python has two advantages. First, it allows GWFAST to fully exploit the capabilities of this language to vectorize the computation on multiple events at a time, which would be impossible if we had to interact with a code written in C such as LAL. Second, it allows for the possibility of using AD (and in particular the library JAX) to compute derivatives; see Section 3.3.

It is also possible to use the waveform module separately from the rest of the code. For example, in order to generate the IMRPhenomD waveform amplitude and phase for a given set of events, it is sufficient to run the following sequence of commands:

Code example 2. Calculation of waveforms in GWFAST and WF4Py. As an illustration, we use the IMRPhenomD model.

This small piece of code shows that, as GWFAST's waveforms are written in pure Python and fully vectorized, our software does not have to rely on for loops over the events, as in a code interacting with C. Note that the order of the entries in the events dictionary is arbitrary.

All of our waveforms also include a routine to compute the time to coalescence as a function of frequency, needed to take into account Earth's rotation in the computation of the strain, which includes terms up to 3.5 PN order, ⁶ called tau_star, and a function to compute the cut frequency for the given waveform, so to properly build the frequency grid, called fcut, as seen in the above example. Waveform objects in GWFAST contain the attribute ParNums giving a dictionary of the form {'name_of_parameter':position}, with 'name_of_parameter' being a string with the parameter name as in Table 1 and position being an int corresponding to the position of the parameter in the Fisher matrix.

Apart from their implementation in GWFAST, which includes some features specific for JAX compatibility, we publicly release a pure numpy and ready-to-use version of the waveform models alone, WF4Py. ⁷ The syntax for using waveforms in this library is the same as in the example above. This module further implements the waveform model IMRPhenomXAS (Pratten et al. 2020), which is a full inspiral-merger-ringdown model tuned for the fundamental mode of BBH systems with aligned spins and mass ratios ranging from 1 to 1000, among the last to be released. ⁸

Finally, all waveform models available in LAL can be accessed in GWFAST through the wrapper class LAL_WF, which can be used as follows:

Code example 3. How to use LAL waveforms in GWFAST. As an illustration, we use the IMRPhenomXPHM model.

where the first entry has to be a string containing the name of the chosen waveform as in the LAL library ⁹ —'IMRPhenomXPHM' in this example—and the Booleans is_tidal, is_HigherModes, is_Precessing, and is_eccentric are used to specify whether the model includes tidal effects, the contribution of higher-order harmonics, precessing spins, or eccentricity (see also footnote 1), respectively.

A signal object can be initialized from the class GWSignal as follows:

Code example 4. Initialization of objects characterizing single detectors in GWFAST for a user-specified location, orientation, shape, and path to the PSD file.

where det_lat and det_long denote the latitude and longitude of the detector in degrees, respectively, xax is the angle between the bisector of the detector arms and the east in degrees, and detector_shape is its shape, which can be either 'L' (for L-shaped detectors) or 'T' (for triangular-shaped detectors). A triangular-shaped detector is defined as three colocated detectors with an opening angle of 60° in a closed-loop configuration.

Other options can be passed when initializing the GWSignal object, in particular: the useEarthMotion Boolean is used to turn on and off the computation of the effect of Earth's rotation; fmin and fmax can be used to set the minimum and maximum of the frequency grid (in hertz) and have default values fmin = 2. 0 and fmax = None (meaning that the grid extends up to the cut frequency of the waveform); DutyFactor can be used to set the duty factor of the detector, i.e., the fraction of time each detector is supposed to be operational, between 0 and 1 (default is None, meaning no duty cycle is considered). For triangular-shaped detectors, the duty factor refers to each of the three components of the triangle separately.

The entry psd_path is the path to the file containing its amplitude spectral density (ASD) or power spectral density (PSD); the flag is_ASD can be used to specify whether the given one is a PSD or ASD. GWFAST contains the following publicly available ASDs in the folder data/psds/:

• 1.  
  the sensitivities from Abbott et al. (2017; last update in January 2020), which can be found at https://dcc.ligo.org/LIGO-T1500293/public, in the folder unofficial_curves_all_dets;
• 2.  
  the representative sensitivities of the LIGO and Virgo detectors during the observing runs O1 and O2 (https://dcc.ligo.org/P1800374/public/), O3a (https://dcc.ligo.org/LIGO-P2000251/public), and O3b (estimated using PyCBC around the times reported in the caption of Figure 2 of Abbott et al. 2021a), in the folder LVC_O1O2O3;
• 3.  
  the sensitivities adopted in Abbott et al. (2020a) for the LIGO, Virgo, and KAGRA detectors during the O3, O4, and O5 runs (https://dcc.ligo.org/LIGO-T2000012/public), in the folder observing_scenarios_paper;
• 4.  
  the official ET–D sensitivity curve, from the document ET-0000A-18.txt (https://apps.et-gw.eu/tds/?content=3&r=14065); and
• 5.  
  the latest sensitivity curves for CE, used in Srivastava et al. (2022), for various detector configurations (https://dcc.cosmicexplorer.org/CE-T2000017/public), in the folder ce_strain.

GWFAST also contains some predefined detector configurations in the module globals. These are listed, together with their acronyms, in Table 2. The locations of the current detectors are taken from Gossan et al. (2022), while the CE sites are taken from Borhanian (2021) as illustrative examples. Predefined detector configurations can be easily imported from globals. In the following, we show how to initialize one signal corresponding to one CE at the Hanford site, with orientation γ = 0:

Code example 5. Initialization of objects characterizing single detectors in GWFAST using predefined detector configurations.

Any other user-defined configuration can easily be added as in Code example 4. With the object of type GWSignal initialized, the user can easily compute all of the quantities characterizing the signal. In particular, from the _PatternFunction function, it is possible to get the pattern functions of the detector; from GWAmplitudes, it is possible to compute the +' and ×' amplitudes of the signal at the detector (i.e., multiplied by the pattern functions and the spherical harmonics), while the full signal strain can be obtained through the function GWstrain.

Table 2. Summary of the Positions, Orientations, Angle between Arms, Shapes, and Acronyms of the Detectors Available in GWFAST

Detector	Latitude λ	Longitude φ	Orientation γ	Arms Aperture ζ	Shape	Name in GWFAST
LIGO Hanford, USA	465	−1194	171°	90°	'L'	'H1'
LIGO Livingston, USA	306	−908	2427	90°	'L'	'L1'
Virgo, Cascina, IT	436	105	1156	90°	'L'	'Virgo'
KAGRA, Hida, JP	364	1373	154	90°	'L'	'KAGRA'
LIGO India, Hingoli, IN	196	770	2874	90°	'L'	'LIGOI'
ET Sardinia, IT	405	94	0°	60°	'T'	'ETS'
ET Meuse–Rhine, EU	507	59	0°	60°	'T'	'ETMR'
CE1 Idaho, USA	438	−1128	−45°	90°	'L'	'CE1Id'
CE2 New Mexico, USA	332	−1065	−105°	90°	'L'	'CE2NM'
CE2 New South Wales, AU	−34°	145°	0°	90°	'L'	'CE2NSW'

Note. Using user-defined configurations is straightforward (see the text). The orientation γ denotes the angle between the bisector of the arms (the first arm in the case of a triangle) and East.

Download table as:  ASCIITypeset image

From more than one GWSignal object, one can define a network, which is composed by multiple detectors. The detectors composing the network have to be inserted into a dictionary, which can then be used to initialize an object from the class DetNet, characterizing the network:

Code example 6. Initialization of the network object in GWFAST from single detector objects.

From both the signal and network objects, it is then possible to compute the S/Ns and Fisher matrices for a set of events. The matched filter S/Ns in a single detector are computed using the definition

$$\begin{eqnarray}&&{{\rm{S}}/{\rm{N}}}_{i}^{2}=4{\int }_{{f}_{min}}^{{f}_{max}}\displaystyle \frac{| {\tilde{h}}_{(i)}(f){| }^{2}}{{S}_{n,i}(f)}df,\end{eqnarray} \tag{ 3 }$$

where ${\tilde{h}}_{(i)}(f)$ denotes the GW strain in Fourier domain at the ith detector, and S_n,i (f) is the noise spectral density of the ith detector. The network S/N is defined as the sum in quadrature of the single detectors' S/Ns:

$$\begin{eqnarray}&&{{\rm{S}}/{\rm{N}}}^{2}={\sum }_{i}{{\rm{S}}/{\rm{N}}}_{i}^{2}.\end{eqnarray} \tag{ 4 }$$

This can be obtained as simply as:

Code example 7. Computation of S/Ns (coded as "SNRs") in GWFAST both for a single detector and for a network.

The output of the methods above is a numpy array of the same length of the number of events.

The FIM elements for a single detector are computed from the definition

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{{ij}}=4{\int }_{0}^{\infty }{d}f\displaystyle \frac{{\partial }_{i}\tilde{h}{\partial }_{j}{\tilde{h}}^{* }}{{S}_{n}(f)},\end{eqnarray} \tag{ 5 }$$

where ∂_i denotes the derivative with respect to the parameter i. The FIM for the network is obtained by summing the individual Fisher matrices (which relies on the fact that different detectors are independent ¹¹ ). The FIM for a single detector or a network can be obtained with a single function call:

Code example 8. Computation of Fisher matrices in GWFAST both for a single detector and for a network.

The FIMs are returned by GWFAST in the form of a numpy array, treated as an array of matrices in the last dimension. For example, an array of FIMs for five BBH events with nine waveform parameters will have dimension (9, 9, 5). In the case of a network, it might be useful to store the S/Ns and Fisher matrices of the single detectors. This can be done passing to the functions SNR and FisherMatr the flag return_all = True. In this case, the output of both functions is a dictionary, with keys corresponding to the detectors (one key for each arm in the case of a triangular-shaped detector) and a key 'net' for the network S/Ns and FIMs.

The default parameters for which GWFAST computes the FIM, in the quasi-circular, nonprecessing, and nontidal case, are, in order, ¹² ${{ \mathcal M }}_{c}$ in units of M_⊙, η, d_L in units of gigaparsecs, θ, ϕ, ι, ψ, t_c in units of seconds, Φ_c , and χ_s , χ_a . In the case of precessing spins, the FIM is computed for the full set χ_1,z , χ_2,z , χ_1,x , χ_2,x , χ_1,y , χ_2,y in this order, and, in the BNS and NSBH case, the tidal parameters $\tilde{{\rm{\Lambda }}}$ and $\delta \tilde{{\rm{\Lambda }}}$ are also included. In the eccentric case, the parameter e₀ is also included, and appears in the Fisher matrix after both spins and tidal parameters. We chose to use the combinations (χ_s , χ_a ) instead of (χ_1,z , χ_2,z ) in the nonprecessing case so to have two orthogonal parameters, but the FIM can be as well computed in terms of the latter quantities passing the flag use_chi1chi2 = True to the FisherMatr function. The choice of the combination $(\tilde{{\rm{\Lambda }}},\delta \tilde{{\rm{\Lambda }}})$ in place of (Λ₁, Λ₂) is due to the fact that the parameter $\tilde{{\rm{\Lambda }}}$ is much better constrained than the two dimensionless tidal deformabilities separately, as it is the combination entering at 5 PN order in the inspiral signal. It is also possible to compute the FIM in terms of the combination (m₁, m₂), i.e., the two component redshifted masses, in units of M_⊙, instead of (M_c , η), by passing to the FisherMatr function the flag use_m1m2 = True. Finally, if the contribution of precessing spins is included, setting the flag use_prec_ang = True, instead of ι and χ_i,c , the FIM will be computed in terms of the parameters θ_JN , χ₁, χ₂, θ_s,1, θ_s,2, ϕ_JL , ϕ_1,2, which are more commonly used in the context of parameter estimation of GW events.

As an example, to access the values of the (d_L , d_L ) elements of the FIM for all of the events in the dictionary, the user just has to run:

Code example 9. How to access specific FIM elements in GWFAST.

Both of the classes GWSignal and DetNet also include a function to compute the optimal coordinates for a signal to be seen by the considered detectors (i.e., the location corresponding to the maximum S/N), as a function of the time of coalescence. This is obtained by maximizing the pattern functions, and can be accessed as:

Code example 10. How to compute the optimal location of a binary for a network of detectors in GWFAST at $\mathrm{GMST}=0\,{\rm{day}}$.

where the time can be provided both as a GMST or as a GPS time, setting the Boolean is_tGPS = True. The syntax to compute the optimal location is equivalent for an object of type GWSignal. ¹³

The computation of the derivatives of the signal with respect to its parameters is the key ingredient of a Fisher code. In its pure Python version, GWFAST is based on AD (Margossian 2018) as implemented in the library JAX (Bradbury et al. 2018). Differently from finite difference techniques, AD exploits the fact that each code function, however complex, is built up from elementary arithmetical operations and simple functions, whose derivatives are well known; thus, by applying the chain-rule repeatedly, it is possible to automatically compute derivatives of arbitrary order near machine precision, with a number of operations comparable to those of the original function. Having a properly written pure Python code, which is a fundamental requirement for this technique to work, it is possible to apply AD to get a fast and accurate evaluation of the GW strain derivatives, in a semianalytic way, despite the complexity of the function, and for multiple events at a time.

We tested the reliability of JAX by comparing the Fisher matrices obtained using the TaylorF2_RestrictedPN waveform model with GWFAST and an independent code written in Wolfram Mathematica, capable of computing derivatives analytically with respect to all parameters. The results for the relative differences of the diagonal elements of the Fisher matrices computed on a sample of 100 events are shown in Figure 4, from which it is possible to see the excellent agreement among the two codes. This test acts on three levels: (1) it proves the absence of bugs, given that the two codes were developed independently, (2) it shows the good behavior of JAX AD against an actual analytical computation, and (3) it verifies that integration is correctly performed, given the consistency of the results obtained with two different programming languages, having completely different integration routines.

Zoom In Zoom Out Reset image size

GWFAST also allows the user to compute analytically the derivatives of the signal with respect to many of the parameters, namely d_L , θ, ϕ, ι, ψ, t_c , and Φ_c to further speed up the calculation. ¹⁴ This can be done passing the flag computeAnalyticalDeriv = True to the function FisherMatr. We checked that, for these parameters, the analytical results and the result obtained by JAX agree at machine precision, i.e., 10⁻¹⁵. Finally, the FIM for a triangular detector is computed by using the fact that, for a closed configuration, the sum of the signals is identically zero for geometrical reasons (Freise et al. 2009), which further reduces the computational time by one-third. ¹⁵ When instead using waveforms coming from LAL, which are written in C, GWFAST will compute derivatives using the library numdifftools, ¹⁶ which relies on finite differences techniques. In this case, the computation is performed using the central differencing scheme with an adaptive computation of the step size. Both of these choices can be controlled though the arguments methodNDT and stepNDT, respectively. The finite difference computation can be used alternatively to AD also when exploiting Python waveforms, passing the flag computeDerivFinDiff = True to the function FisherMatr. Note that, also when using finite differences techniques, derivatives with respect to the parameters d_L , θ, ϕ, ι, ψ, t_c , and Φ_c can be performed analytically.

Tools for obtaining the covariance and analyzing the reliability of the inversion are contained in the module fisherTools. Note that all of the functions described here assume that the input FIM is an array of matrices in the last dimension, as described in Section 3.2 The conditioning of the FIM can be checked in GWFAST via the function CheckFisher, which returns the eigenvalues, eigenvectors, and condition number of the matrix. The inversion of the FIM to yield the covariance is done with the function CovMatr, as:

Code example 11. Computation of covariance matrices from Fisher matrices in GWFAST.

By default, each row and column is normalized to the square root of the diagonal of the FIM before inversion, so that the resulting matrix has adimensional entries with ones on the diagonal and the remaining elements in the interval [−1, 1] (Harms et al. 2022). ¹⁷ The inverse transformation is applied after inversion to yield the inverse of the original matrix. GWFAST also implements a variety of possibilities for the technique used to find the inverse. The Python library for precision arithmetic mpmath is used for the inversion. The inversion is not performed if the condition number is larger than a threshold that can be specified by the user via the argument condNumbMax. Its default value is 10⁵⁰ (so the code will try to invert every matrix irrespective of the conditioning). The available possibilities are listed below, and can be specified by the argument invMethodIn:

• 'inv': inverse computed by mpmath;
• 'cho': inverse computed by means of the Cholesky decomposition, i.e., the (Hermitian, positive-definite) FIM is expressed as a product of a lower triangular matrix and its conjugate transpose, and the latter is inverted. This is the default option in GWFAST; ¹⁸
• 'svd': the singular-value decomposition (SVD) of the FIM is used to invert the matrix. In this case, there is the additional option of truncating the smallest singular values to the minimum allowed numerical precision, which can help regularize poorly conditioned matrices. This can be required by setting the Boolean truncate = True. In this case, for each singular value s, if the ratio of its absolute value to the absolute value of the largest singular value, $\max {s}_{i}$, is smaller than a threshold λ, the singular value s is replaced with $\lambda \times \max ({s}_{i})$. The value of the threshold λ can be specified with the argument svals_thresh, which is set by default to 10¹⁵;
• 'svd_reg': the SVD of the FIM is used to invert the matrix, and eigenvalues smaller than the threshold specified by the argument svals_thresh are not included in the inversion. This ensures that the error on poorly constrained parameters is not propagated to the other ones (Harms et al. 2022). However, it might result in underestimating the uncertainty for parameters whose eigenvalues are excluded, and the effect should be carefully checked.
• 'lu': inversion is done by means of the lower–upper decomposition, i.e., the factorization of the FIM into the product of one lower triangular matrix and one upper triangular matrix. This can be a useful option since, as for the Cholesky decomposition, the inversion of a triangular matrix is easier than the one of a full matrix. Differently from the Cholesky decomposition, however, the original matrix does not have to be hermitian and positive-definite, which can make this method more stable against numerical noise for poorly conditioned matrices.

The error on the inversion is computed in GWFAST by the function compute_inversion_error with the definition $\epsilon ={| | {\rm{\Gamma }}\cdot {{\rm{\Gamma }}}^{-1}-{\mathbb{1}}| | }_{\max }={\max }_{{ij}}| {({\rm{\Gamma }}\cdot {{\rm{\Gamma }}}^{-1}-{\mathbb{1}})}_{{ij}}|$, where 1 denotes the identity matrix, Γ is the FIM, and Γ⁻¹ is its inverse as computed by the code:

Code example 12. Computation of the inversion errors in GWFAST.

Two other utilities to check the quality of the inversion are available in GWFAST in the module fisherTools. The function check_covariance computes the inversion error, and prints the difference between Γ · Γ⁻¹ and the identity on the diagonal, and the off-diagonal elements of Γ · Γ⁻¹ exceeding a given threshold specified with the argument tol. Second, the function perturb_Fisher adds random perturbations to the FIM to a specified decimal (given by the argument eps, whose default is 10⁻¹⁰), and checks if the inversion remains stable.

While the square root of the diagonal elements of the covariance matrix give the expected marginalized 1σ errors on the parameters, a useful metric for GW parameter estimation is the sky localization region at some given confidence level. This is computed by (Barack & Cutler 2004; Wen & Chen 2010)

$$\begin{eqnarray}&&\begin{array}{l}{\rm{\Delta }}{{\rm{\Omega }}}_{{X} \% }=-2\pi | \sin \theta | \sqrt{{\left({{\rm{\Gamma }}}^{-1}\right)}_{\theta \theta }{\left({{\rm{\Gamma }}}^{-1}\right)}_{\phi \phi }-{\left({{\rm{\Gamma }}}^{-1}\right)}_{\theta \phi }^{2}}\\ \quad \times \,\mathrm{ln}\left(1-{X}/100\right).\end{array}\end{eqnarray} \tag{ 6 }$$

The function compute_localization_region computes the sky localization region, in square degrees or steradians, according to the previous definition. The desired units can be specified through the units key, which can have values 'SqDeg' and 'Sterad', and the confidence level is specified by the optional argument perc_level (with default 90%). An example of usage is presented in Code example 13.

The Fisher approach allows us to treat straightforwardly some common situations encountered in parameter estimation, which we summarize here together with a description of their implementation in GWFAST. All functions described in the following belong to the module fisherTools:

• 1.  
  In order to fix some parameters to their fiducial values, one has to remove from the FIM (before inverting it) the corresponding rows and columns. This is done with the function fixParams, which takes the following arguments (in the order they are listed here): the original matrix, the dictionary specifying the position of each parameter in the FIM (accessible from the waveform object; see the end of Section 2.2 for an explanation), and a list of string with names of the parameters to be fixed, with the same names as in Table 1. The function returns the new matrix, and a dictionary of the same form of the input dictionary, with the keys corresponding to the fixed parameters removed, and the remaining rescaled;
• 2.  
  In order to add a Gaussian prior on some parameters, one has to add to the FIM a prior matrix P_ij corresponding to the inverse covariance of the prior. For the moment, GWFAST supports the addition of a diagonal prior matrix. This can be done with the function addPrior, which takes as input the original matrix, a list of values to be added on the diagonal of the Fisher (representing thus the inverse covariance of the prior on the corresponding parameter), the dictionary specifying the position of each parameter in the FIM, and the list of names of parameters on which the prior should be added;
• 3.  
  In order to marginalize over some parameters, one has to remove from the covariance matrix (after the inversion of the FIM) the corresponding rows and columns. This can be done again with the function fixParams described in the first point.

Code example 13. Example of manipulations of the Fisher matrix: fix the spins to their fiducial values, add a Gaussian prior on the angles with standard deviation 2π, compute the corresponding covariance, compute the forecasted 90% localization region in square degrees.