# Neural Astrophysical Wind Models

Dustin Nguyen

###### Abstract

The bulk kinematics and thermodynamics of hot supernovae-driven galactic winds
is critically dependent on both the amount of swept up cool clouds and non-
spherical collimated flow geometry. However, accurately parameterizing these
physics is difficult because their functional forms are often unknown, and
because the coupled non-linear flow equations contain singularities. We show
that deep neural networks embedded as individual terms in the governing
coupled ordinary differential equations (ODEs) can robustly discover both of
these physics, without any prior knowledge of the true function structure, as
a supervised learning task. We optimize a loss function based on the Mach
number, rather than the explicitly solved-for 3 conserved variables, and apply
a penalty term towards near-diverging solutions. The same neural network
architecture is used for learning both the hidden mass-loading and surface
area expansion rates. This work further highlights the feasibility of neural
ODEs as a promising discovery tool with mechanistic interpretability for non-
linear inverse problems.

Machine Learning, ICML

##  1 Introduction

Mathematics and the physical laws that have followed suit are unreasonably
effective in creating models with predictive power for the natural sciences
(Wigner, 1960). In physics, understanding implies that the governing
principles are written in the form of explicit differential equations of
motion (Meiss, 2007). However, when modelling real world systems, we do not
always know what the true underlying physics is.

In the context of explaining observations of nearby galactic superwinds,
recent works have suggested the inclusion of additional physics can help
overcome the so-called “cloud-crushing problem” (Klein et al., 1994), which
challenges the existence of observed high velocity cool clouds. These efforts
include incorporation of efficient multi-phase mixing (Gronke & Oh, 2018,
2020), cosmic ray acceleration (Quataert et al., 2022b, a), among many other
mechanisms (Faucher-Giguere & Oh, 2023). While incorporation of new physics
into semi-analytic models and 3D simulations are a vibrant and active field of
research (e.g., Schneider & Robertson, 2015; Kim et al., 2020; Pandya et al.,
2022; Wibking & Krumholz, 2022; Smith et al., 2023; Nguyen et al., 2023; Tan &
Fielding, 2023, and many others), the development of optimization methods to
systematically discover and characterize the structure of physics from
observational data has not received similar attention. Namely, descriptions of
superwinds from nearby starburst prototype M82 have relied on parameter
estimation of classic theoretical models (Chevalier & Clegg, 1985; Strickland
& Heckman, 2009; Lopez et al., 2020; Nguyen & Thompson, 2021; Yuan et al.,
2023). Methods for developing new data-driven models have not yet been
explored for such systems.

Deep neural networks (LeCun et al., 2015; Goodfellow et al., 2016) can be used
as a universal approximator for unknown functions and operators (Hornik et
al., 1990; Lu et al., 2021). Recently, deep neural networks have been used as
numerical discretization schemes of ODEs/PDEs (Chen et al., 2019) based on the
adjoint sensitivity method. This idea has been expanded to a hybrid model
termed the Universal Differential Equation (UDE, Rackauckas et al., 2020) that
allows one to freely augment ODEs/PDEs with universal approximators such as
neural networks. Embedded neural networks are guaranteed to obey conservation
laws, as they can only influence the dynamics of the system term by term. In
this paper, neural ODEs are ODEs with a neural network as a part of the
equations and are solved using standard numerical methods (unlike Physics-
Informed Neural Networks, where the numerical solver is a “black box”). Neural
networks embedded within ODEs/PDEs are a new and active field of research and
have been used to predict COVID-19 pandemic waves (Kuwahara & Bauch, 2023),
and various other applications within biology, chemistry, mathematics, and
physics (Vortmeyer-Kley et al., 2021; Gelbrecht et al., 2021; Keith et al.,
2021; Fronk & Petzold, 2023; Stepaniants et al., 2023; Yin et al., 2023;
Santana et al., 2023).

We present the first analysis of neural coupled ODEs, within the UDE
framework, that contain singularities. To illustrate generality, we consider
two experiments where the embedded neural network learns two different types
of physics: 1. mass-loading of swept up clouds, and, 2. the surface area
expansion rate attributed to flow geometry. We define a custom loss function
based solely on the Mach number, rather than the 3 individual conserved
variables, and adopt linear scaling towards early integration steps.
Additionally, this loss function penalizes solutions that approach the sonic
point, which is infamously known to lead to numerical instability (Lamers &
Cassinelli, 1999). We will show that regression on such a loss function is
sufficient to characterize the underlying $v$, $\rho$, and $P$ profiles and
discover the hidden physics embedded within the training data.

##  2 Methods

The hydrodynamic equations for a non-radiative highly supersonic steady-state
hot flow moving in the $x$ direction are (Cowie et al., 1981; Nguyen &
Thompson, 2021)

$\displaystyle\frac{1}{A}\frac{d}{dx}(\rho vA)=\dot{\mu},$ (1) $\displaystyle
v\frac{dv}{dx}=-\frac{1}{\rho}\frac{dP}{dx}-\frac{\dot{\mu}v}{\rho},$ (2)
$\displaystyle
v\frac{d\epsilon}{dx}-\frac{vP}{\rho^{2}}\frac{d\rho}{dx}=\frac{\dot{\mu}}{\rho}\bigg{(}\frac{v^{2}}{2}-\epsilon-\frac{P}{\rho}\bigg{)},$
(3)

where $v$, $\rho$, $P$, $\epsilon$, $\dot{\mu}$, and $A(x)$ are the bulk
velocity, density, pressure, specific internal energy, volumetric mass-loading
rate, and flow area, respectively. By substituting the equations into each
other, the derivatives will each contain the singularity:
$(\mathcal{M}^{2}-1)^{-1}$ (i.e., is a pole at the sonic point where
$\mathcal{M}=1$. See Lamers & Cassinelli, 1999). We consider two experiments
of training a neural network to learn: 1. the volumetric mass-loading rate,
$\dot{\mu}$, and 2. the surface area expansion rate $d\,\ln A/dx$. In each,
the training dataset is calculated with the analytic functional form of either
$\dot{\mu}(x)$ or $d\,\ln A/dx$. After generating the training dataset, we
forget the true physical variable and replace it with neural network $\Phi(x)$
that has a total of 5 layers (See Fig. 1 for details).

Figure 1: Architecture of neural network $\Phi(x)$. The input is a single
position $x$, which then gets passed through a total of 5 Dense layers where
the three hidden layers contain skip connections and 64 nodes. Each layer uses
swish as the activation function, with the output layer using relu to ensure
positivity.

The dimensions of each simulation is $0.3\,\mathrm{kpc}\leq x\leq
2\mathrm{kpc}$ and number of steps $n_{x}=500$. We have assumed the flow to
have already left the starburst volume $R=0.3\,$kpc (Chevalier & Clegg, 1985).
The initial conditions are then $v_{0}=1500\,\mathrm{km\,s^{-1}}$,
$n_{0}=10\,\mathrm{cm^{-3}}$, and $T_{0}=10^{7}\,\mathrm{K}$, which
corresponds to a highly supersonic $(\mathcal{M}_{0}\sim 3)$ non-radiative
wind that has asymptotic velocities roughly within the allowed range for
starburst galaxy M82 (Strickland & Heckman, 2009).

We calculate the loss function using the Mach number $\mathcal{M}=v/(\gamma
P/\rho)^{1/2}$. We find this method to outperform summing Mean-Squared-Error
(MSE) of the individual conserved variables $v$, $\rho$, and $P$, even when
using min-max normalization. The loss function is

$\mathcal{L}=\sum_{i}^{n_{x}}\bigg{[}\mathcal{W}_{i}\times\bigg{(}\mathcal{M}_{i}-\hat{\mathcal{M}}_{i}\bigg{)}^{2}\bigg{]}+\zeta_{\mathrm{penalty}}(\mathcal{M})$
(4)

The first term is a weighted MSE, where the weights $\mathcal{W}_{i}$ linearly
scale the MSE as a function of $n_{x}$, where $\mathcal{W}_{0}=1$ and
$\mathcal{W}_{n_{x}}=\kappa<1$. We do not include division by $n_{x}$, as it
does not impact training. We vary $\kappa$ during training for each
optimization algorithm (details below). This scaling increases sensitivity to
early solutions, which is important for non-linear problems. The latter term
$\zeta_{\mathrm{penalty}}(\mathcal{M})$ penalizes solutions with Mach numbers
close to $1$, as to prevent the neural network from sampling diverging
solutions ($\mathcal{M}=1$) and is a function of the total loss
$\mathcal{L}_{j}$ and $\mathcal{M}_{i}$ as:

$\displaystyle\zeta_{\mathrm{penalty},j}=\mathcal{L}_{j}\times\omega\sum_{i}^{nx}\bigg{[}1-\bigg{(}1-\hat{\mathcal{M}}_{i}\bigg{)}^{2}\bigg{]}$
(5)
$\displaystyle\quad\quad\quad\quad\quad\quad(\mathcal{M}_{\mathrm{min}}\leq\hat{\mathcal{M}}_{i}\leq\mathcal{M}_{\mathrm{max}}),$
(6)

where $j$ is the training iteration, $\omega$ is the penalty weight,
$\mathcal{M}_{\mathrm{min}}=1$, and $\mathcal{M}_{\mathrm{max}}=1.5$ sets the
threshold for penalization. We summarize the training process in the pseudo-
code below.

Algorithm 1 Training Psuedo-Code

Input: Initial Condition Array $(v_{0},\rho_{0},P_{0})$ size $1\times 3$

1\.  Calculate true solution using $\dot{\mu}(x)$ or $A(x)$ in Eqs. 1,2,3. The
solution is then a matrix $y$ of size $n_{x}\times 3$. Use $y$ to calculate
$\mathcal{M}$ which is size $n_{x}\times 1$.

2\.  Calculate prediction using neural network $\Phi(x)_{p_{0}}$ within neural
equations. $p_{0}$ is initial parameters array for constructing
$\Phi(x)_{p_{0}}$. The solution is a matrix $\hat{y}_{0}$ of size $n_{x}\times
3$. Use $\hat{y}_{0}$ to calculate $\hat{\mathcal{M}_{0}}$ which is size
$n_{x}\times 1$.

for $j=0$ to $n_{\mathrm{iters}}$ do

Reconstruct neural net $\Phi(x)_{p_{j}}$ using parameters $p_{j}$.

Calculate $\hat{y}_{j}$, then $\hat{\mathcal{M}}_{j}$.

