# Rapid Parameter Estimation for Merging Massive Black Hole Binaries Using
ODE-Based Generative Models

Bo Liang${}^{1}\,{}^{2}\,{}^{3}\,{}^{4}$, Minghui Du*1, He
Wang*${}^{4}\,{}^{6}$, Yuxiang Xu${}^{1}\,{}^{2}\,{}^{3}\,{}^{4}$, Chang Liu7,
Xiaotong Wei1, Peng Xu${}^{1}\,{}^{2}\,{}^{4}\,{}^{5}$, Li-e Qiang7, Ziren
Luo${}^{1}\,{}^{2}\,{}^{4}\,{}^{6}$ 1 Center for Gravitational Wave
Experiment, National Microgravity Laboratory, Institute of Mechanics, Chinese
Academy of Sciences, Beijing 100190, China  2 Key Laboratory of Gravitational
Wave Precision Measurement of Zhejiang Province, Hangzhou Institute for
Advanced Study, UCAS, Hangzhou 310024, China  3 Shanghai Institute of Optics
and Fine Mechanics, Chinese Academy of Sciences, Shanghai 201800, China  4
Taiji Laboratory for Gravitational Wave Universe (Beijing/Hangzhou),
University of Chinese Academy of Sciences (UCAS), Beijing 100049, China  5
Lanzhou Center of Theoretical Physics, Lanzhou University, Lanzhou 730000,
China  6 International Centre for Theoretical Physics Asia-Pacific (ICTP-AP),
University of Chinese Academy of Sciences (UCAS), Beijing 100049, China  7
National Space Science Center, Chinese Academy of Sciences, Beijing 100190,
China<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Detecting the coalescences of massive black hole binaries (MBHBs) is one of
the primary targets for space-based gravitational wave observatories such as
LISA, Taiji, and Tianqin. The fast and accurate parameter estimation of
merging MBHBs is of great significance for both astrophysics and the global
fitting of all resolvable sources. However, such analyses entail significant
computational costs. To address these challenges, inspired by the latest
progress in generative models, we proposed a novel artificial intelligence
(AI) based parameter estimation method called Variance Preserving Flow
Matching Posterior Estimation (VPFMPE). Specifically, we utilize triangular
interpolation to maintain variance over time, thereby constructing a transport
path for training continuous normalization flows. Compared to the simple
linear interpolation method used in flow matching to construct the optimal
transport path, our approach better captures continuous temporal variations,
making it more suitable for the parameter estimation of MBHBs. We also
constructed a novel network architecture for gravitational wave feature
extraction that outperforms existing methods in terms of training speed and
accuracy. Additionally, we creatively introduce a parameter transformation
method based on the symmetry in the detector’s response function. This
transformation is integrated within VPFMPE, allowing us to train the model
using a simplified dataset, and then perform parameter estimation on more
general data, hence also acting as a crucial factor in improving the training
speed. In conclusion, for the first time, within a comprehensive and
reasonable parameter range, we have achieved a complete and unbiased
11-dimensional rapid inference for MBHBs in the presence of astrophysical
confusion noise using ODE-based generative models. In the experiments based on
simulated data, our model produces posterior distributions comparable to those
obtained by nested sampling.

##  1 Introduction

Gravitational wave (GW) astronomy has made significant progress due to the
breakthroughs in ground-based detections led by the LIGO-Virgo-KAGRA network
[1, 2, 3]. Scheduled in the upcoming decade, space-based detectors such as the
Laser Interferometer Space Antenna (LISA) [4, 5], Taiji [6, 7, 8], and Tianqin
[9, 10] aim to detect GWs in the 0.1 mHz - 1 Hz frequency band associated with
enormous astrophysical and cosmological sources [11]. As one of the
representative projects, Taiji consists of a triangle of three spacecraft
(S/Cs) with a baseline separation of 3 million kilometers and aims to detect
low-frequency GWs emitted by sources such as compact galactic binaries (GBs)
[12], massive black hole binaries (MBHBs) [13], extreme mass ratio inspirals
(EMRIs) [14], as well as the stochastic gravitational wave background of
astrophysical or cosmological origins.

The merger of MBHB with component masses $10^{4}\sim 10^{7}M_{\odot}$ is one
of the primary observation targets of space-based GW detectors, since they may
shed light upon the growth and merger history of massive black holes, the
dynamic behavior of curved spacetime, and the nature of gravity, etc. Due to
the presence of the surrounding gas environment, MBHBs may produce detectable
electromagnetic (EM) signals in addition to GWs before, during, and after the
mergers [15, 16, 17, 18]. Joint EM and GW observation of the same source may
unveil the mystery of the evolution of MBHBs and their host galaxies, and
provide a new measurement of the cosmic expansion rate, hence shedding light
upon the Hubble tension problem in cosmology [19, 20]. For these purposes, it
is crucial that the estimation of MBHB parameters, especially the sky location
(ecliptic longitude, latitude), should be conducted with high precision and
low latency, to guide the search of EM observatories.

Besides, the parameter estimation of merging MBHB must also be robust against
the contamination of the confusion foreground formed by tens of millions of
Galactic binaries. Indeed, in the area of space-based GW detection, data
processing faces the challenge caused by the presence of a large number of
astrophysical sources in the measurement band [21]. The non-trivial overlap of
individual signals in the time and frequency domains increases the complexity
of the task. The accurate analysis of up to $\mathcal{O}(10^{4})$ overlapping
sources in the data requires the use of global fit techniques [22, 23, 24],
namely jointly fitting all sources and noise in the data, which imposes a high
challenge on both computational time and resources. The prototype global fit
pipeline proposed by Ref. [24] currently requires as many as
$\mathcal{O}(10^{3})$ CPUs and several days to process the 12-month Lisa Data
Challenge (LDC) data, and the issue of computational resources would become
even more critical for the multi-year datasets in the future. Due to these
challenges and problems, the aim of this paper is to develop an accelerated
and robust method to estimate the parameters of merging MBHBs.

When discussing the efficiency and speed of parameter estimation, the
advancement in machine learning technology is revolutionizing our
understanding of this area. Currently, most machine learning methods applied
to the merging of MBHB use neural posterior estimation(NPE) [25, 26, 27] of
discrete normalized flows. NPE has been widely used in fields such as
atmospheric retrievals (AR) [28] and ground-based gravitational wave
parameters estimation [29, 30]. NPE was also used by the winning entry to the
2023 edition of the ARIEL data challenge [31].

The mission lifetimes of LISA, Taiji, and Tianqin are planned to last for
years, during which the payloads, orbits, and noise characteristics will vary
over time. Accordingly, the models used for parameter estimation may also need
to be regularly updated. As a result, the training speed of the model will be
linked to the timeliness of scientific output. However, a significant
challenge to the current NPE model in MBHB inference is the extensive training
duration, which hampers its ability to conduct real-time data analysis.
Additionally, previous research [32] in this area shows that NPE falls short
in accurately reproducing the posterior distribution obtained by the benchmark
Nested Sampling method. This prolonged training period not only delays the
application of the model but also limits its scientific utility in rapidly
evolving research contexts. To address these challenges, we applied the latest
generative model to MBHBs and made further improvements and optimizations to
better solve the problem. Compared to NPE, Flow Matching Posterior Estimation
(FMPE) offers more flexibility in designing non-diffusion paths, such as
optimal transport paths, and enables direct access to density functions [33].
Moreover, this method can achieve faster and more accurate results than NPE.
Nevertheless, FMPE relies solely on simple linear interpolants when
constructing optimal transport paths. This reliance impacts the extraction and
modeling of MBHB waveform features, which in turn affects the efficiency and
accuracy of parameter estimation.