Calculate loss $\mathcal{L}_{j}$ using MSE between $\mathcal{M}$ and
$\hat{\mathcal{M}}_{j}$ (batched over all positions $n_{x}$)

if
$\mathcal{M}_{\mathrm{min}}\leq\hat{\mathcal{M}}_{j}\leq\mathcal{M}_{\mathrm{max}}$
(batch) then

$\mathcal{L}_{j}+=\mathcal{L}_{j}*\mathrm{Penalty\ Term}$

end if

Optimize for new set of parameters $p_{j}$ that minimize $\mathcal{L}_{j}$
using ADAM or BFGS optimization algorithms.

end for if $j=n_{\mathrm{iters}}$ or gradient tolerance met

We solve the equations (both classical and neural) using standard numerical
methods involving 5th order RK4 (Tsitouras, 2011) with adaptive step size
control (Rackauckas & Nie, 2017). We calculate gradients using forward-mode
automatic differentiation (Revels et al., 2016). All of the work is done using
packages contained within Julia’s Scientific Machine Learning (SciML)
ecosystem (Rackauckas et al., 2020).

##  3 Results

Learning Mass-loading: We test the neural network’s ability to learn mass-
loading $\dot{\mu}(x)=\dot{\mu}(\lambda/x)^{\Delta}$. This represents intense
cloud entrainment immediately after the wind leaves the host galaxy
($\lambda=0.3\,\mathrm{kpc}$, $\Delta=4$, and
$\dot{\mu}_{0}=500\,\mathrm{M_{\odot}\,yr^{-1}\,kpc^{-3}}$). Here we take the
flow geometry to be spherical (i.e., $A(x)\propto x^{2}$). After calculating
the training (true) solution with this function, we “forget” $\dot{\mu}(x)$
and replace it with a neural network $\Phi(x)$ (See Fig. 1) such that
Equations 1, 2, and 3 are re-written as neural wind equations:

$\displaystyle\frac{1}{x^{2}}\frac{d}{dx}(\rho vx^{2})=\Phi,$ (7)
$\displaystyle v\frac{dv}{dx}=-\frac{1}{\rho}\frac{dP}{dx}-\frac{\Phi
v}{\rho},$ (8) $\displaystyle
v\frac{d\epsilon}{dx}-\frac{vP}{\rho^{2}}\frac{d\rho}{dx}=\frac{\Phi}{\rho}\bigg{(}\frac{v^{2}}{2}-\epsilon-\frac{P}{\rho}\bigg{)}.$
(9)

We initialize the neural network with an output of approximately 0 (using
Flux.jl’s default Glorot initialization). This implies that the first
prediction solution is identical to an adiabatic spherical wind with no mass-
loading (i.e., $\dot{\mu}\simeq 0$). We train for $n_{\mathrm{iters}}=300$
iterations using ADAM with a learning rate of $0.01$ and exponential decay
rate of $0.98$, followed by additional iterations using BFGS. We use
$\kappa=0.1$ for scaling solutions. We use penalty term $\omega=0.1$ and
$\omega=0.75$ for ADAM and BFGS, respectively. Optimization with BFGS
converges to the set gradient tolerance of $10^{-9}$ after 304 iterations.

In the left three panels of Figure 2, we plot three different solutions for
velocity, number density ($n=\rho/\mu m_{p})$, and temperature ($T=P/nk_{B})$
radial profile predictions from: the true solution (solid black), the
untrained model (orange dashed), the ADAM trained model (purple dashed), and
the ADAM \+ BFGS trained (blue dashed). Although we explicitly train on only
the Mach number, the three conserved variables ($v$, $\rho$, and $P$) end up
getting well-matched (blue dashed lines). In the right most panel of Figure 2,
we plot various volumetric mass-loading rates calculated with the true
function $\dot{\mu}=\dot{\mu}_{0}(\lambda/x)^{\Delta}$ (black solid) and the
output of the neural networks (colored dashed). In Figure 3, we plot the
learning curve. The bumps in the ADAM portion of the learning curve (first 300
iterations), demonstrate the optimizer sampling solutions within the penalty
threshold, and then updating the mass-loading rate, encoded by $\Phi(x)$, for
new predictions that move away from such penalties. Training decreases the
normalized loss by nearly 7 orders of magnitude. The predictions of the neural
ODEs agree with the true solutions, and the output of the neural network
indicate that it has learned the true mass-loading function.

Figure 2: Density, temperature, and velocity profiles for solutions (left 3
panels) calculated with the mass-loading (right panel) given by the true
function (black solid) and the output of neural networks (colored dashed).
After BFGS and ADAM training, the three kinematic and thermodynamic profiles
are well-matched (blue dashed on black lines), and the underlying mass-loading
function was learned. Figure 3: Learning curve for training neural network
$\Phi(x)$ to learn underlying volumetric mass-loading rate $\dot{\mu}$ from
the solved wind solutions. The bumps in the loss function during ADAM training
epochs are attributed to penalties for near-diverging solutions.

Learning Surface Area Expansion: We now test the neural network’s ability to
learn properties of the flow geometry $A(x)=A_{0}(1+(x/\eta)^{2})$, which
describes a flared-cylinder flow tube that has previously been used for solar
coronal hole flux tubes (Kopp & Holzer, 1976) and Milky-Way disk winds
(Everett et al., 2008). Qualitatively, this function implies a roughly
constant flow area up until $\eta$ when the flow begins to undergo spherical
expansion. For illustration, Equation 1 can be expanded into differential form
as: $d\ln\rho/dx+d\ln v/dx+d\ln A/dx=\dot{\mu}/\rho v$. This means that the
effects of flow geometry are governed specifically by the surface area
expansion rate $d\,\ln A/dx$, and not the magnitude of $A(x)$ (e.g., flows
expanding into $4\pi x^{2}$ or $2\pi x^{2}$ behave the same). Setting
$\dot{\mu}=0$ we generate the training dataset using
$A(x)=A_{0}(1+(x/\eta)^{2})$ so that $d\ln A/dx=2x/(\eta^{2}+x^{2})$, where
$A_{0}=0.25\,\mathrm{kpc}$ and $\eta=1.0\,$kpc. We now “forget” $d\ln A/dx$,
replacing it with neural network $\Phi(x)$. Again, we initialize the neural
network with output of approximately 0. This implies that the first prediction
solution is identical to an adiabatic planar wind (i.e.,
$A(x)\simeq\mathrm{constant}$). We train for 300 iterations using ADAM with a
learning rate of $0.01$ and exponential decay rate of $0.98$, followed by
additional iterations using BFGS. We take $\kappa=0.1$ and $\omega=0.1$.
Optimization with BFGS converges to the set gradient tolerance after 47
iterations.

Figure 4: Left: Mach number solutions calculated with the true area function
(black solid), predictions from the neural ODEs (colored dashed lines). Right:
The true surface area expansion rate for a flared cylinder $d\log A/dx$ (black
solid) and the output from neural network $\Phi(x)$ for various degrees of
training.

In the left panel of Figure 4, we plot the Mach number solution calculated
with the true area function (black line) and from predictions of the neural
ODEs with varied degrees of training (colored dashed lines). We see that the
best fitting line is, again, the model that underwent both ADAM and BFGS
optimization. In the right panel we plot the surface area expansion rate
calculated for the true surface area function (black), and the direct output
of neural networks (colored dashed), finding agreement between the two (blue
dashed and solid black).

##  4 Conclusion and Future Work

In this work we show that deep neural networks embedded as individual terms
within astrophysical wind equations can discover the true functional form of
both mass-loading and flow geometry expansion, without any prior knowledge.
Additionally, to our knowledge at the time of writing, this is the first use
of the UDE framework for coupled ODEs that contain singularities (i.e., poles
at $\mathcal{M}=1$). To avoid diverging solutions, we introduced a penalty
term that activates when the predicted solution falls within a threshold. This
method is shown to be effective, and guides the optimization away from
diverging solutions (Fig. 3). We require BFGS optimization after ADAM to
converge to the functional form of the true solution. BFGS is a quasi-Newton
method that uses second-order derivative information that outperforms ADAM
here at low values of normalized loss. We note that using ADAM at the
beginning of training is still beneficial, as it is faster and can overcome
local minima better. We demonstrate the same neural network architecture can
be used to learn two very different types of physics. We show that the outputs
of the neural network are interpretable (Figs. 2 and 4). We find that training
on the Mach number $\mathcal{M}$ outperformed training on the conserved
variables $v$, $\rho$, and $P$. Beyond winds, we expect neural ODEs to be
useful in discovering physics in other non-linear systems.

In future work we would like to explore the predictive power of the trained
neural networks, beyond the training dataset. Here, we specifically studied
the ability of a neural network embedded within coupled ODEs to accurately
capture physics during supervised training, and showcase the interpretability
of the result. One possible direction is to use the SINDy (Brunton et al.,
2016) algorithm to extract symbolic representations of the trained neural
networks which has been shown to be useful in creating hybrid models capable
of making accurate long-term predictions for partially observed non-linear
systems (Rackauckas et al., 2020; Santana et al., 2023; Fronk & Petzold,
2023). However, this approach requires the creation of a basis function
library, which may not be generally applicable, especially when complicated
non-polynomial functions are involved. If we continue to use the direct output
of the trained neural network, careful considerations will need to be made
towards: 1. more advanced neural network architectures that better capture
long-term dependencies such as LSTMs (as opposed to the simple forward-feed
neural network used here), and 2. overfitting. Additionally, future work will
investigate how noisy data may affect learning outcomes. We note that even in
cases where predictive power beyond the training dataset is limited, the use
of neural ODEs as a discovery tool for a contained dataset is still useful, as
we demonstrate here.

## Acknowledgements

DN thanks Todd Thompson for useful discussions. DN also thanks Drummond
Fielding, Peng Oh, Brent Tan, Viraj Pandya, Ryan Farber, Evan Schneider, and
Max Gronke for insightful conversations at the “Modelling of Multi-phase Gas
in Astrophysical Media 2023” conference. DN thanks Chris Rackauckas for useful
pointers on using Julia. DN acknowledges funding from NASA 21-ASTRO21-0174.