Therefore, to capture the continuous temporal variations and enhance the
parameter estimation of MBHBs, we aim to explore more effective methods for
constructing transport paths and MBHB wave source feature extraction, which
are particularly beneficial and targeted for extracting the continuous
waveform characteristics of MBHBs. We consider utilizing trig interpolants to
design transport paths [34, 35] that allow for the continuous evolution of
MBHB waveform features along the time axis for modeling. This modeling
accurately bridges any two probability density functions within a finite time,
which belongs to a broad category of continuous-time stochastic processes. We
call this method Variance Preserving Flow Matching Posterior Estimation
(VPFMPE). Compared to standard Monte Carlo methods, the VPFMPE model achieves
comparable performance.

Figure 1:  The Diagrams for the training and sampling phases of our model. On
the lower left, the training phase is displayed, where $t$ denotes time and
$x_{1}$ represents the data from the first day of observed MBHBs. In this
phase, a feature extractor, as detailed in section 3.1, processes the data to
produce a final vector field. This vector field enables the network
architecture to learn a transformation from a base distribution ($\theta_{0}$,
sampled from a Gaussian distribution) to a posterior distribution
($\theta_{1}$, sampled from the posterior distribution). On the lower right,
the sampling phase is shown. Here, the data with noise ($x_{t}$), which
includes both confusion noise and instrumental noise from MBHB observations at
any time, is mapped back to the first day’s data $x_{1}$ using the coordinate
transformation described in section 3.2. The Gaussian distribution is then
converted into a posterior distribution by solving the ordinary differential
equation (ODE) with the vector field learned during the training phase. The
top of the image displays simulation data from both the training and sampling
phases.

The challenges of parameter estimation for MBHBs not only involve the need for
sophisticated algorithms capable of advanced feature extraction and inference,
but also lie in the excessively wide prior range of parameters and the
resulting complexity of data. Recent studies [36, 32] adopted priors that do
not fully cover the complete physical ranges of parameters. This compromise
was due to that setting broad prior ranges during model training can reduce
the effectiveness of the model. This underscores the importance of finding a
careful balance between the breadth of prior ranges and the accuracy of
estimation. To date, no research has managed to achieve unbiased parameter
estimation across the entire range of physical priors. Our method is proven to
overcome this limitation. In this paper, we expanded the prior ranges of
parameters to improve the applicability model.

In order to further enhance the efficiency and effectiveness of model
training, we tackled a common complexity faced by space-based gravitational
wave detectors, namely the time-varying response function. Utilizing the
coordinate transformation of parameters, we can essentially conduct both the
model training and the parameter estimation in a coordinate frame fixed to the
detector, in which the response function is time-invariant. This strategy
significantly simplifies the training task. Consequently, the integration of
coordinate transformation and advanced generative models enhances the speed
and accuracy of parameter estimation for MBHBs across an extensive range of
priors.

In addition, we have developed a novel network architecture for gravitational
wave feature extraction. This architecture illustrated in Figure 4, named
TResMlp, introduces the innovation of encoding time $t$, resulting in superior
training speed and accuracy compared to existing methods [33, 32, 36].

The training of the model can be finished in about two days, as shown in
Figure 1, which is roughly three times faster than the training duration in
our previous work [32]. Additionally, its performance on space-based
gravitational wave data is on par with that of traditional Monte Carlo
methods. This demonstrates the significant potential of integrating artificial
intelligence with VPFMPE techniques to achieve results that are comparable to
those obtained with standard Monte Carlo methods.

##  2 Data Generation and Preprocessing

In this section, we briefly describe the simulation of data used in model
training and validation. The IMRPhenomD model is adopted to describe the
source-frame waveforms of MBHBs. Given the widely recognized necessity of
using Time Delay Interferometry (TDI) to suppress laser frequency noise,
Taiji’s response to the MBHB signals is calculated in the form of two noise-
orthogonal TDI channels $(A,E)$. Specifically, the fast Fourier domain
transfer function derived by Ref. [37] is employed to project the IMRPhenomD
waveform to the $(A,E)$ channels, which can be written in a concise form as

$\tilde{A}(f)=\sum_{\alpha}\mathcal{T}^{A}_{\alpha}(f,t_{f})\tilde{h}_{\alpha}(f),\quad\tilde{E}(f)=\sum_{\alpha}\mathcal{T}^{E}_{\alpha}(f,t_{f})\tilde{h}_{\alpha}(f),\quad\alpha\in\\{+,\times\\}.$
(1)

As can be seen, the transfer function has a close dependence on the orbits of
the detector (through the $t_{f}$ term). For the sake of speed and simplicity,
the motions of Taiji’s satellites are modeled based on an equal-arm analytic
orbit with nominal arm-length $L=3\times 10^{9}{\rm m}$ and an orbital period
of 1 year.

The TDI response of each MBHB signal can be fully specified by 11 parameters:
the chirp mass $\mathcal{M}_{c}$, mass ratio $q\equiv m_{2}/m_{1}$, the
dimensionless spins along the $z-$axis $\chi_{z1}$ and $\chi_{z2}$, the time
$t_{c}$ and phase $\varphi_{c}$ at coalescence, the luminosity distance
$d_{L}$, the inclination angle $\iota$, polarization angle $\psi$, and the
ecliptic coordinate $\\{\lambda,\beta\\}$, being the longitude and latitude
respectively. The priors on these parameters are listed in Table 1. Similar to
the LDC-sangria data, the prior range of $\mathcal{M}_{c}$ covers an order of
magnitude. In our previous work [32], considering the complexity of data, we
limited the prior ranges of spins to positive values. On this simplified
parameter estimation task, NPE has been demonstrated to achieve unbiased
estimations. However, its performance is subpar when applied to MBHB data with
a complete spin prior. An advantage of the models proposed in this paper is
that they are applicable to the complete physical prior range of spins (i.e.
$[-1,1]$). Parameter $t_{c}$ is defined relative to a reference time $t_{\rm
ref}$. Since the mergers of MBHBs manifest as prominent peaks in the time
domain, we consider a prior range of $\sim 2000$ s ( $\sim 0.02$ day) around
$t_{\rm ref}$. Notably, to simplify the dataset and reduce the difficulty of
training, $t_{\rm ref}$ is fixed to the $1^{\rm th}$ day in the training set.
The response function varies with the motion of detectors with a 1-year
period, thus a training set with $t_{c}\in 1\pm 0.01$ day is far from being
sufficient to cover all the observable MBHBs. While, thanks to the method
described in Sec. 3.2, our model can be generalized to make inferences
throughout the entire task cycle. Besides, the range of $d_{L}$ corresponds to
the redshift range of $z\in[1,10]$, which is sufficient to encompass the
majority of MBHBs observable by space-based GW detectors. The angles $\iota$
and $\beta$ are uniformly distributed in terms of $\cos\iota$ and $\sin\beta$
to represent that the sources are isotropically distributed on the sphere.

Table 1: The prior ranges of parameters. The last column is the alias of
parameters used in the code as well as the corner plot.