## References

  * Brunton et al. (2016) Brunton, S. L., Proctor, J. L., and Kutz, J. N.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  _Proceedings of the National Academy of Sciences_ , 113(15):3932–3937, 2016.  doi: 10.1073/pnas.1517384113.
  * Chen et al. (2019) Chen, R. T. Q., Rubanova, Y., Bettencourt, J., and Duvenaud, D.  Neural ordinary differential equations, 2019.
  * Chevalier & Clegg (1985) Chevalier, R. A. and Clegg, A. W.  Wind from a starburst galaxy nucleus.  _Nature_ , 317(6032):44–45, Sep 1985.  doi: 10.1038/317044a0.
  * Cowie et al. (1981) Cowie, L. L., McKee, C. F., and Ostriker, J. P.  Supernova remnant revolution in an inhomogeneous medium. I - Numerical models.  _The Astrophysical Journal_ , 247:908–924, Aug 1981.  doi: 10.1086/159100.
  * Everett et al. (2008) Everett, J. E., Zweibel, E. G., Benjamin, R. A., McCammon, D., Rocks, L., and Gallagher, III, J. S.  The Milky Way’s Kiloparsec-Scale Wind: A Hybrid Cosmic-Ray and Thermally Driven Outflow.  _The Astrophysical Journal_ , 674:258–270, February 2008.  doi: 10.1086/524766.
  * Faucher-Giguere & Oh (2023) Faucher-Giguere, C.-A. and Oh, S. P.  Key Physical Processes in the Circumgalactic Medium.  _arXiv e-prints_ , art. arXiv:2301.10253, January 2023.  doi: 10.48550/arXiv.2301.10253.
  * Fronk & Petzold (2023) Fronk, C. and Petzold, L.  Interpretable polynomial neural ordinary differential equations.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ , 33(4):043101, 2023.  doi: 10.1063/5.0130803.
  * Gelbrecht et al. (2021) Gelbrecht, M., Boers, N., and Kurths, J.  Neural partial differential equations for chaotic systems.  _New Journal of Physics_ , 23(4):043005, April 2021.  doi: 10.1088/1367-2630/abeb90.
  * Goodfellow et al. (2016) Goodfellow, I., Bengio, Y., and Courville, A.  _Deep Learning_.  The MIT Press, Cambridge, MA, 2016.  ISBN 9780262035613.
  * Gronke & Oh (2018) Gronke, M. and Oh, S. P.  The growth and entrainment of cold gas in a hot wind.  _Monthly Notices of the Royal Astronomical Society_ , 480(1):L111–L115, October 2018.  doi: 10.1093/mnrasl/sly131.
  * Gronke & Oh (2020) Gronke, M. and Oh, S. P.  How cold gas continuously entrains mass and momentum from a hot wind.  _Monthly Notices of the Royal Astronomical Society_ , 492(2):1970–1990, February 2020.  doi: 10.1093/mnras/stz3332.
  * Hornik et al. (1990) Hornik, K., Stinchcombe, M., and White, H.  Universal approximation of an unknown mapping and its derivatives using multilayer feedforward networks.  _Neural Networks_ , 3(5):551–560, 1990.  ISSN 0893-6080.  doi: https://doi.org/10.1016/0893-6080(90)90005-6.  URL https://www.sciencedirect.com/science/article/pii/0893608090900056.
  * Keith et al. (2021) Keith, B., Khadse, A., and Field, S. E.  Learning orbital dynamics of binary black hole systems from gravitational wave measurements.  _Physical Review Research_ , 3(4):043101, November 2021.  doi: 10.1103/PhysRevResearch.3.043101.
  * Kim et al. (2020) Kim, C.-G., Ostriker, E. C., Fielding, D. B., Smith, M. C., Bryan, G. L., Somerville, R. S., Forbes, J. C., Genel, S., and Hernquist, L.  A Framework for Multiphase Galactic Wind Launching Using TIGRESS.  _The Astrophysical Journal Letters_ , 903(2):L34, November 2020.  doi: 10.3847/2041-8213/abc252.
  * Klein et al. (1994) Klein, R. I., McKee, C. F., and Colella, P.  On the Hydrodynamic Interaction of Shock Waves with Interstellar Clouds. I. Nonradiative Shocks in Small Clouds.  _The Astrophysical Journal_ , 420:213, January 1994.  doi: 10.1086/173554.
  * Kopp & Holzer (1976) Kopp, R. A. and Holzer, T. E.  Dynamics of coronal hole regions. I. Steady polytropic flows with multiple critical points.  _solphys_ , 49(1):43–56, July 1976.  doi: 10.1007/BF00221484.
  * Kuwahara & Bauch (2023) Kuwahara, B. S. and Bauch, C. T.  Predicting covid-19 pandemic waves with biologically and behaviorally informed universal differential equations.  In _medRxiv_ , 2023.
  * Lamers & Cassinelli (1999) Lamers, H. J. G. L. M. and Cassinelli, J. P.  _Introduction to Stellar Winds_.  1999\.
  * LeCun et al. (2015) LeCun, Y., Bengio, Y., and Hinton, G.  _Deep Learning_.  Nature, New York, NY, 2015.  doi: 10.1038/nature14539.  URL https://www.nature.com/articles/nature14539.
  * Lopez et al. (2020) Lopez, L. A., Mathur, S., Nguyen, D. D., Thompson, T. A., and Olivier, G. M.  Temperature and Metallicity Gradients in the Hot Gas Outflows of M82.  _The Astrophysical Journal_ , 904(2):152, December 2020.  doi: 10.3847/1538-4357/abc010.
  * Lu et al. (2021) Lu, Q., Mao, Z., Zuo, W., and Dong, B.  Learning nonlinear operators via deeponet based on the universal approximation theorem of operators.  _Nature Machine Intelligence_ , 3:270–278, 2021.  doi: 10.1038/s42256-021-00302-5.  URL https://www.nature.com/articles/s42256-021-00302-5.
  * Meiss (2007) Meiss, J. D.  _Differential Dynamical Systems_.  Society for Industrial and Applied Mathematics, 2007.  doi: 10.1137/1.9780898718232.  URL https://epubs.siam.org/doi/abs/10.1137/1.9780898718232.
  * Nguyen & Thompson (2021) Nguyen, D. D. and Thompson, T. A.  Mass-loading and non-spherical divergence in hot galactic winds: implications for X-ray observations.  _Monthly Notices of the Royal Astronomical Society_ , October 2021.  doi: 10.1093/mnras/stab2910.
  * Nguyen et al. (2023) Nguyen, D. D., Thompson, T. A., Schneider, E. E., Lopez, S., and Lopez, L. A.  Dynamics of hot galactic winds launched from spherically-stratified starburst cores.  _Monthly Notices of the Royal Astronomical Society_ , 518(1):L87–L91, January 2023.  doi: 10.1093/mnrasl/slac141.
  * Pandya et al. (2022) Pandya, V., Fielding, D. B., Bryan, G. L., Carr, C., Somerville, R. S., Stern, J., Faucher-Giguere, C.-A., Hafen, Z., and Angles-Alcazar, D.  A unified model for the co-evolution of galaxies and their circumgalactic medium: the relative roles of turbulence and atomic cooling physics.  _arXiv e-prints_ , art. arXiv:2211.09755, November 2022.  doi: 10.48550/arXiv.2211.09755.
  * Quataert et al. (2022a) Quataert, E., Jiang, Y.-F., and Thompson, T. A.  The physics of galactic winds driven by cosmic rays - II. Isothermal streaming solutions.  _Monthly Notices of the Royal Astronomical Society_ , 510(1):920–945, February 2022a.  doi: 10.1093/mnras/stab3274.
  * Quataert et al. (2022b) Quataert, E., Thompson, T. A., and Jiang, Y.-F.  The physics of galactic winds driven by cosmic rays I: Diffusion.  _Monthly Notices of the Royal Astronomical Society_ , 510(1):1184–1203, February 2022b.  doi: 10.1093/mnras/stab3273.
  * Rackauckas & Nie (2017) Rackauckas, C. and Nie, Q.  Differentialequations.jl–a performant and feature-rich ecosystem for solving differential equations in julia.  _Journal of Open Research Software_ , 5(1):15, 2017.
  * Rackauckas et al. (2020) Rackauckas, C., Ma, Y., Martensen, J., Warner, C., Zubov, K., Supekar, R., Skinner, D., and Ramadhan, A.  Universal differential equations for scientific machine learning.  _arXiv preprint arXiv:2001.04385_ , 2020.
  * Revels et al. (2016) Revels, J., Lubin, M., and Papamarkou, T.  Forward-mode automatic differentiation in Julia.  _arXiv:1607.07892 [cs.MS]_ , 2016.  URL https://arxiv.org/abs/1607.07892.
  * Santana et al. (2023) Santana, V. V., Costa, E., Rebello, C. M., Mafalda Ribeiro, A., Rackauckas, C., and Nogueira, I. B. R.  Efficient hybrid modeling and sorption model discovery for non-linear advection-diffusion-sorption systems: A systematic scientific machine learning approach.  _arXiv e-prints_ , art. arXiv:2303.13555, March 2023.  doi: 10.48550/arXiv.2303.13555.
  * Schneider & Robertson (2015) Schneider, E. E. and Robertson, B. E.  CHOLLA: A New Massively Parallel Hydrodynamics Code for Astrophysical Simulation.  _The Astrophysical Journal Supplement Series_ , 217(2):24, April 2015.  doi: 10.1088/0067-0049/217/2/24.
  * Smith et al. (2023) Smith, M. C., Fielding, D. B., Bryan, G. L., Kim, C.-G., Ostriker, E. C., Somerville, R. S., Stern, J., Su, K.-Y., Weinberger, R., Hu, C.-Y., Forbes, J. C., Hernquist, L., Burkhart, B., and Li, Y.  Arkenstone I: A Novel Method for Robustly Capturing High Specific Energy Outflows In Cosmological Simulations.  _arXiv e-prints_ , art. arXiv:2301.07116, January 2023.  doi: 10.48550/arXiv.2301.07116.
  * Stepaniants et al. (2023) Stepaniants, G., Hastewell, A. D., Skinner, D. J., Totz, J. F., and Dunkel, J.  Discovering dynamics and parameters of nonlinear oscillatory and chaotic systems from partial observations.  _arXiv e-prints_ , art. arXiv:2304.04818, April 2023.  doi: 10.48550/arXiv.2304.04818.
  * Strickland & Heckman (2009) Strickland, D. K. and Heckman, T. M.  Supernova Feedback Efficiency and Mass Loading in the Starburst and Galactic Superwind Exemplar M82.  _The Astrophysical Journal_ , 697:2030–2056, June 2009.  doi: 10.1088/0004-637X/697/2/2030.
  * Tan & Fielding (2023) Tan, B. and Fielding, D. B.  Cloud Atlas: Navigating the Multiphase Landscape of Tempestuous Galactic Winds.  _arXiv e-prints_ , art. arXiv:2305.14424, May 2023.  doi: 10.48550/arXiv.2305.14424.
  * Tsitouras (2011) Tsitouras, C.  Runge–kutta pairs of order 5(4) satisfying only the first column simplifying assumption.  _Computers & Mathematics with Applications_, 62(2):770–775, 2011.  ISSN 0898-1221.  doi: https://doi.org/10.1016/j.camwa.2011.06.002.  URL https://www.sciencedirect.com/science/article/pii/S0898122111004706.
  * Vortmeyer-Kley et al. (2021) Vortmeyer-Kley, R., Nieters, P., and Pipa, G.  A trajectory-based loss function to learn missing terms in bifurcating dynamical systems.  _Scientific Reports_ , 11:21181, Oct 2021.  doi: 10.1038/s41598-021-99609-x.  URL https://www.nature.com/articles/s41598-021-99609-x.
  * Wibking & Krumholz (2022) Wibking, B. D. and Krumholz, M. R.  QUOKKA: a code for two-moment AMR radiation hydrodynamics on GPUs.  _Monthly Notices of the Royal Astronomical Society_ , 512(1):1430–1449, May 2022.  doi: 10.1093/mnras/stac439.
  * Wigner (1960) Wigner, E. P.  The unreasonable effectiveness of mathematics in the natural sciences.  _Communications on Pure and Applied Mathematics_ , 13(1):1–14, 1960.
  * Yin et al. (2023) Yin, S., Wu, J., and Song, P.  Optimal control by deep learning techniques and its applications on epidemic models.  _Journal of Mathematical Biology_ , 87:2121–2148, Apr 2023\.  doi: 10.1007/s00285-023-01873-0.  URL https://link.springer.com/article/10.1007/s00285-023-01873-0.
  * Yuan et al. (2023) Yuan, Y., Krumholz, M. R., and Martin, C. L.  The observable properties of cool winds from galaxies, AGN, and star clusters - II. 3D models for the multiphase wind of M82.  _Monthly Notices of the Royal Astronomical Society_ , 518(3):4084–4105, January 2023.  doi: 10.1093/mnras/stac3241.


# Adapting the PyCBC pipeline to find and infer the properties of
gravitational waves from massive black hole binaries in LISA

Connor R. Weaving1∗ , Laura K. Nuttall1 , Ian W. Harry1 , Shichao Wu2 and
Alexander Nitz3,2 0009-0008-2697-2998 0000-0002-8599-8791 0000-0002-5304-9372
0000-0002-9188-5435 0000-0002-1850-4587 1 University of Portsmouth,
Portsmouth, PO1 3FX, United Kingdom  2 Max-Planck-Institut für
Gravitationsphysik (Albert-Einstein-Institut), D-30167 Hannover, Germany  3
Department of Physics, Syracuse University, Syracuse NY 13244, USA  ∗
Corresponding author

###### Abstract

The Laser Interferometer Space Antenna (LISA), due for launch in the mid
2030s, is expected to observe gravitational waves (GW)s from merging massive
black hole binaries (MBHB)s. These signals can last from days to months,
depending on the masses of the black holes, and are expected to be observed
with high signal to noise ratios (SNR)s out to high redshifts. We have adapted
the PyCBC software package to enable a template bank search and inference of
GWs from MBHBs. The pipeline is tested on the LISA data challenge (LDC)’s
Challenge 2a (“Sangria”), which contains MBHBs and thousands of galactic
binaries (GBs) in simulated instrumental LISA noise. Our search identifies all
6 MBHB signals with more than $92\%$ of the optimal SNR. The subsequent
parameter inference step recovers the masses and spins within their $90\%$
confidence interval. Sky position parameters have 8 high likelihood modes
which are recovered but often our posteriors favour the incorrect sky mode. We
observe that the addition of GBs biases the parameter recovery of masses and
spins away from the injected values, reinforcing the need for a global fit
pipeline which will simultaneously fit the parameters of the GB signals before
estimating the parameters of MBHBs.

##  1 Introduction

In preparation for the Laser Interferometer Space Antenna (LISA) [1] launch in
the mid 2030’s, considerable efforts are being made to design and develop the
data analysis tools needed to extract and examine the numerous different types
of gravitational wave (GW) signals expected to be found in the LISA band. LISA
will probe the mHz GW spectrum, allowing us to explore sources such as massive
black hole binaries (MBHB) with total mass between $10^{5}-10^{8}M_{\odot}$
[2]; intermediate-mass black hole binaries (IMBHB) with total masses in the
range $10^{2}-10^{5}M_{\odot}$ [3]; extreme mass-ratio and intermediate mass-
ratio inspirals (EMRI and IMRI) where coalescences have mass ratios of
$10^{-6}-10^{-3}$ and $10^{-3}-10^{-1}$ respectively, with total masses in the
range $10^{3}-10^{7}M_{\odot}$ [4, 5], and more. Each of these systems, if
detected and analysed properly, will help constrain and test a wide range of
theories. For example, the effect of dark matter on a GW propagating from a
MBHB or constraining the deviation from the Kerr metric using EMRIs (for
example see Ref [6]). Depending on the model of MBHB populations, the merger
rate is estimated to be from O(2) - O(100) $yr^{-1}$ [7, 8, 9, 10].

In order to ensure that we can detect and analyse these signals, the LISA data
challenge (LDC) [11, 12] generates mock LISA data with various GW signals
injected into simulated LISA noise. The series of challenges will eventually
replicate all of the complexities that the true LISA data will contain, such
as gaps in the data; noise artefacts found within the data that will affect
the ability to infer a system’s parameters; realistic LISA spacecraft orbits;
a time varying estimate of the noise to account for longer, GW signals. All of
these issues will need to be addressed before LISA data can be analysed.
Several groups are working on different aspects of LISA data analysis and
their work is highlighted in Refs [13, 14, 15, 16, 17, 18, 19, 20, 21, 22].

In this paper, we present our pipeline which was developed to analyse the
second LDC, referred to as “Sangria” [23, 24]. The challenge consists of a
training and blind data set where both the signals and the Gaussian
instrumental noise are generated in the same way. We used the training data
set to test and develop our pipeline as the number of signals with their
injected parameters were given with the data. We then ran our analysis on the
blind data, which contains an unknown number of signals. Our results for the
blind data challenge were submitted to the LDC in October 2022. The blind
challenge has now finished and the injected parameter values have been
released publicly and can be accessed in Ref [25].

There are two types of GW signals injected into Sangria. The first type of
signal is galactic white dwarf binaries (GBs) [26] which are expected to be
the most abundant source of GW signals present in LISA data. The accumulation
of these signals forms a “confusion noise” which dominates the $0.5-3$ mHz
range of the frequency spectrum. A subset of these are known as verification
binaries (VBs) [27], which are signals that are known in advance due to
electromagnetic observations. The second type of signal in the data are MBHBs.
These are generated using the IMRPhenomD waveform model [28, 29]. IMRPhenomD
produces two GW polarizations which encode the effect a GW has on the relative
distance between particles, the derivation of which is well known and can be
found in Ref [30]. These polarizations are then projected onto the LISA
constellation to produce an arm response. The arm response from LISA is known
to be dominated by laser frequency noise and therefore Time Delay
Interferometry (TDI) [31, 32] techniques must be used in order to extract GW
information from the raw data streams. The Michelson TDI variables X, Y and Z
are provided in this data challenge.

The data set spans a year’s worth of data sampled at 5 seconds, which is the
expected sample rate for LISA. Figure 1 shows a time series of the X TDI data
stream of the Sangria training data set. An overview of previous LDCs can be
found in Refs [33, 34, 35].

Figure 1: Timeseries of X TDI variable from the Sangria training data set. The
data consists of a year’s worth of data sampled at 5s, with 15 MBHB signals
injected into the data. This data also contains Gaussian instrumental noise
and 30 million GBs, of which 17 are VBs. This image was taken from the LDC
website for the Sangria challenge in Ref [23].

The tools used in our analysis have been developed within PyCBC [36]. PyCBC is
an open source GW data analysis software package and is one of the main
pipelines used to detect and characterise compact binary mergers from LIGO,
Virgo and KAGRA (LVK) detectors [37, 38, 39, 40]. PyCBC has been used to
detect and analyse $\mathcal{O}(100)$ GW signals from compact binary mergers
to date, an overview of which can be found in Refs [41, 42, 43]. We have
adapted the techniques developed for LVK searches and applied them to LISA-
like signals, in this case MBHBs.

It is important to note that there is far more work to do to achieve a
reliable analysis of MBHBs in realistic conditions. The tools made available
through this paper are a first attempt at the bigger picture solution, which
will require many more iterations before being considered “production ready”.

The paper is laid out as follows. In Section 2, we introduce PyCBC and the
components which make up our pipeline. We give the basic details needed to
understand how PyCBC approaches detecting and recovering GW system parameters.
In Section 3, we outline how we have adapted PyCBC in order to detect and
infer MBHBs. In Section 4, we present the results our search and each stage of
inference leading to the analysis of the blind Sangria data. Finally in
Section 5 we discuss the challenges in developing these tools, future work and
our conclusions.

##  2 PyCBC: Search and inference

We can split the PyCBC analysis pipeline up into two distinct sections, the
search and the inference. The search portion of the pipeline is responsible
for identifying candidate GW signals within the data, and the inference
portion returns a statistical estimate of the parameters of the GW signals
identified. Below we outline which tools within PyCBC have been used and
adapted for MBHB signals.

###  2.1 The search

The main component of the search is the application of matched filtering. The
matched filter signal-to-noise ratio (SNR) is designed to indicate when a
particular signal we are searching for could be present in data [44]. For data
in the form $d(t)=h(t)+n(t)$, where $h(t)$ is a modelled signal we are looking
for and $n(t)$ is the detector noise, we can choose a linear filter, $F(t)$,
such that the output SNR is maximised for the signal $h(t)$. We can represent
this filtering acting in the time or frequency domain, as shown in Equation 1,

$\phi=\int^{\infty}_{-\infty}F(t)d(t)dt=\int^{\infty}_{-\infty}\tilde{F}^{*}(f)\tilde{d}(f)df$
(1)

where tilde represents the Fourier transform. The optimal SNR defined as
$\rho_{opt}=S/N$, where S is the value of $\phi$ where a signal $h(t)$ is
present, and N is the root mean squared value of $\phi$ with no signal.

If we assume the noise to be Gaussian and stationary, we can define the the
one-sided power spectral density (PSD) via $S_{n}$ in Equation 2,