Parameter Description Prior Lower Bound Prior Upper Bound Alias in Code
$\mathcal{M}_{c}$ Chirp mass of the binary system $1.0\times 10^{5}M_{\odot}$
$1.0\times 10^{6}M_{\odot}$ chirp_mass $q$ Ratio of masses of the two black
holes 0.1 1.0 mass_ratio $\chi_{z1}$ Spin of the 1st black hole along the
$z$-axis -1 1 spin_1 $\chi_{z2}$ Spin of the 2nd black hole along the $z$-axis
-1 1 spin_2 $t_{c}$ -0.01 Time of coalescence relative to the reference time
111The reference time $t_{\rm ref}$ is set to the $1^{\rm th}$ day., $t_{\rm
ref}$ day 0.01 day coalescence_time $\varphi_{c}$ Phase at the moment of
coalescence 0 rad $2\pi$ rad coalescence_phase $d_{L}$ Luminosity distance to
the binary system 6000 Mpc 100000 Mpc luminosity_distance $\iota$ Angle of
inclination of the binary orbit 0 rad $\pi$ rad inclination $\beta$ Ecliptic
latitude of the binary system $-\frac{\pi}{2}$ rad $\frac{\pi}{2}$ rad
latitude $\lambda$ Ecliptic longitude of the binary system 0 rad $2\pi$ rad
longitude $\psi$ Polarization angle of the gravitational wave 0 rad $\pi$ rad
psi

The majority of MBHB’s signal-to-noise ratio(SNR) is contributed by a short
duration including the merger phase (typically several hours) [38]. To
reasonably avoid the potential overlap of other signals, the time duration for
each sample is set to one day, with the merger positioned towards the end of
this window. Furthermore, the sampling frequency is set as $1/15$ Hz,
indicating a Nyquist frequency of $0.033$ Hz. This frequency exceeds the
highest instantaneous frequency of the waveforms, ensuring comprehensive
capture of the signal.

Figure 2: A representative data in the test set. The orange and green curves
represent the Fourier transforms of the TDI-A, E data streams, which include
GW signals, instrumental noises, and confusion noise. To clearly illustrate
the relative strengths of the signal and noises, the blue curve depicts the
Fourier transform of the confusion noise, while the gray dashed curve
indicates the theoretical amplitude of the instrumental noises.

Along with the detector responses, we also need the PSD of the noise in each
channel. For instrumental noises, according to the current design of the Taiji
mission, the one-sided ASDs of the two top-level instrumental noises, namely
optical measurement system noise and test mass acceleration noise, are  [39]

$\displaystyle\sqrt{S_{\mathrm{OMS}}(f)}=$ $\displaystyle 8\times
10^{-12}\sqrt{1+\left(\frac{2\mathrm{mHz}}{f}\right)^{4}}\frac{{\mathrm{m}}}{\sqrt{\mathrm{Hz}}},$
(2) $\displaystyle\sqrt{S_{\mathrm{ACC}}(f)}=$ $\displaystyle 3\times
10^{-15}\sqrt{1+\left(\frac{0.4\mathrm{mHz}}{f}\right)^{2}}\sqrt{1+\left(\frac{f}{8\mathrm{mHz}}\right)^{4}}\frac{{\mathrm{m/s^{2}}}}{\sqrt{\mathrm{Hz}}},$
(3)

The total instrumental noise PSDs of the $A,E$ channels $S_{n}^{A,E}(f)$ are
the combinations of these two components (see Ref. [40, 41, 42] for their
specific forms). For simplicity, we treat $S_{n}(f)$ as constant over the
whole mission duration. Our training data is a combination of “clean” TDI
responses and stationary instrumental noises. On a more realistic ground, the
long-term variations and short-term jitters of instrumental noise, along with
the typical anomalies such as glitches and gaps should be considered in future
work.

To test the robustness and generalizability of the model, during the inference
stage, we consider a non-stationary noise caused by the unresolvable
foreground originating from the overlapping GWs of GBs. This confusion noise
is generated with the FastGB algorithm [43] from a catalog of $\sim 3\times
10^{7}$ GBs provided by the LDC dataset, from which $\sim 2\times 10^{4}$
bright sources with SNR greater than 7 were subtracted. This choice of SNR
threshold follows the convention of [44, 45, 46, 47]. Shown in Figure 2 is
representative data in the test set, including GW signal, instrumental noises,
and confusion noise.

##  3 Methodological Framework

The application of machine learning techniques in the field of gravitational
waves has seen widespread success and remarkable achievements [48, 49]. Among
these techniques, deep learning has demonstrated exceptional prowess,
especially in processing ground-based gravitational wave data [50, 51, 52, 53,
54, 55, 56, 57, 58].

NPE with normalizing flows. Neural Posterior Estimation (NPE) combined with
Normalizing Flows (NF) is used to directly fit the posterior distribution
$p(\theta|x)$. In this approach, an invertible mapping $\psi_{x}$ transforms
the base distribution $p_{0}(\mathbf{z})$ into the target distribution
$q(\theta|x)$:

$q(\theta|x)=p_{0}(\psi_{x}^{-1}(\theta))\det\left|\left(\frac{\partial\psi_{x}^{-1}(\theta)}{\partial\theta}\right)\right|.$
(4)

NPE is trained using a maximum likelihood objective and simplifies the
expectation using Bayes’ theorem.

$\displaystyle L_{\mathrm{NPE}}=-E_{p(\theta)p(x|\theta)}\log q(\theta|x).$
(5)

Once training is complete, NPE can efficiently perform inference on new
observations $x$ and provide exact density evaluations of $q(\theta|x)$. This
method retains the expressive power and computational efficiency of
normalizing flows while achieving posterior fitting.

Continuous Normalizing Flows for parameter estimation. Continuous Normalizing
Flows (CNFs) transform a simple base distribution into a more complex
distribution by describing this transformation through continuous time
$t\in[0,1]$. For each $t$, the flow is defined by a vector field $v_{t,x}$ on
the sample space. This vector field represents the trajectory velocity of the
samples at any given time $t$.

$\frac{d}{dt}\psi_{t,x}(\theta)=v_{t,x}(\psi_{t,x}(\theta)).$ (6)

The conversion between the base distribution (at $t=0$) and the target
distribution (at $t=1$) is achieved by integration. A straightforward approach
to training Continuous Normalizing Flows (CNFs) involves starting with a given
base distribution and using an ODE solver to obtain the target distribution.
This target distribution is then constrained to match the real data
distribution by minimizing a divergence measure, such as the KL divergence.

Because many intermediate trajectories are unknown, inferring the distribution
(either through sampling or likelihood computation) requires repeatedly
simulating the ODE, leading to a significant computational burden.

However, Flow matching provides alternative training objectives for CNFs, Dax
et al. [33] apply flow matching [59] to simulation-based inference. The loss
function of FMPE is expressed as:

$L_{FMPE}=E_{t\sim p(t),\theta_{1}\sim p(\theta),x\sim
p(x|\theta_{1}),\theta_{t}\sim p_{t}(\theta_{t}|\theta_{1})}\parallel
v_{t,x}(\theta_{t})-u_{t}(\theta_{t}|\theta_{1})\parallel^{2}.\ \ $ (7)

Here, $v_{t,x}(\theta_{t})$ represents the vector field that generates the
path to the target probability distribution, and
$u_{t}(\theta_{t}|\theta_{1})$ is the sample-conditioned vector field. FMPE
introduces a Gaussian path with time-dependent means $t\theta_{1}$ and
standard deviations $1-(1-\sigma)t$.

$p_{t}(\theta|\theta_{1})=N(t\theta_{1},1-(1-\sigma)t).$ (8)

FMPE defines the sample-conditioned vector field is

$u_{t}(\theta|\theta_{1})=\frac{\theta_{1}-(1-\sigma)\theta_{0}}{1-(1-\sigma)t}.$
(9)

In our work, we reference the studies by Albergo et al. [34, 35] and utilize
triangular interpolation to define the time-dependent mean and standard
deviation, thus constructing the Gaussian path and sample-conditioned vector
field as follows:

$\displaystyle
p_{t}(\theta|\theta_{1})=N(\sin(\frac{\pi}{2t})\theta_{1}+\cos(\frac{\pi}{2t})\theta_{0},\sigma).$
(10)

$\displaystyle
u_{t}(\theta|\theta_{1})=\frac{\pi}{2(\cos(\frac{\pi}{2t})\theta_{1}-\sin(\frac{\pi}{2t})\theta_{0})}.$
(11)

We call this method Variance Preserving Flow Matching Posterior Estimation and
use it to train continuous normalization flows. We show the differences
between the two methods in Table 2.

Table 2: A table with constructing transport paths differences. FMPE uses linear interpolation to construct transport paths, while VPFMPE uses trig interpolants to construct transport paths. | Interpolar Method | Mean of the probability path | Standard deviation of the probability path |
---|---|---|---|---
FMPE | linear | $t\theta_{1}$ | $1-(1-\sigma)t$ |
VPFMPE | trig | $\sin(\frac{\pi}{2t})\theta_{1}+\cos(\frac{\pi}{2t})\theta_{0}$ | $\sigma$ |
Figure 3: Comparison of the posterior distributions of the two methods. Blue
represents the posterior distribution of the VPFMPE method on the injected
MBHB signal. Orange represents the posterior distribution of the FMPE method
on the injected MBHB signal. The red line represents the true value of the
injected MBHB signal.

Comparison of methods. To demonstrate the performance of the transport path
constructed by VPFMPE, we randomly injected an MBHB signal and compared the
FMPE method with the VPFMPE method. Both models were trained using the same
learning rate, optimizer, and gravitational wave feature extractor to ensure
fairness. We show in Figure 3 a comparison of the posterior estimates of the
four parameters of the injected MBHB signal by the two methods.

Discussion. The NPE network just only one-third allocated to the feature
extractor network for GW data. However, the ODE-based generative model only
needs to parameterize a simple vector field, so it can devote its full network
capacity to interpreting the higher dimensions [33].

###  3.1 Network architecture

A key challenge in space-based gravitational wave parameter estimation lies in
managing the high-dimensional nature of its data. To address this, it is
essential to develop a feature extraction network capable of transforming the
high-dimensional input data into a more manageable, low-dimensional
representation [29, 32, 60]. This transformation is critical for adjusting the
dimensionality of the final feature vector to match that of the parameter
vector $\theta$, enabling the network to effectively parameterize the
conditional vector field $(t,x,\theta)\rightarrow v_{t,x}(\theta)$ [33].
Achieving accurate parameterization of this conditional vector field is
fundamental to the construction of an effective network structure.

Figure 4:  TResMlp Network Architecture: Initially, the frequency vector $f$
is used to encode the position of time $t$. Sine and cosine functions are then
applied to extract periodic features from time $t$, resulting in the “Time
Embedding”. Subsequently, a “Linear” layer maps the two channels of the space-
based GW data to produce the “Data Embedding”. Finally, all obtained
embeddings are input into a Multi-Layer Perceptron (MLP) with a residual
structure to generate the final output $v_{(t,x_{1})}(\theta_{t})$.

In recent years, DeResNet networks have been commonly used as the backbone
network for feature extraction in ground-based gravitational wave parameter
estimations or space-based gravitational wave parameter estimation [57, 32,
33, 53], and Gated Linear Unit(GLU)-conditioning [61] can be optionally used
for information aggregation of $(\theta,t)$.

However, if information about the actual parameter $\theta$ is not properly
handled or retained during the process of aggregating $(\theta,t)$ information
using GLU, it could potentially lead to a loss of information about $\theta$.
To encode the time $t$ more effectively and ensure that our feature extraction
network captures its dynamic nature, we introduce a frequency vector $f$. This
frequency vector $f$ adapts dynamically to changes in the parameters of the
wave source, allowing for a more precise encoding of time $t$. To capture the
periodic features of time, we employ sine and cosine functions for encoding.
Our network architecture, referred to as TResMlp, is shown in Figure 4.

We trained the CNFs of MBHBs using two different network architectures. To
ensure fairness, the two models had identical parameters such as the optimizer
and learning rate, differing only in their network structures. Additionally,
both architectures used the same hidden layer configuration. We show in Figure
5 a comparison of the posterior estimates of the four parameters of the
injected MBHBs signal by the two models.

Figure 5: Comparison of the posterior distributions of the two network
architectures. The blue represents the posterior distribution plot of CNFs
trained using the TResMlp network. The orange represents the posterior
distribution plot of CNFs trained using the DeResNet network. The red line
represents the true value of the MBHBs.

Compared to the DeResNet model that use GLU-conditioning for aggregating
information about $(\theta,t)$ [57, 32, 33, 53], our network architecture,
TResMlp, demonstrates faster training speed and higher precision with the same
hidden layers. More experiments and comparisons of related models are
described in detail in Section 4.

###  3.2 Coordinate Transformation and Model Generalization Strategy

To simplify the dataset and reduce the difficulty of training, the coalescence
times of MBHBs are set within $1\pm 0.02$ day in the training set (see Sec.
2). While the response function of detector to GW signals varies in time with
a one-year period, thus in principle the model trained on this dataset can
only be applied to MBHBs that merge within $1\pm 0.02$ day. While, thanks to
the generalization strategy in this section, the application scope of the
model can be extended to the whole mission lifetime.

In this paper, the time of coalescence (dubbed $t_{c}$) is defined relative to
a reference time $t_{\rm ref}$. The complexity of space-based GW data partly
lies in the annual modulation of detector response due to the orbital motion
of detector constellation around the sun. In other words, the mapping between
GW data and MBHB parameters varies with $t_{\rm ref}$. This significantly
increases the complexity of the training set, which would severely damage the
training efficiency and the performance of the model. To tackle this
challenge, we exploit the symmetry (invariance) of the system, fix the
training data around a certain $t_{\rm ref0}$, and then the trained model is
generalized via our transformation rule to make inferences at any time. The
idea behind this is a paradigm that combines physical information with neural
networks, and methods belonging to this family, such as the group-equivariant
neural posterior estimation [62, 60], have shown broad application prospects.

Figure 6: Comparison of the effect of the model generalization strategy on
model performance. The left panel shows the P-P plot without the model
generalization strategy, while the right panel shows the P-P plot with the
model generalization strategy. Both plots are for a set of 4000 injections
with added instrument noise and reference times ranging from 1 to 365 days.
Ideally, the cumulative distribution function of these percentiles for each
parameter should align closely with the diagonal line, indicating optimal
network performance. In such a case, the percentiles would be uniformly
distributed between 0 and 1. The grey regions on the plot represent the
$1\sigma$ and $2\sigma$ confidence intervals.

We describe the transformation rule as follows. We define two coordinate
systems, namely the Solar System Barycenter (SSB), in which physical
quantities are indexed with “S”, and the Taiji detector reference frame, in
which physical quantities are indexed with “T”. The parameters which vary with
the coordinate systems are $\\{t_{c},\lambda,\beta,\psi\\}$, where $t_{c}$
should be interpreted as the time of arrival at the origin of the coordinate
system (relative to $t_{\rm ref}$). Since the data in consideration last only
for one day, for preliminary parameter estimation tasks, the location of MBHB
remains stationary relative to the Taiji reference frame, and the mapping
between source parameters $\bm{\theta}_{T}$ and TDI response is independent of
$t_{\rm ref}$. This is precisely the symmetry based on which we can achieve
model generalization. The $t_{\rm ref}$-dependent coordinate transformation
from Taiji frame to SSB frame is denoted as $f_{t_{\rm
ref}}:\bm{\theta}_{T}\rightarrow\bm{\theta}_{S,{t_{\rm ref}}}$ (see Ref. [63]
for the expression of $f$). The parameters are firstly inferred with the model
trained around $t_{\rm ref0}$, hence dubbed $\bm{\theta}_{S,t_{\rm ref0}}$,
and the targeted parameters $\bm{\theta}_{S,t_{\rm ref}}$ can be obtained
following two steps:

$\bm{\theta}_{T}=f_{t_{\rm ref0}}^{-1}\left(\bm{\theta}_{S,t_{\rm
ref0}}\right),\quad\bm{\theta}_{S,t_{\rm ref}}=f_{t_{\rm
ref}}\left(\bm{\theta}_{T}\right),$ (12)

or in terms of posterior distributions:

$p\left(\bm{\theta}_{S,t_{\rm ref}}|d\right)=p\left(\bm{\theta}_{S,t_{\rm
ref0}}|d\right)J\left(f^{\ }_{t_{\rm ref0}}\right)J\left(f^{-1}_{t_{\rm
ref}}\right),$ (13)

with $J\left(f^{\ }_{t_{\rm ref0}}\right)$ and $J\left(f^{-1}_{t_{\rm
ref}}\right)$ being the Jacobian determinants corresponding to $f_{t_{\rm
ref0}}$ and $f^{-1}_{t_{\rm ref}}$, respectively. Formula (13) looks similar
to the transformations in normalizing flows, except that the mapping
$f_{t_{\rm ref}}$ is constructed based on our understanding of the system,
instead of training a neural network. Training a network to represent this
transformation is unnecessary since we possess comprehensive knowledge of its
form and can effortlessly compute it. As a result, The embedding of this
physical information greatly simplifies training tasks.

Below, we describe how models can be combined with the model generalization
strategy. Firstly, during the dataset generation, the merger times of MBHBs
are limited to 0.02 days centered at a specific reference time (the 1st day).
This significantly simplifies the dataset’s complexity and reduces the
challenge of MBHB parameter estimation. Subsequently, we employ VPFMPE as the
loss function for training CNFs. We utilize the novel network architecture
TResMlp as the gravitational wave feature extractor, labeling this approach as
VPFMPE-TResMlp.

In the inference phase, we use the model generalization strategy to map the
reference time data for inference to the reference time data fixed during
training. This allows the model to be trained on only one day of data but
tested and applied to data spanning an entire year. We use the P-P plot in
Figure 6 to demonstrate the importance of using a model generalization
strategy.

##  4 Experiments

To ensure unbiased parameter estimation, we examine our algorithm on injected
GW signals and noise in this section. For the parameter estimation of the
signals, we employed the four models mentioned above: VPFMPE-TResMlp, FMPE-
TResMlp, FMPE-DeResNet [33], and NPE-DeResNet [32]. We applied the model
generalization strategy 3.2 to all four models by default. Otherwise, due to
the complexity of the dataset, the models would lack inference capability.

Table 3: In the 850 simulations, we compared the correlation coefficients between the true values and the recovered values for the four models. In addition, in the last column of the table, we have included the training time for the three models. | $\mathcal{M}_{c}$ | $q$ | $\chi_{z1}$ | $\chi_{z2}$ | $d_{L}$ | Training time
---|---|---|---|---|---|---
VPFMPE-TResMlp | 0.9978 | 0.9790 | 0.9723 | 0.8163 | 0.9408 |  $\approx$ 2.0 days
FMPE-TResMlp | 0.9978 | 0.9787 | 0.9719 | 0.8139 | 0.9400 |  $\approx$ 2.0 days
FMPE-DeResNet | 0.9970 | 0.9762 | 0.9563 | 0.7174 | 0.9397 |  $\approx$ 2.3 days
NPE-DeResNet | 0.9261 | 0.8647 | 0.6620 | 0.1681 | 0.8634 |  $\approx$ 6.5 days

To demonstrate the effectiveness of our algorithm and validate that the
recovered probability distribution accurately represents the confidence we
should have in the signal parameters, we conducted extensive testing on
simulated signals and noise. We injected 850 waveforms drawn from the prior,
along with residual instrumental and foreground confusion noise. We will take
the median of the recovered values as the point estimate and calculate the
correlation coefficient with the true values of the parameters. The
correlation coefficient reflects the degree of linear correlation between the
true value distribution and the estimated values. A coefficient value of 1
indicates perfect consistency between the estimated and true values. Note that
the reference times of these testing signals are randomly distributed
throughout the year, rather than being fixed to the 1st day.

Our results are presented in Table 3. In the final column of the Table 3, we
have included the approximate training times for the three models. The VPFMPE-
TResMlp model achieved the best results and had the fastest training time. For
$t_{c}$, $\varphi_{c}$, $\iota$, $\beta$, $\lambda$, and $\psi$, due to the
presence of multimodality (i.e. there is more than one peak in the posterior
distributions), using the median as a point estimate is not applicable.

Figure 7:  Corner plot of parameter posterior distributions for gravitational
wave signals injected with instrumental and confusion noise. The blue contours
illustrate the two-dimensional joint posteriors obtained from our (VPFMPE-
TResMlp) model, with peaks selected based on $t_{c}$ that align with true
values. In contrast, the orange contours represent the corresponding
posteriors obtained from benchmark analysis using Nested Sampling. For each
model, the contour boundaries respectively encompass 95% and 60% confidence
levels. The true parameter values of the simulated signal are indicated by red
vertical and horizontal lines.

To date, existing research has not yet succeeded in providing unbiased
estimates of MBHBs across the entire parameter space. For example, Ref. [36]
only focused on four parameters, and Ref. [32] have limited to positive spins.
Also, we use machine learning techniques for complete parameter estimation of
MBHBs and achieve comparable results to nested sampling. Therefore, in
subsequent comparisons, we will only compare the results of VPFMPE-TResMlp to
those of the MCMC method.

The aim of our ODE-based generative models for MBHBs is to reproduce the
comparable posterior distribution as the MCMC method. Thus, by comparing the
results of Bayesian inference on the same event, the performance of the model
can be better assessed. In this case, we use Nested Sampling [64, 65, 66] as a
Monte Carlo technique for comparison experiments, as shown in Figure 7.

An illustration demonstrating an example using the same prior probabilities
outlined in Table 1 is provided in Figure 7, where a test injection is
subjected to confusion noise on the 30th reference day. The injected GW signal
originates from an MBHB with a total signal-to-noise ratio of 1543, and its
parameters are detailed in Table 4. Compared to the standard nested sampling
technique, our model exhibited results similar to the traditional nested
sampling, demonstrating consistency in the sense of the degeneracy among
extrinsic parameters 4. This consistency aligns with the crucial physical
characteristics of MBHBs. Note that compared to nested sampling, our machine
learning process always produces similar posterior distributions on $t_{c}$,
$d_{L}$, $\iota$, $\beta$, $\lambda$, and $\psi$ 7.

We have achieved accurate estimates of the sky position parameters. This is
most advantageous since accurate and reliable estimates of the sky position
parameters (longitude, latitude) are useful in guiding the search for
electromagnetic observatories.

We have collated the true values of the injections alongside the results
generated by both the model and the nested sampling process in Table 4. This
provides a concise overview of the performance of each method.