$\langle\tilde{n}(f)\tilde{n}^{*}(f^{{}^{\prime}})\rangle=\frac{1}{2}S_{n}(|f|)\delta(f-f^{{}^{\prime}}).$
(2)

where $f\geq 0$ and the square brackets represents the ensemble average. By
introducing the inner product of two real time series $a(t)$ and $b(t)$ in
Equation 3

$(a|b)=4\Re\int^{\infty}_{0}\frac{\tilde{a}^{*}(f)\tilde{b}(f)}{S_{n}(f)}df,$
(3)

we can express the optimal SNR as

$\rho_{opt}=\frac{(u|h)}{(u|u)^{1/2}},\quad\textrm{where}\quad\tilde{u}(f)=\frac{1}{2}\tilde{F}(f)S_{n}(f).$
(4)

In order to maximise $\rho_{opt}$ it is clear that $\tilde{F}(f)$ should be
proportional to $h(t)$

$\tilde{F}(f)=A\frac{\tilde{h}(f)}{S_{n}(f)},$ (5)

where A is some arbitrary constant. This results in the optimal SNR taking the
form

$\rho_{opt}=(h|h)^{1/2}.$ (6)

In order to filter the data for a given signal, we would compute the matched
filter SNR in Equation 7 which is distinct from $\rho_{opt}$

$\rho_{mf}=\frac{(d|h)}{(h|h)^{1/2}}.$ (7)

To generate an SNR timeseries, we can inverse Fourier transform the matched
filter as shown in Equation 8,

$(d|h)(t)=\int^{\infty}_{-\infty}\frac{\tilde{d}(f)|\tilde{h}(f)|^{*}}{S_{n}(f)}\exp[-2\pi
ift]df,$ (8)

where $t$ is time. Equation 8 shifts the waveform $h(t)$ to the defined time.
By computing across the entire time of the data, the SNR timeseries returned
indicates at which times the signal $h(t)$ is likely to be present by peaks in
the SNR timeseries.

We use a model for signals, $h(t)$, which we will now refer to as templates,
that represents the response of a detector given a GW passing through the
instrument. PyCBC can project the GW polarizations produced by varying
waveform models onto different ground-based detectors to return an arm
response. For some GW emitting systems, many parameters are required to define
the full signal response. Given such a broad parameter space of values to
explore, we require a large number of templates to search the data in order to
potentially detect a GW. PyCBC generates template banks which can contain
$\mathcal{O}(100,000)$ templates for a LVK-like search to explore the
parameters which affect GW emission from compact binary coalescing systems.
These parameters are split into intrinsic, which define the masses of both
compact objects and their spins, and extrinsic, which describes the system’s
sky position, the luminosity distance to the binary and the orientation of the
binary relative to the detector [45]. Details of the various methods used to
generate template banks can be found in Refs [46, 47, 48, 49, 50].

When generating a template bank, the match, as defined in Equation 9, is
computed between templates.

$\mathcal{M}(h_{1},h_{2})=\mathrm{max}_{\phi_{c},t_{c}}\left[(h_{1}|h_{2})\right],$
(9)

given two templates $h_{1}$ and $h_{2}$, the match tells us how similar to
each other these templates are whilst maximising over the coalescence time
$t_{c}$ and the phase at coalescence $\phi_{c}$. Templates are first
normalised such that the inner product gives $(h_{1}|h_{1})=1$. Using the
normalised templates in the match returns a number in the range [0,1] where 1
states that the two normalised templates are identical.

The matched filter is then applied to all templates and the data we are
searching over, giving an SNR timeseries per template. Each timeseries is then
passed to a clustering algorithm, to return the highest SNR template per time
in the data. A trigger is defined by a peak in the SNR timeseries reaching, or
surpassing, a certain threshold. Additional steps are then applied to search
for coincidences over multiple detectors and calculate their significance [37,
51].

###  2.2 The inference

PyCBC can infer signal parameters through sampling the posterior probability
defined by Bayes theorem

$p(\underline{\theta}|d)=\frac{\mathcal{L}(d|\underline{\theta})p(\underline{\theta})}{p(d)}$
(10)

where $\mathcal{L}(d|\underline{\theta})$ is the likelihood of the data $d$
given the parameters $\underline{\theta}$, $p(\underline{\theta})$ is the
prior of the parameters and $p(d)$ is the evidence and $\underline{\theta}$
represents an integer $N$ number of parameters
$\underline{\theta}=\\{\theta_{1},...,\theta_{N}\\}$. The main component of
Equation 10 to compute is the likelihood which under the assumptions of
Gaussianity and stationarity takes the form:

$\log\mathcal{L}(d,\underline{\theta})\propto(d|h_{\theta})-\frac{1}{2}(h_{\theta}|h_{\theta})-\frac{1}{2}(d|d),$
(11)

where $h_{\theta}\equiv h(\underline{\theta})$.

We can speed up likelihood calculations by enforcing some assumptions about
our GW.

$\tilde{h}(f;\underline{\theta},\phi)=A(f,\underline{\theta})\exp[i(\Psi(f,\underline{\theta})+\phi)].$
(12)

We use IMRPhenomD (2,2) mode only signals, and express our GW strain,
$\tilde{h}(f)$, as a product of smoothly varying amplitude and phase factors,
where $f$ is frequency, $A(f,\underline{\theta})$ is the amplitude of the GW,
$\Psi$ is the phase of the waveform and $\phi$ is some arbitrary phase
constant.

With the GW strain expressed in this form, we use the application of
heterodyning (sometimes referred to as relative binning) which assumes a
smooth signal which varies in amplitude and phase. Heterodyning is a signal
processing technique used to separate data into slowly and rapidly varying
frequency components. Our inference application of heterodyning allows us to
greatly speed up likelihood calculations by first finding a high likelihood
template with our template bank search. The corresponding template is then
used as a reference to explore the region of high likelihood. A derivation of
this method applied to GWs can be found in Refs [52, 53, 54].

In order to explore the parameter space, there are several samplers in PyCBC
which implement different methods, including Markov Chain Monte Carlo (MCMC)
and Nested Sampling. An overview of PyCBC’s inference setup can be found in
Ref [55].

##  3 Adapting PyCBC’s search and inference for MBHBs

###  3.1 Templates

In order to perform a matched filter search, we need to be able to generate
templates that represent LISA’s arm response given a GW signal passing through
the instrument. The data we are searching over is given in terms of TDI
variables; we therefore must put the templates in this format.

A vital part of this development is the use and implementation of the open
source software BBHx waveform generator [56, 57, 58, 59, 60, 61, 62, 63],
which we use to generate the expected strain from MBHB systems. BBHx generates
the GW polarizations, projects them onto LISA to form an arm response and then
creates the TDI channels. BBHx implements the IMRPhenomHM [64, 65], an aligned
spin waveform model which produces GW polarizations based on IMRPhenomD, but
includes other sub-dominant modes. In this work, we set the IMRPhenomHM to
generate just the (2,2) mode in order to match the Sangria convention. We
transform both the data and the TDI templates into noise independent data
streams referred to as A, E and T [66] via the transforms in Equation 13.

$A=(Z-X)/\sqrt{2},\quad E=(X-2Y+Z)/\sqrt{6},\quad T=(X+Y+Z)/\sqrt{3}.$ (13)

The data streams A and E have equal SNR across the entire LISA frequency band.
This is not true for the T TDI stream as at lower frequencies the SNR
contribution from T does not compete with A and E. However, T does contribute
SNRs that can be greater than A or E at higher frequencies in the LISA band
and is therefore included in the inference analysis [67]. BBHx generates
templates A, E and T in the same convention as in Equation 13.

###  3.2 Template bank

This work explores whether a template bank based search is a viable approach
to finding MBHBs in LISA data, so to reduce complexity we only perform the
search on the A TDI stream. We generated a template bank for the A TDI data
stream, which was used for the filtering stage of the search on just the A TDI
data. As the A and E TDI streams have equal SNR across the LISA frequency
band, searching through just the A or just E data streams would be equivalent.
We do lose power when searching over one instead of both of the TDI streams,
but due to the SNRs of these signals being so large, we find that the
information lost is negligible. Due to the SNRs of MBHBs being so large, we
populated a less dense template bank when compared with LIGO template bank
searches to balance accuracy of the signal trigger time versus computation
cost of producing the bank.

For the bank generation we use stochastic placement [68, 69, 70]. Specifically
we use the PyCBC algorithm “brute bank” with the parameters defined in Table 1
and a minimal match of 0.9 which resulted in a bank size of 50 templates. The
minimal match defines how densely the parameter space is sampled. For example,
with a minimal match of 0.9, if any newly generated template matches with any
other template in the bank above the minimal match, it is rejected and another
template is generated. The process repeats until a pre defined number of
templates have been rejected from the bank placement.

Parameter | Lower Bound | Upper Bound
---|---|---
$M_{T}$ ($M_{\odot}$) | $2\times 10^{5}$ | $2\times 10^{7}$
$q$ | 1 | 4
$a_{1}$, $a_{2}$ | -0.99 | 0.99
$\beta$ (rad) | $-\pi/2$ | $\pi/2$
$\lambda$ (rad) | 0 | $2\pi$
$\psi$, $\iota$ (rad) | 0 | $\pi$
Table 1:  The parameter ranges used to generate the template bank. This covers
the full range of potential systems that are defined within the conventions of
the Sangria challenge. $M_{T}$ is detector frame total mass of the black
holes; $q$ is mass ratio where $m_{1}>m_{2}$ and $m_{1}$ is the primary mass
and $m_{2}$ the secondary; $a_{1}$ and $a_{2}$ are the aligned spins of the
two black holes; $\beta$ is ecliptic latitude; $\lambda$ is ecliptic
longitude; $\psi$ is polarization angle; $\iota$ is inclination.

The luminosity distance $D_{L}$ is marginalised over in the matched filter
calculations and is set to an arbitrary value when generating templates. The
time of coalescence $t_{c}$ is sampled with one month from the beginning and
end of the year removed. This is to ensure no edge corruption artefacts are
present in the data when performing matched filtering, which could cause false
triggers.

When computing the match, maximising Equation 9 over time of coalescence
$t_{c}$ implies that we do not take into consideration the varying response
function from LISA over time. LISA’s response function varies over time due to
the antennas sky position changing within the duration of the MBHB signal.
However, because the majority of the SNR from merging MBHBs is within the
final $\mathcal{O}(hour)$, and these systems have SNRs $\mathcal{O}(1000)$, we
are safe to assume LISA is approximately stationary during this time period.