Table 4: This table compares the parameters injected in the prior with the parameters recovered by the Nested Sampling and VPFMPE-TResMlp methods. The GW signal is emitted by a merging MBHB with a total SNR of 1543. The recovered values are accompanied by their 1 $\sigma$ confidence regions. Parameter | Injected value | VPFMPE-TResMlp | Nest
---|---|---|---
$\mathcal{M}_{c}$ [$M_{\odot}$] | 531067.4 | $528444.8000_{-1206.2360}^{+1141.3900}$ | $531038.3974_{-107.2206}^{+125.0852}$
$q$ | 0.6355 | $0.5994_{-0.0220}^{+0.0236}$ | $0.6405_{-0.0070}^{+0.0060}$
$\chi_{z1}$ | 0.2930 | $0.4051_{-0.1042}^{+0.0858}$ | $0.2636_{-0.0259}^{+0.0318}$
$\chi_{z2}$ | 0.2726 | $0.0273_{-0.1552}^{+0.1807}$ | $0.3233_{-0.0517}^{+0.0409}$
$t_{c}$ [day] | 0.00028 | $0.0004_{-0.0001}^{+0.0001}$ | $0.0003_{-0.0000}^{+0.0000}$
$\varphi_{c}$ [rad] | 0.5412 | $3.1765_{-1.8449}^{+1.7528}$ | $3.6714_{-3.0978}^{+1.6062}$
$d_{L}$ [Mpc] | 47054.1 | $48698.8670_{-1486.6328}^{+1568.7870}$ | $46334.8753_{-891.2652}^{+882.4137}$
$\iota$ [rad] | 0.8207 | $0.8004_{-0.0207}^{+0.0191}$ | $0.8298_{-0.0109}^{+0.0098}$
$\beta$ [rad] | 0.0459 | $0.1405_{-0.0317}^{+0.0340}$ | $0.0401_{-0.0188}^{+0.0188}$
$\lambda$ [rad] | 3.3544 | $3.3024_{-0.0322}^{+0.0317}$ | $3.3350_{-0.0167}^{+0.0188}$
$\psi$ [rad] | 1.5634 | $1.5680_{-1.5388}^{+0.0442}$ | $1.5692_{-1.5613}^{+0.0188}$

A noteworthy phenomenon is that our machine learning pipeline produces a
posterior similar to the nested sampling approach in terms of extrinsic
parameters, and the posterior distributions of both converge in the same order
of magnitude.

In the final part of the results, we explore the computational efficiency of
the algorithm. Table 5 shows the speed of generating posterior samples using
nested sampling compared to the VPFMPE method. It is important to note that we
used nested sampling to search near the true intrinsic parameter values of the
injected MBHB. In practice, the time consumed by nested sampling is
significantly longer than the times listed in the table. Even under these
conditions, our VPFMPE method remains two orders of magnitude faster than the
nested sampling method.

Table 5: Comparison of computational time required by different posterior sampling approaches to generate their respective samples. | Number of Posterior Samples | Runtime (seconds) | Time per Sample |
---|---|---|---|---
Nested Sampling | 11215 | 3000 | 0.26750 |
VPFMPE | 14336 | 240 | 0.01674 |

Our model not only prioritizes efficiency but also generates a posterior
distribution similar to nested sampling. Moreover, it is well-suited for the
majority of MBHB sources anticipated by future Taiji observations. This
contribution will significantly enhance global fit analysis.

##  5 Conclusions

In this study, we have demonstrated the inaugural application of ODE-based
generative models in the field of space-based GW detection, taking the
parameter estimation of merging MBHBs as a representative example. We
recommend the use of ODE-based generative models, particularly the VPFMPE
method, as the mainstream machine learning inference approach.

Our method provides comprehensive statistical inference for the 11-dimensional
parameters of MBHBs, establishing a new benchmark in this area. This approach,
as opposed to likelihood-based inference methods, can lead to both faster
training and more accurate estimation.

Furthermore, we propose the integration of model generation strategy based on
the symmetry of the system with ODE-based generative models. This enables us
to simplify the training task, while still obtaining a model that is
applicable to the entire mission period. Note that the coordinate
transformation and ODE-based generative models are theoretically applicable to
other space-based missions such as LISA.

It is also noteworthy that our method exhibits robustness against challenges
such as confusion noise and the complexities of varying reference times,
highlighting its reliability and potential for broader applications. The
significant speed of our method, surpassing traditional techniques by orders
of magnitude, positions it as a valuable tool for preprocessing in global
fitting. Moreover, the precise and reliable estimation of sky location
parameters (longitude, latitude) proves beneficial in guiding the search for
electromagnetic observatories.

This acceleration in data processing and analysis not only streamlines
workflows but also broadens the horizons for gravitational wave data
exploration. Overall, these advancements herald a transformative era in space-
based GW detection, fortifying the capabilities of missions like Taiji and
bringing us closer to unraveling the mysteries of the universe.

##  6 Acknowledgements

This study is supported by the National Key Research and Development Program
of China (Grant No. 2021YFC2201901, Grant No. 2021YFC2203004, Grant No.
2020YFC2200100 and Grant No. 2021YFC2201903). International Partnership
Program of the Chinese Academy of Sciences, Grant No. 025GJHZ2023106GC.

## References

  * [1] The LIGO Scientific Collaboration, J Aasi, et al.  Advanced ligo.  Classical and Quantum Gravity, 32(7):074001, mar 2015.
  * [2] Fet al Acernese, M Agathos, K Agatsuma, Damiano Aisa, N Allemandou, Aea Allocca, J Amarni, Pia Astone, G Balestri, G Ballardin, et al.  Advanced virgo: a second-generation interferometric gravitational wave detector.  Classical and Quantum Gravity, 32(2):024001, 2014.
  * [3] T. Akutsu et al.  Kagra: 2.5 generation interferometric gravitational wave detector.  Nature Astronomy, 3(1):35–40, 2019.
  * [4] Pau Amaro-Seoane, Heather Audley, Stanislav Babak, John Baker, Enrico Barausse, Peter Bender, Emanuele Berti, Pierre Binetruy, Michael Born, Daniele Bortoluzzi, et al.  Laser interferometer space antenna, 2017.
  * [5] John Baker, Jillian Bellovary, Peter L Bender, Emanuele Berti, Robert Caldwell, Jordan Camp, John W Conklin, Neil Cornish, Curt Cutler, Ryan DeRosa, et al.  The laser interferometer space antenna: unveiling the millihertz gravitational wave sky, 2019.
  * [6] Wen-Rui Hu and Yue-Liang Wu.  The Taiji Program in Space for gravitational wave physics and the nature of gravity.  National Science Review, 4(5):685–686, 10 2017.
  * [7] Ziren Luo, Yan Wang, Yueliang Wu, Wenrui Hu, and Gang Jin.  The Taiji program: A concise overview.  Progress of Theoretical and Experimental Physics, 2021(5):05A108, 07 2020.
  * [8] Taiji Scientific Collaboration, Yue-Liang Wu, et al.  China’s first step towards probing the expanding universe and the nature of gravity using a space borne gravitational wave antenna.  Communications Physics, 4(1):34, December 2021.
  * [9] Jun Luo, Li-Sheng Chen, Hui-Zong Duan, Yun-Gui Gong, Shoucun Hu, Jianghui Ji, Qi Liu, Jianwei Mei, Vadim Milyukov, Mikhail Sazhin, et al.  Tianqin: a space-borne gravitational wave detector.  Classical and Quantum Gravity, 33(3):035010, jan 2016.
  * [10] Jun Luo, Li-Sheng Chen, Hui-Zong Duan, Yun-Gui Gong, Shoucun Hu, Jianghui Ji, Qi Liu, Jianwei Mei, Vadim Milyukov, Mikhail Sazhin, Cheng-Gang Shao, Viktor T Toth, Hai-Bo Tu, Yamin Wang, Yan Wang, Hsien-Chi Yeh, Ming-Sheng Zhan, Yonghe Zhang, Vladimir Zharov, and Ze-Bing Zhou.  Tianqin: a space-borne gravitational wave detector.  Classical and Quantum Gravity, 33(3):035010, jan 2016.
  * [11] Gerhard Heinzel, Claus Braxmaier, Roland Schilling, Albrecht Rüdiger, David Robertson, M Te Plate, Vinzenz Wand, K Arai, Ulrich Johann, and Karsten Danzmann.  Interferometry for the lisa technology package (ltp) aboard smart-2.  Classical and Quantum Gravity, 20(10):S153, 2003.
  * [12] Xue-Hao Zhang, Shao-Dong Zhao, Soumya D Mohanty, and Yu-Xiao Liu.  Resolving galactic binaries using a network of space-borne gravitational wave detectors.  Physical Review D, 106(10):102004, 2022.
  * [13] Antoine Klein, Enrico Barausse, Alberto Sesana, Antoine Petiteau, Emanuele Berti, Stanislav Babak, Jonathan Gair, Sofiane Aoudia, Ian Hinder, Frank Ohme, et al.  Science with the space-based interferometer elisa: Supermassive black hole binaries.  Physical Review D, 93(2):024003, 2016.
  * [14] Andrea Maselli, Nicola Franchini, Leonardo Gualtieri, and Thomas P Sotiriou.  Detecting scalar fields with extreme mass ratio inspirals.  Physical Review Letters, 125(14):141101, 2020.
  * [15] Bence Kocsis, Zsolt Frei, Zoltán Haiman, and Kristen Menou.  Finding the electromagnetic counterparts of cosmological standard sirens.  The Astrophysical Journal, 637(1):27, jan 2006.
  * [16] Alessandra De Rosa, Cristian Vignali, Tamara Bogdanović, Pedro R. Capelo, Maria Charisi, Massimo Dotti, Bernd Husemann, Elisabeta Lusso, Lucio Mayer, Zsolt Paragi, Jessie Runnoe, Alberto Sesana, Lisa Steinborn, Stefano Bianchi, Monica Colpi, Luciano del Valle, Sándor Frey, Krisztina É. Gabányi, Margherita Giustini, Matteo Guainazzi, Zoltan Haiman, Noelia Herrera Ruiz, Rubén Herrero-Illana, Kazushi Iwasawa, S. Komossa, Davide Lena, Nora Loiseau, Miguel Perez-Torres, Enrico Piconcelli, and Marta Volonteri.  The quest for dual and binary supermassive black holes: A multi-messenger view.  New Astronomy Reviews, 86:101525, 2019.
  * [17] Milos Milosavljevic and E. S. Phinney.  The Afterglow of massive black hole coalescence.  Astrophys. J. Lett., 622:L93–L96, 2005.
  * [18] M. Dotti, R. Salvaterra, A. Sesana, M. Colpi, and F. Haardt.  On the search of electromagnetic cosmological counterparts to coalescences of massive black hole binaries.  Monthly Notices of the Royal Astronomical Society, 372(2):869–875, 09 2006.
  * [19] Nicola Tamanini, Chiara Caprini, Enrico Barausse, Alberto Sesana, Antoine Klein, and Antoine Petiteau.  Science with the space-based interferometer eLISA. III: Probing the expansion of the Universe using gravitational wave standard sirens.  JCAP, 04:002, 2016.
  * [20] Maria Dainotti, Biagio De Simone, Giovanni Montani;, Tiziano Schiavone;, and Gaetano Lambiase.  The Hubble constant tension: current status and future perspectives through new cosmological probes.  PoS, CORFU2022:235, 2023.
  * [21] Jean-Baptiste Bayle, Béatrice Bonga, Chiara Caprini, Daniela Doneva, Martina Muratore, Antoine Petiteau, Elena Rossi, and Lijing Shao.  Overview and progress on the laser interferometer space antenna mission.  Nature Astronomy, 6(12):1334–1338, 2022.
  * [22] Neil J. Cornish and Jeff Crowder.  Lisa data analysis using markov chain monte carlo methods.  Phys. Rev. D, 72:043005, Aug 2005.
  * [23] Tyson B. Littenberg, Neil J. Cornish, Kristen Lackeos, and Travis Robson.  Global analysis of the gravitational wave signal from galactic binaries.  Phys. Rev. D, 101:123021, Jun 2020.
  * [24] Tyson B. Littenberg and Neil J. Cornish.  Prototype global analysis of lisa data with multiple source types.  Phys. Rev. D, 107:063004, Mar 2023.
  * [25] George Papamakarios and Iain Murray.  Fast $\epsilon$-free inference of simulation models with bayesian conditional density estimation.  In Neural Information Processing Systems, 2016.
  * [26] Jan-Matthis Lueckmann, Pedro J. Gonçalves, Giacomo Bassetto, Kaan Öcal, Marcel Nonnenmacher, and Jakob H. Macke.  Flexible statistical inference for mechanistic models of neural dynamics.  In Neural Information Processing Systems, 2017.
  * [27] David S. Greenberg, Marcel Nonnenmacher, and Jakob H. Macke.  Automatic posterior transformation for likelihood-free inference.  In International Conference on Machine Learning, 2019.
  * [28] Timothy D. Gebhard, Jonas Wildberger, Maximilian Dax, Daniel Angerhausen, Sascha P. Quanz, and Bernhard Scholkopf.  Inferring atmospheric properties of exoplanets with flow matching and neural importance sampling.  ArXiv, abs/2312.08295, 2023.
  * [29] Maximilian Dax, Stephen R. Green, Jonathan R. Gair, Jakob H. Macke, Alessandra Buonanno, and Bernhard Scholkopf.  Real-time gravitational-wave science with neural posterior estimation.  Physical review letters, 127 24:241103, 2021.
  * [30] Stephen R. Green, Christine M. Simpson, and Jonathan R. Gair.  Gravitational-wave parameter estimation with autoregressive neural network flows.  ArXiv, abs/2002.07656, 2020.
  * [31] Mayeul Aubin, Carolina Cuesta-Lázaro, Ethan Tregidga, Javier Viana, Cecilia Garraffo, Iouli E. Gordon, Mercedes López-Morales, Robert J. Hargreaves, Vladimir Yu. Makhnev, Jeremy J. Drake, Douglas P. Finkbeiner, and Phillip A. Cargile.  Simulation-based inference for exoplanet atmospheric retrieval: Insights from winning the ariel data challenge 2023 using normalizing flows.  ArXiv, abs/2309.09337, 2023.
  * [32] Minghui Du, Bo-Hua Liang, He Wang, Peng Xu, Ziren Luo, and Yueliang Wu.  Advancing space-based gravitational wave astronomy: Rapid detection and parameter estimation using normalizing flows.  2023\.
  * [33] Maximilian Dax, Jonas Wildberger, Simon Buchholz, Stephen R. Green, Jakob H. Macke, and Bernhard Scholkopf.  Flow matching for scalable simulation-based inference.  ArXiv, abs/2305.17161, 2023.
  * [34] Michael S Albergo, Nicholas M. Boffi, and Eric Vanden-Eijnden.  Stochastic interpolants: A unifying framework for flows and diffusions.  ArXiv, abs/2303.08797, 2023.
  * [35] Alexander Tong, Kilian Fatras, Nikolay Malkin, Guillaume Huguet, Yanlei Zhang, Jarrid Rector-Brooks, Guy Wolf, and Yoshua Bengio.  Improving and generalizing flow-based generative models with minibatch optimal transport, 2024.
  * [36] Wen-Hong Ruan, He Wang, Chang Liu, and Zong-Kuan Guo.  Rapid search for massive black hole binary coalescences using deep learning.  Physics Letters B, 841:137904, 2023.
  * [37] Sylvain Marsat and John G. Baker.  Fourier-domain modulations and delays of gravitational-wave signals, 6 2018.
  * [38] Geraint Pratten, Antoine Klein, Christopher J. Moore, Hannah Middleton, Nathan Steinle, Patricia Schmidt, and Alberto Vecchio.  Lisa science performance in observations of short-lived signals from massive black hole binary coalescences.  Phys. Rev. D, 107:123026, Jun 2023.
  * [39] Ziren Luo, ZongKuan Guo, Gang Jin, Yueliang Wu, and Wenrui Hu.  A brief analysis to taiji: Science and technology.  Results in Physics, 16:102918, 2020.
  * [40] Michele Vallisneri, Jeff Crowder, and Massimo Tinto.  Sensitivity and parameter-estimation precision for alternate lisa configurations.  Classical and Quantum Gravity, 25(6):065005, mar 2008.
  * [41] Gang Wang, Wei-Tou Ni, Wen-Biao Han, Shu-Cheng Yang, and Xing-Yu Zhong.  Numerical simulation of sky localization for LISA-TAIJI joint observation.  Phys. Rev. D, 102(2):024089, 2020.
  * [42] Gang Wang and Wei-Tou Ni.  Revisiting time delay interferometry for unequal-arm LISA and TAIJI.  Phys. Scripta, 98(7):075005, 2023.
  * [43] Neil J. Cornish and Tyson B. Littenberg.  Tests of bayesian model selection techniques for gravitational wave astronomy.  Phys. Rev. D, 76:083006, Oct 2007.
  * [44] Gang Wang, Zhen Yan, Bin Hu, and Wei-Tou Ni.  Investigating galactic double white dwarfs for the sub-mHz gravitational wave mission ASTROD-GW.  Phys. Rev. D, 107(12):124022, 2023.
  * [45] Valeriya Korol, Elena M Rossi, and Enrico Barausse.  A multimessenger study of the Milky Way’s stellar disc and bulge with LISA, Gaia, and LSST.  Monthly Notices of the Royal Astronomical Society, 483(4):5518–5533, 12 2018.
  * [46] Chang Liu, Wen-Hong Ruan, and Zong-Kuan Guo.  Confusion noise from galactic binaries for taiji.  Phys. Rev. D, 107:064021, Mar 2023.
  * [47] Xue-Hao Zhang, Shao-Dong Zhao, Soumya D. Mohanty, and Yu-Xiao Liu.  Resolving Galactic binaries using a network of space-borne gravitational wave detectors.  Phys. Rev. D, 106(10):102004, 2022.
  * [48] M. I. Jordan and T. M. Mitchell.  Machine learning: Trends, perspectives, and prospects.  Science, 349(6245):255–260, 2015.
  * [49] Elena Cuoco, Jade Powell, Marco Cavaglià, Kendall Ackley, Michał Bejger, Chayan Chatterjee, Michael Coughlin, Scott Coughlin, Paul Easter, Reed Essick, Hunter Gabbard, Timothy Gebhard, Shaon Ghosh, Leïla Haegel, Alberto Iess, David Keitel, Zsuzsa Márka, Szabolcs Márka, Filip Morawski, Tri Nguyen, Rich Ormiston, Michael Pürrer, Massimiliano Razzano, Kai Staats, Gabriele Vajente, and Daniel Williams.  Enhancing gravitational-wave science with machine learning.  Machine Learning: Science and Technology, 2(1):011002, dec 2020.
  * [50] Hunter Gabbard, Chris Messenger, Ik Siong Heng, Francesco Tonolini, and Roderick Murray-Smith.  Bayesian parameter estimation using conditional variational autoencoders for gravitational-wave astronomy.  Nature Physics, 18(1):112–117, 2022.
  * [51] Chayan Chatterjee, Linqing Wen, Kevin Vinsen, Manoj Kovalam, and Amitava Datta.  Using deep learning to localize gravitational wave sources.  Phys. Rev. D, 100:103025, Nov 2019.
  * [52] Stephen R. Green, Christine Simpson, and Jonathan Gair.  Gravitational-wave parameter estimation with autoregressive neural network flows.  Phys. Rev. D, 102:104057, Nov 2020.
  * [53] Stephen R Green and Jonathan Gair.  Complete parameter inference for gw150914 using deep learning.  Machine Learning: Science and Technology, 2(3):03LT01, jun 2021.
  * [54] Arnaud Delaunoy, Antoine Wehenkel, Tanja Hinderer, Samaya Nissanke, Christoph Weniger, Andrew R Williamson, and Gilles Louppe.  Lightning-fast gravitational wave parameter inference through neural amortization, 2020.
  * [55] Plamen G. Krastev, Kiranjyot Gill, V. Ashley Villar, and Edo Berger.  Detection and parameter estimation of gravitational waves from binary neutron-star mergers in real ligo data using deep learning.  Physics Letters B, 815:136161, 2021.
  * [56] Hongyu Shen, E A Huerta, Eamonn O’Shea, Prayush Kumar, and Zhizhen Zhao.  Statistically-informed deep learning for gravitational wave parameter estimation.  Machine Learning: Science and Technology, 3(1):015007, nov 2021.
  * [57] Maximilian Dax, Stephen R. Green, Jonathan Gair, Jakob H. Macke, Alessandra Buonanno, and Bernhard Schölkopf.  Real-time gravitational wave science with neural posterior estimation.  Phys. Rev. Lett., 127:241103, Dec 2021.
  * [58] Maximilian Dax, Stephen R. Green, Jonathan Gair, Michael Pürrer, Jonas Wildberger, Jakob H. Macke, Alessandra Buonanno, and Bernhard Schölkopf.  Neural importance sampling for rapid and reliable gravitational-wave inference.  Phys. Rev. Lett., 130:171403, Apr 2023.
  * [59] Yaron Lipman, Ricky TQ Chen, Heli Ben-Hamu, Maximilian Nickel, and Matt Le.  Flow matching for generative modeling.  arXiv preprint arXiv:2210.02747, 2022.
  * [60] Maximilian Dax, Stephen R. Green, Jonathan Gair, Jakob H. Macke, Alessandra Buonanno, and Bernhard Schölkopf.  Real-time gravitational wave science with neural posterior estimation.  Phys. Rev. Lett., 127:241103, Dec 2021.
  * [61] Yann N. Dauphin, Angela Fan, Michael Auli, and David Grangier.  Language Modeling with Gated Convolutional Networks.  arXiv e-prints, page arXiv:1612.08083, December 2016.
  * [62] Maximilian Dax, Stephen R. Green, Jonathan Gair, Michael Deistler, Bernhard Schölkopf, and Jakob H. Macke.  Group equivariant neural posterior estimation.  11 2021.
  * [63] Michele Vallisneri and Chad R. Galley.  Non-sky-averaged sensitivity curves for space-based gravitational-wave observatories.  Class. Quant. Grav., 29:124015, 2012.
  * [64] John Veitch and Alberto Vecchio.  Bayesian coherent analysis of in-spiral gravitational wave signals with a detector network.  Physical Review D, 81(6):062003, 2010.
  * [65] John Skilling.  Nested sampling for general bayesian computation.  2006\.
  * [66] Gregory Ashton, Moritz Hübner, Paul D Lasky, Colm Talbot, Kendall Ackley, Sylvia Biscoveanu, Qi Chu, Atul Divakarla, Paul J Easter, Boris Goncharov, et al.  Bilby: A user-friendly bayesian inference library for gravitational-wave astronomy.  The Astrophysical Journal Supplement Series, 241(2):27, 2019.