As $t_{c}$ is directly related to LISA’s orbital position in the case of BBHx,
we ensure that the template bank can generate templates at different orbital
positions by varying $t_{c}$. The actual recovery of $t_{c}$ is done via the
filter output time series from the inner product, shown in Equation 8, as
$t_{c}$ is not an explicit template bank parameter.

###  3.3 Sampling, priors and sky folding

We sample the heterodyne likelihood using the Nested Sampling tool Dynesty
[71]. Table 2 outlines the parameters used with their corresponding
distributions in the inference. $M_{c}$ is the chirp mass of the system [72].

Parameter | Lower Bound | Upper Bound
---|---|---
$M_{c}$ ($M_{\odot}$) | 87055 | 8705505
$q$ | 1 | 4
$a_{1}$, $a_{2}$ | -0.99 | 0.99
$t_{c}$ (s) |  $t_{ref}$ \- 2 days |  $t_{ref}$ \+ 2 days
$\beta$, $\lambda$, $\psi$ (rad) | 0 | $\pi/2$
$\iota$ (rad) | 0 | $\pi$
$D_{L}$ (Mpc) | $10^{3}$ | $10^{6}$
Table 2: A list of all parameters and their ranges given to the prior,
$p(\underline{\theta})$. $M_{c}$, $q$, $a_{1}$, $a_{2}$, $t_{c}$, $\psi$ and
$D_{L}$ are distributed uniformly. $\beta$, $\lambda$ and $\iota$ are
distributed isotropically. $t_{ref}$ is the coalescence time found in the
search.

The importance of multi-modality and degeneracies when considering extrinsic
parameters are summarised in detail in Ref [60] and are vital to understand
for our application of sampling we introduce here. Considering a (2,2) mode
signal within the LISA frame, and applying the frozen low-f approximation
which neglects both time and frequency dependence in the instrument response,
there are 8 exact degeneracies for sky location. Four longitudinal at
$\\{\lambda_{L}+\frac{\pi}{2}(0,1,2,3),\psi_{L}+\frac{\pi}{2}(0,1,2,3)\\}$ and
two latitudinal at $\\{\pm\beta_{L},\pm\cos{\iota_{L}},\pm\cos{\psi_{L}}\\}$,
where the subscript L represents the extrinsic parameters in the LISA
reference frame.

$\eqalign{\lambda^{(r)}_{L}&=\lambda_{L},\\\ \beta^{(r)}_{L}&=-\beta_{L},\\\
\psi^{(r)}_{L}&=\pi-\psi_{L},\\\ \iota^{(r)}_{L}&=\pi-\iota_{L},\\\
\phi^{(r)}_{L}&=\phi_{c},}$ (14)

The BBHx waveform model that we use in this work, uses the frozen
approximation and frequency-dependent instrument response which neglects
LISA’s motion. This breaks the 8-fold degeneracy into an approximate 2-fold
degeneracy in the reflected mode shown in Equation 14.

Instead of sampling the entire parameter space, we restrict our sky location
parameters to one of eight octants of the sky. After drawing a sample from our
defined octant, we compute the likelihood at the other 7 modes. We marginalise
over the 8 sky modes by summation during the inference process. We then
normalise them to sum to one. A post-processing step is then implemented to
“unfold” the sky positions by randomly selecting one of the 8 sky modes based
on the likelihood of each. This effectively samples the entire parameter space
without having to sample other parts of the sky as we compute the likelihood
directly at the analytical location of the other sky modes.

All samples for inference are taken within the LISA reference frame, not the
solar system barycentre (SSB) frame, even though the injected parameters in
Sangria are all with reference to the SSB frame. The techniques described
above for sampling the sky position are used in every inference run in this
work.

###  3.4 PSDs

When GBs are included in the data, we assume them to be part of the noise
$n(t)$. As there are several thousands of these signals of varying amplitude
and phase, we can justify their accumulation to be Gaussian (but not
stationary) by the central limit theorem. When constructing the PSDs for all
three TDI variables, we implemented Welch’s method over 31 segments of the
year’s worth of data. This gave frequency bins with a resolution of $10^{-6}$
Hz, which were narrow enough to model the structure of the PSD accurately. The
same number of segments is used to generate all other PSDs for the different
inference stages and the search. We generate a PSD per TDI channel at each
inference stage.

By constructing the PSDs in this way we have ignored the fact that the PSDs
are time-dependent and therefore cannot be stationary. This is also true in
the case where the GBs are not included in the noise due to LISA’s motion
across the sky. We prioritised the LISA high-frequency information over the
entire year’s worth of data by reading all of the data and using 31 segments
instead of a quasi-stationary PSD which is estimated on smaller chunks of
data.

.

Figure 2: A figure showing the different PSDs we use in the results section.
The red curve is the A TDI PSD of the Sangria data and the black curve is the
A TDI PSD of Sangria without GBs.

##  4 Results

In this section we present the results of our search and multiple inference
stages. The associated data release which contains all of the code and
instructions on how to reproduce the results presented can be found in Ref
[73]. It contains the results for all analysis of the search and inference in
all of their various stages. Documentation for ongoing developments in PyCBC
for LISA related searches and inference can be found in Ref [74].

###  4.1 Search results

To quantify how effective the template bank search was, we computed the match
between the highest SNR template per trigger against the injected waveform
parameters. All matches are greater than 0.92 and are presented in Table 3 for
each candidate.

Signal index | Match
---|---
0 | 0.94510
1 | 0.92666
2 | 0.97344
3 | 0.94091
4 | 0.94562
5 | 0.92439
Table 3: The match defined in Equation 9 between the highest SNR template per
trigger and the waveforms injected value of the template bank search over the
Sangria blind data.

The template parameters from the search and the corresponding injected
waveform parameters for signal 0 can be found in Table 4. We attempted to
produce template banks at higher minimal matches but the time to complete was
$\mathcal{O}(days)$ instead of $\mathcal{O}(minutes)$ in the 0.9 minimal match
case. A template bank search as we have applied it worked exactly how we
needed it when considering the balance of computational cost verses accuracy
of the returned template.

Parameters | $\csvcolii$ | $\csvcoliii$
---|---|---

Table 4: Comparison with template parameters found from the search verses the
injected waveform parameters of signal 0 in the Sangria blind data.

###  4.2 Inference results

In this section, we present our results from the inference. All of the runs
have the reference parameters for the heterodyne likelihood set to the
injected parameter values. Each inference run uses all three TDI data streams.
Running the heterodyne likelihood with the reference values found in the
search gave no significant bias when comparing to the reference parameters
being set to the injected signal values.

When building the inference tools for this challenge, there were several
stages of additional complexities which were successively built upon to ensure
that each step performed as expected. We performed inference with increasing
complexity as follows:

  * •

Injections with BBHx in zero noise.

  * •

Injections with BBHx in Sangria training noise without GBs.

  * •

Injections with BBHx in Sangria training noise with Sangria training GBs.

  * •

Sangria blind data.

.

Figure 3: A corner plot showing the recovered intrinsic parameters posteriors
of signal 0 in zero noise. The blue curves represent a year’s worth of data
containing only signal 0. The orange curves represent data containing all 6
signals. The red lines represent the injected parameter values of the signal.
The tiles at the top of each column display the median and $90\%$ range for
each inference. The contour in the 2D tiles represents the $90\%$ confidence
intervals.

When generating data with BBHx injections for the with and without GB cases,
we used the Sangria training data set. The Gaussian LISA noise and GBs for
Sangria blind are generated in the same way as the training data set, but are
unique to the blind data set. The Sangria data has waveform injections that
were created with the LDC code base and not injected with BBHx.

We will now present the results of each stage of this development. The
parameters used to create the BBHx injections are identical to that of the
signals in the Sangria blind data. We ran analysis on all signals injected
with BBHx and Sangria blind, but will concentrate our discussion on the first
signal (from time zero) which we refer to as signal 0. We will distinguish the
Sangria blind injection of the first signal as “blind 0” as it is not produced
with BBHx and is in the Sangria blind noise and GBs data. Quoted parameters
found in the inference are the maximum log-likelihood sample from the
posterior.

All of the BBHx inference runs were performed with the heterodyne parameters
set to the true injected parameters. The final Sangria runs reference
parameters were set to the parameters found in the search section of the
pipeline. Figure 2 shows the PSDs used for the first 3 stages before analysing
the Sangria blind data. The same sampler settings are used for each analysis
to ensure a fair comparison.

.

Figure 4: Extrinsic parameter posterior corner plot of injected BBHx signal in
zero noise for signal 0. Red lines represent the injected parameter values of
the signal. The tiles at the top of each column display the median and $90\%$
range for each inference. The contour in the 2D tiles represents the $90\%$
confidence intervals.

For each stage of development, we performed two sets of inference. One set of
inference contained only signal 0 injected into the data, shown in blue in
subsequent Figures, and the other contained all 6 signals from Sangria, shown
in orange. All injections are created in a year’s worth of data sampled at 5s
intervals to match the convention of the Sangria challenge.

Note that all of these parameters are given in the detector frame. This is the
standard return of BBHx. For the interest of those wanting to inspect
parameters in the source frame, the redshift for signal 0 injections is
$z=3.72$. Finally, we highlight the degeneracy,
$h(\underline{\theta},\psi)=h(\underline{\theta},\psi+\pi)$, which can be
observed throughout the results section.

####  4.2.1 Zero noise, BBHx injection -

Figure 3 presents the recovered intrinsic parameter posteriors for this
analysis, and Figure 4 the extrinsic parameters, along with their $90\%$
confidence intervals. The black PSD curve from Figure 2, was used to scale the
signals, to give SNRs of comparable sizes to subsequent analyses. We can see
for the single and multiple injection cases that the intrinsic parameters are
recovered within the $90\%$ confidence interval of the injected values.
However, most extrinsic parameters do not fall within the $90\%$ confidence
interval; the sky position $\beta$ has found the correct sky mode but the
injected value falls just outside the peak of both single and multiple signal
runs. The correct mode for $\lambda$ is also recovered, however, in both the
single and multiple signal cases the $90\%$ confidence intervals do not
capture the injected value.

.

Figure 5: Intrinsic parameter posterior corner plot of injected BBHx signal(s)
in Gaussian stationary noise without GBs for signal 0. Red lines represent the
injected parameter values of the signal. The tiles at the top of each column
display the median and $90\%$ range for each inference. The contour in the 2D
tiles represents the $90\%$ confidence intervals.

####  4.2.2 Sangria noise without GBs, BBHx injections -

These signals are injected into Gaussian, stationary noise using the Sangria
training data. The black curve in Figure 2 is the PSD used for this stage.

As Figure 5 shows, the inclusion of Gaussian noise when comparing to the zero
noise case for both multiple and single injections has no negative bias to the
inferred intrinsic parameters. All intrinsic parameters fall within the $90\%$
confidence interval in both cases. The extrinsic parameter recovery, shown in
Figure 6, shows no notable differences when comparing to the zero noise case.
We do start to see some preference to other likelihood modes in sky position
as in two cases (signals 2 and 4) the inference preferred the other high
likelihood mode of $\beta$ instead of the true injected mode. This is true in
both the single injection case and the multiple injection case.

.

Figure 6: Extrinsic parameter posterior corner plot of injected BBHx signal(s)
in Gaussian stationary noise without GBs for signal 0. Red lines represent the
injected parameter values of the signal. The tiles at the top of each column
display the median and $90\%$ range for each inference. The contour in the 2D
tiles represents the $90\%$ confidence intervals.

####  4.2.3 Sangria noise with GBs, BBHx injections -

The red curve in Figure 2 is the PSD used for this stage and contains all of
the noise from the GBs. We see that this significantly changes the shape of
the PSD in the $10^{-3}-10^{-2}Hz$ region of frequency space.

.

Figure 7: Intrinsic parameter posterior corner plot of injected BBHx signal(s)
in Gaussian stationary noise with GBs for signal 0. Red lines represent the
injected parameter values of the signal. The tiles at the top of each column
display the median and $90\%$ range for each inference. The contour in the 2D
tiles represents the $90\%$ confidence intervals.

The posteriors in Figures 7 and 8 are broader which is expected given the
noise curve with GBs present. Figure 7 demonstrates that the intrinsic
parameters still fall within the $90\%$ confidence interval.

Signal | Sangria noise without GBs | Sangria noise with GBs | Sangria blind
---|---|---|---
0 | 78% | 34% | 93%
1 | 41% | 84% | 98%
2 | 20% | 44% | 99%
3 | 2.7% | 84% | 100%
4 | 36% | 84% | 81%
5 | 18% | 6.7% | 86%
Table 5: Each percentage represents the amount of samples whose likelihood
value is greater than the likelihood of the corresponding injected waveform.

GBs have had the most significant effect on inference thus far as the peaks in
the posterior likelihood have shifted from the injected value when compared to
the no GBs case. We also see a lot more structure in the posteriors that
includes GBs. In Table 5 we have calculated the percentage of samples from the
posterior with a higher likelihood than that of the likelihood of the injected
waveform. If the assumptions of Gaussianity and stationarity are valid, these
values for a collection of posteriors from different signals/injections should
be uniformly distributed. We see for the case of Gaussian stationary noise
without GBs that the distribution of values does appear to be uniform. The
test performed here would have benefited from $\mathcal{O}(100)$ posteriors to
properly conclude that the distribution is uniform. When including GBs, we can
see that the distribution is no longer uniform. The effect of GBs being non-
stationary therefore has biased the inferred parameters away from the injected
value. Even with only six signals we can clearly see the bias when the data
contains GBs and when inspecting the full Sangria data. We have presented the
other 5 signal parameter injections in Appendix A. We can clearly see that
Figures A1-A5 show a bias in intrinsic parameter recovery as signal 0.

.

Figure 8: Extrinsic parameter posterior corner plot of injected BBHx signal(s)
in Gaussian stationary noise with GBs for signal 0. Red lines represent the
injected parameter values of the signal. The tiles at the top of each column
display the median and $90\%$ range for each inference. The contour in the 2D
tiles represents the $90\%$ confidence intervals.

As for the extrinsic parameters for signal 0, the inferred values represent
the injected waveform well. This is not true for the other signals in the data
as most struggle to find the correct sky mode for $\beta$ and/or $\lambda$.

Blind 0
---
Parameters | $\csvcolii$ | $\csvcoliii$

Table 6: The injected and recovered parameters for signal 0 in the blind
Sangria dataset. The recovered values are stated with the 90% confidence
regions.

.

Figure 9: Intrinsic parameter posterior corner plot of Sangria blind data for
signal 0. Red lines represent the injected parameter values of the signal. The
tiles at the top of each column display the median and $90\%$ range for each
inference. The contour in the 2D tiles represents the $90\%$ confidence
intervals.

####  4.2.4 Sangria blind -

The full parameter recovery of blind 0 in Sangria is presented in Table 6,
which also displays the injected parameter values. The intrinsic parameters in
Figure 9 fall within the $90\%$ confidence interval. However, this is not true
for all signals in the data set. There is one case where all intrinsic
parameters peak at roughly the $88\%$ confidence interval. Overall, the peaks
in likelihood are far less prominent when compared to the previous cases.

Extrinsic parameters also have far less prominent peaks than previous cases
and suffer from the same sky position degeneracy issues where the wrong mode
is commonly favoured. For half of the cases, $D_{L}$ is recovered within the
$90\%$ interval and for all cases, $t_{c}$ and $\iota$ are recovered in the
$90\%$ confidence interval.

.

Figure 10: Extrinsic parameter posterior corner plot of Sangria blind data for
signal 0. Red lines represent the injected parameter values of the signal. The
tiles at the top of each column display the median and $90\%$ range for each
inference. The contour in the 2D tiles represents the $90\%$ confidence
intervals.

Figure 10 shows the extrinsic recovery of our reference signal and is the only
example out of the six signals that recovers the correct mode for $\lambda$.
Extrinsic parameter recovery would be greatly improved when analysing signals
with higher order modes [60]. As this data challenge uses just (2,2) mode
injections, it was expected that the extrinsic parameters would be the most
difficult to recover correctly.

##  5 Conclusion

The application of a template bank search was able to identify all MBHB
signals in the Sangria mock data challenge. Although the search for MBHB
signals is relatively trivial due to the high SNR, being able to provide
reference parameters for the heterodyne likelihood is the main advantage when
taking this approach. We can produce a template bank $\mathcal{O}(minutes)$ to
fulfill a minimal match of 0.9. Given that the SNR recovered is greater than
$92\%$, we find that this technique is a viable approach to quickly search
through LISA data for MBHBs. At the same time, the templates returned from the
search provide appropriately high likelihood reference parameters for
heterodyning which do not bias parameter estimation. To expand and improve
upon the search, we could consider a coincident search over both the A and E
TDI data streams. This would filter each TDI stream separately and then check
for coincident triggers.

We have shown that the introduction of different forms of noise can have an
effect on both intrinsic and extrinsic parameter recovery, although the
extrinsic parameter recovery suffers more than the intrinsic. Something to
consider about this data set is how the injected signals only contain the
dominant (2,2) mode. Other works have shown that sky position recovery is
vastly improved as other sub dominant modes are introduced to parameter
estimation [60]. A further study could be to investigate the inference of a
multi-modal signal, in the presence of LISA instrumental noise. We could then
go further by introducing non-stationarities in the data, with the inclusion
of GBs.

The intrinsic parameter recovery of MBHBs in the presence of confusion noise
caused by GBs can bias the result away from the true value. Although the
injected values in our work did fall within the $90\%$ confidence intervals,
the peak likelihood has moved away from the injected value when compared to
the no GB inference case, see Figures 5 and 7. This motivates the need for a
full global fit pipeline, that can identify and infer the properties of the
GBs and then remove them from the data before running inference on MBHBs. This
was first described in Ref [13].

An area of further investigation concerns how we estimate PSDs. As stated
previously, in this work we estimate our PSDs by reading in the entire year’s
worth of data. Due to LISA’s motion across the sky, this makes the PSDs non-
stationary. One possible method is to construct our PSD by assuming quasi-
stationarity on smaller time chunks of the year’s worth of data, and assemble
these into one PSD using the time/frequency track of the reference signal in
question. It would also be interesting to investigate what effects the PSD
being non-stationary has on the search and inference of such high SNR signals.

##  6 Acknowledgements

We would like to thank Stas Babak, Quentin Baghi, Maude Le Jeune, Jean-
Baptiste Bayle, and Natalia Korsakova for helpful discussions regarding LISA
data analysis. We also thank Michael Katz for help with using BBHx and for
making his code publicly available for use. We also thank the GW group at the
University of Portsmouth for their discussion and support during the project.
CW would like to particularly thank Gareth Cabourn Davies, Arthur Tolley and
Ellie Wadman for their advice, feedback and support. We are grateful to the
computing teams from the University of Portsmouth and AEI Hannover for their
significant technical support. CW is supported by a STFC studentship. LKN
acknowledges support from the UKRI Future Leaders Fellowship through grant
MR/T01881X/1. IWH and LKN acknowledge support from the UK Space Agency through
grant ST/X002225/1. Numerical computations were done on the Sciama High
Performance Compute (HPC) cluster which is supported by the ICG, SEPNet and
the University of Portsmouth.

## References

## References

  * [1] Amaro-Seoane P et al (LISA) 2017 (Preprint 1702.00786)
  * [2] Wise J H 2023 (Preprint 2304.09311)
  * [3] Strokov V, Fragione G and Berti E 2023 (Preprint 2303.00015)
  * [4] Barack L and Pound A 2019 Rept. Prog. Phys. 82 016904 (Preprint 1805.10385)
  * [5] Gair J R, Babak S, Sesana A, Amaro-Seoane P, Barausse E, Berry C P L, Berti E and Sopuerta C 2017 J. Phys. Conf. Ser. 840 012021 (Preprint 1704.00009)
  * [6] Barausse E et al. 2020 Gen. Rel. Grav. 52 81 (Preprint 2001.09793)
  * [7] Li K, Bogdanović T, Ballantyne D R and Bonetti M 2022 Astrophys. J. 933 104 (Preprint 2201.11088)
  * [8] Li K, Bogdanović T, Ballantyne D R and Bonetti M 2022 (Preprint 2207.14231)
  * [9] Steinle N, Middleton H, Moore C J, Chen S, Klein A, Pratten G, Buscicchio R, Finch E and Vecchio A 2023 (Preprint 2305.05955)
  * [10] Fang Y and Yang H 2022 (Preprint 2209.14509)
  * [11] 2023 (the new) LISA data challenge - Available online at: https://lisa-ldc.lal.in2p3.fr/, last accessed on 13.06.2023
  * [12] Baghi Q (LDC Working Group) 2022 The LISA Data Challenges 56th Rencontres de Moriond on Gravitation (Preprint 2204.12142)
  * [13] Cornish N J and Crowder J 2005 Phys. Rev. D 72 043005 (Preprint gr-qc/0506059)
  * [14] Vallisneri M 2009 Class. Quant. Grav. 26 094024 (Preprint 0812.0751)
  * [15] Karnesis N, Katz M L, Korsakova N, Gair J R and Stergioulas N 2023 (Preprint 2303.02164)
  * [16] Bayle J B, Hartwig O, Lilley M, Hees A, Chapman-Bird C, Woan G and Wolf P 2023 End-to-end simulation and analysis pipeline for LISA 57th Rencontres de Moriond on Gravitation (Preprint 2305.09702)
  * [17] Cornish N J and Hellings R W 2003 Class. Quant. Grav. 20 4851–4860 (Preprint gr-qc/0306096)
  * [18] Rubbo L J, Cornish N J and Poujade O 2004 Phys. Rev. D 69 082003 (Preprint gr-qc/0311069)
  * [19] Petiteau A, Auger G, Halloin H, Jeannin O, Plagnol E, Pireaux S, Regimbau T and Vinet J Y 2008 Phys. Rev. D 77 023002 (Preprint 0802.2023)
  * [20] Littenberg T B and Cornish N J 2023 Phys. Rev. D 107 063004 (Preprint 2301.03673)
  * [21] Bayle J B and Hartwig O 2023 Phys. Rev. D 107 083019 (Preprint 2212.05351)
  * [22] Digman M C and Cornish N J 2022 Astrophys. J. 940 10 (Preprint 2206.14813)
  * [23] 2023 LISA data challenge 2a: Sangria - Available online at: https://lisa-ldc.lal.in2p3.fr/challenge2a, last accessed on 13.06.2023
  * [24] 2023 LISA data challegne 2a, Sangria manual - Available online at: https://lisa-ldc.lal.in2p3.fr/static/data/pdf/LDC-manual-Sangria.pdf, last accessed on 13.06.2023
  * [25] LISA data challenge - sangria results and data files - Available online at: https://lisa-ldc.lal.in2p3.fr/accounts/login/?next=/file
  * [26] Nelemans G, Yungelson L R and Portegies Zwart S F 2001 Astron. Astrophys. 375 890–898 (Preprint astro-ph/0105221)
  * [27] Finch E, Bartolucci G, Chucherko D, Patterson B G, Korol V, Klein A, Bandopadhyay D, Middleton H, Moore C J and Vecchio A 2023 Mon. Not. Roy. Astron. Soc. 522 5358–5373 (Preprint 2210.10812)
  * [28] Husa S, Khan S, Hannam M, Pürrer M, Ohme F, Jiménez Forteza X and Bohé A 2016 Phys. Rev. D 93 044006 (Preprint 1508.07250)
  * [29] Khan S, Husa S, Hannam M, Ohme F, Pürrer M, Jiménez Forteza X and Bohé A 2016 Phys. Rev. D 93 044007 (Preprint 1508.07253)
  * [30] Maggiore M 2008 Gravitational waves: Volume 1: Theory and experiments vol 1 (Oxford university press)
  * [31] Shaddock D A, Tinto M, Estabrook F B and Armstrong J W 2003 Phys. Rev. D 68 061303 (Preprint gr-qc/0307080)
  * [32] Tinto M, Estabrook F B and Armstrong J W 2004 Phys. Rev. D 69 082001 (Preprint gr-qc/0310017)
  * [33] Arnaud K A et al. 2006 AIP Conf. Proc. 873 619–624 (Preprint gr-qc/0609105)
  * [34] Babak S et al. (Mock LISA Data Challenge Task Force) 2010 Class. Quant. Grav. 27 084009 (Preprint 0912.0548)
  * [35] Arnaud K A et al. 2007 Class. Quant. Grav. 24 S551–S564 (Preprint gr-qc/0701170)
  * [36] PyCBC Development Team and the LIGO / Virgo Collaborations 2023 PyCBC - Available online at: https://pycbc.org/, last accessed on 13.06.2023
  * [37] Usman S A et al. 2016 Class. Quant. Grav. 33 215004 (Preprint 1508.02357)
  * [38] Biwer C M, Capano C D, De S, Cabero M, Brown D A, Nitz A H and Raymond V 2019 Publ. Astron. Soc. Pac. 131 024503 (Preprint 1807.10312)
  * [39] Nitz A, Harry I, Brown D, Biwer C M, Willis J, Canton T D, Capano C, Dent T, Pekowsky L, Williamson A R, De S, Cabero M, Machenschalk B, Macleod D, Kumar P, Reyes S, dfinstad, Pannarale F, Kumar S, Massinger T, Tápai M, Singer L, Davies G S C, Khan S, Fairhurst S, Nielsen A, Singh S, Chandra K, shasvath and veronica villa 2022 gwastro/pycbc: v2.0.2 release of pycbc URL https://doi.org/10.5281/zenodo.6324278
  * [40] Abbott B P et al. (LIGO Scientific, Virgo) 2016 Phys. Rev. D 93 122003 (Preprint 1602.03839)
  * [41] Abbott R et al. (LIGO Scientific, VIRGO) 2021 (Preprint 2108.01045)
  * [42] Abbott R et al. (LIGO Scientific, VIRGO, KAGRA) 2021 (Preprint 2111.03606)
  * [43] Nitz A H, Kumar S, Wang Y F, Kastha S, Wu S, Schäfer M, Dhurkunde R and Capano C D 2023 Astrophys. J. 946 59 (Preprint 2112.06878)
  * [44] Helstrom C W 1968 Statistical theory of signal detection (Pergamon Press)
  * [45] Allen B, Anderson W G, Brady P R, Brown D A and Creighton J D E 2012 Phys. Rev. D 85 122006 (Preprint gr-qc/0509116)
  * [46] Owen B J 1996 Phys. Rev. D 53 6749–6761 (Preprint gr-qc/9511032)
  * [47] Owen B J and Sathyaprakash B S 1999 Phys. Rev. D 60 022002 (Preprint gr-qc/9808076)
  * [48] Babak S, Balasubramanian R, Churches D, Cokelaer T and Sathyaprakash B S 2006 Class. Quant. Grav. 23 5477–5504 (Preprint gr-qc/0604037)
  * [49] Keppel D 2013 Phys. Rev. D 87 124003 (Preprint 1303.2005)
  * [50] Privitera S, Mohapatra S R P, Ajith P, Cannon K, Fotopoulos N, Frei M A, Hanna C, Weinstein A J and Whelan J T 2014 Phys. Rev. D 89 024003 (Preprint 1310.5633)
  * [51] Davies G S, Dent T, Tápai M, Harry I, McIsaac C and Nitz A H 2020 Phys. Rev. D 102 022004 (Preprint 2002.08291)
  * [52] Cornish N J 2010 (Preprint 1007.4820)
  * [53] Zackay B, Dai L and Venumadhav T 2018 (Preprint 1806.08792)
  * [54] Cornish N J 2021 Phys. Rev. D 104 104054 (Preprint 2109.02728)
  * [55] Biwer C M, Capano C D, De S, Cabero M, Brown D A, Nitz A H and Raymond V 2019 Publ. Astron. Soc. Pac. 131 024503 (Preprint 1807.10312)
  * [56] Katz M L, Marsat S, Chua A J K, Babak S and Larson S L 2020 Phys. Rev. D 102 023033 (Preprint 2005.01827)
  * [57] Katz M L 2022 Phys. Rev. D 105 044055 (Preprint 2111.01064)
  * [58] Katz M 2021 mikekatz04/bbhx: First official public release 10.5281/zenodo.5730688
  * [59] Marsat S and Baker J G 2018 (Preprint 1806.10734)
  * [60] Marsat S, Baker J G and Dal Canton T 2021 Phys. Rev. D 103 083011 (Preprint 2003.00357)
  * [61] London L et al 2018 Phys. Rev. Lett. 120 161102 (Preprint 1708.00404)
  * [62] Husa S, Khan S, Hannam M, Pürrer M, Ohme F, Jiménez Forteza X and Bohé A 2016 Phys. Rev. D 93 044006 (Preprint 1508.07250)
  * [63] Khan S, Husa S, Hannam M, Ohme F, Pürrer M, Jiménez Forteza X and Bohé A 2016 Phys. Rev. D 93 044007 (Preprint 1508.07253)
  * [64] London L et al 2018 Phys. Rev. Lett. 120 161102 (Preprint 1708.00404)
  * [65] Abbott B P et al. (LIGO, VIRGO) 2017 Observation of Gravitational Waves from a Binary Black Hole Merger pp 291–311
  * [66] Tinto M and Dhurandhar S 2021 Living Rev. Rel 24 (Preprint 10.1007/s41114-020-00029-6)
  * [67] Tinto M and Dhurandhar S V 2005 Living Rev. Rel. 8 4 (Preprint gr-qc/0409034)
  * [68] Babak S 2008 Class. Quant. Grav. 25 195011 (Preprint 0801.4070)
  * [69] Harry I W, Allen B and Sathyaprakash B S 2009 Phys. Rev. D 80 104014 (Preprint 0908.2090)
  * [70] Ajith P, Fotopoulos N, Privitera S, Neunzert A and Weinstein A J 2014 Phys. Rev. D 89 084041 (Preprint 1210.6666)
  * [71] Speagle J S 2020 Mon. Not. Roy. Astron. Soc. 493 3132–3158 (Preprint 1904.02180)
  * [72] Blanchet L, Damour T, Iyer B R, Will C M and Wiseman A G 1995 Phys. Rev. Lett. 74 3515–3518 (Preprint gr-qc/9501027)
  * [73] Weaving, C R and Wu, S. GitHub to mbhbs-with-pycbc. https://github.com/gwastro/mbhbs-with-pycbc
  * [74] 2023 LISA parameter estimation for simulated massive black hole binary - Available online at: http://pycbc.org/pycbc/latest/html/inference/examples/lisa_smbhb_ldc_pe.html, last accessed on 13.06.2023

##  Appendix A Intrinsic parameter inference plots with GBs

.

Figure 11: Inference for BBHx injection of signal 1 from Sangria blind over
intrinsic parameters.

.

Figure 12: Inference for BBHx injection of signal 2 from Sangria blind over
intrinsic parameters.

.

Figure 13: Inference for BBHx injection of signal 3 from Sangria blind over
intrinsic parameters.

.

Figure 14: Inference for BBHx injection of signal 4 from Sangria blind over
intrinsic parameters.

.

Figure 15: Inference for BBHx injection of signal 5 from Sangria blind over
intrinsic parameters.