{'fraction_non_alphanumeric': 0.033964728935336384, 'fraction_numerical': 0.014470180374817867, 'mean_word_length': 4.7927823050058205, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Normalizing flows as an avenue to study overlapping gravitational wave signals

Due to its speed after training, machine learning is often envisaged as a solution to a manifold of the issues faced in gravitational-wave astronomy. Demonstrations have been given for various applications in gravitational-wave data analysis. In this work, we focus on a challenging problem faced by third-generation detectors: parameter inference for overlapping signals. Due to the high detection rate and increased duration of the signals, they will start to overlap, possibly making traditional parameter inference techniques difficult to use. Here, we show a proof-of-concept application of normalizing flows to perform parameter estimation on overlapped binary black hole systems.


INTRODUCTION
Over the last few years, the improved sensitivity of the LIGO [1] and Virgo [2] detectors has made the detection of gravitational waves (GWs) originating from compact binary coalescences (CBCs) more and more common, with over 90 detections reported after the third observation run [3]. This new information channel has had major impact in fundamental physics [4], astrophysics [5], and cosmology [6]. Soon, the upgrade of the current detectors and the addition of KAGRA [7][8][9][10][11] and LIGO India [12] to the network of ground-based interferometers will lead to even more detections. In addition, the passage from second-generation (2G) to third-generation (3G) detectors (Einstein Telescope (ET) [13,14] and Cosmic Explorer (CE) [15][16][17]) will lead to an important increase in the number of observed CBCs. These detectors are also projected to have a reduced lower frequency cutoff [18], leading to a longer duration for the signals. Therefore, CBC signals will likely overlap in 3G detectors [19][20][21][22][23].
It has been established that analyzing one of the overlapping signals without accounting for the presence of the other can lead to biases in the recovered posteriors, especially when the merger times of the two events are close [20][21][22][23][24]. These could impact any direct science case for CBCs (e.g. tests of general relativity [25]), but also indirectly related ones such as the hunt for primordial black holes, where subtraction of the foreground sources is required [26][27][28][29][30][31]. In Ref. [32], the authors demonstrate on two overlapped binary black holes (BBHs) how using adapted Bayesian inference can help reduce the biases. They use two methods: (i) hierarchical subtraction, where one analyzes the dominant signal first and subtracts it before analyzing the second event, and (ii) joint parameter estimation, where the two signals are analyzed jointly. They show that the second method is less prone to biases, but computationally heavier. Their analysis also suffers some instabilities, and further upgrades are needed for it to be entirely reliable. An issue also mentioned in this work is the computational time. Indeed, more than 10 5 CBC mergers are expected in the 3G era [20]. So having analyses taking several weeks to complete is not a realistic alternative.
While it is true that traditional methods can be spedup [33][34][35][36], or that quantum computing [37] could potentially be used in the future, the development of frameworks capable of doing complete analyses in very short timescales is crucial for the development of 3G detectors. Therefore, in this work, we propose the first step in that direction, showing how overlapping BBHs can be analyzed with machine learning, and more specifically, a normalizing flow [38][39][40] approach. We start by presenting our machine-learning model. Then we show the setup of our analysis and the results obtained. Finally, we give our conclusions and outline some prospects.

MACHINE LEARNING FOR OVERLAPPING GRAVITATIONAL WAVES
The use of machine learning in GW data analysis has been growing over the last years, having a wide range of applications. See Ref. [41] for a review. A subset of these methods fall under the umbrella of simulation-based inference [42], and are being developed to perform parameter estimation for CBCs [43][44][45][46][47][48][49]. [46][47][48] use normalizing flows to get posterior distributions for BBH parameters. Such methods have been shown to have results close to those from MCMC and nested sampling. Our approach is somewhat similar to theirs, with some notable differences explained below.
Our approach uses continuous conditional normalizing flows [50,51], which is a variant of normalizing flows suited for probabilistic modeling and Bayesian inference problems. Due to the recursive and continuous nature of these models their memory footprint can be quite small [52]. These qualities allow for extensive training on home-grade GPUs while retaining the ability to capture complex distributions. We cover them in more detail below.
Normalizing flows are a method in machine learning through which a neural network can learn the mapping from some simple base distribution p u (u) (a Gaussian, for example) to a more complex final distribution q(θ). This is done through a series of invertible and differentiable transformations, summarized by a function g(θ). However, in our case, it is not sufficient to go from one distribution to the other. We also need to do this conditionally on the GW data that we wish to analyze. To account for this, we use conditional normalizing flows [53], where the transformation functions are dependent on the data d (hence, g = g(θ, d)). A major difference with [53] is that our base distributions are kept static; experiments on toymodels did not show any benefits in having conditional priors. Thus our model g(θ, d) is a trainable conditional bijective function that transforms a simple 30-D Gaussian into a 30-D complex distribution. The bijectivity allows us to express and sample q(θ|d) in terms of g(θ, d) and p u (u) via the change of variable theorem: where det(J g −1 (θ, d)) is the determinant of the Jacobian J g −1 (θ, d) of the transformation. We train the model by minimizing the forward KL-divergence, which is equivalent to maximum likelihood estimation [40,54]. As noted by [48], q(θ|d) should cover the actual (Bayesian) posterior p(θ|d), and asymptotically approach it as training progresses due to the mode-covering nature of the forward KL divergence. A distinctive choice of our method is the continuous nature of the flow, which is linked to the transformation function itself. Neural ordinary differential equations (neural ODEs) [52] are the foundation of continuous normalizing flows; neural ODEs are not represented by a stack of discrete layers but a hypernetwork [55]. Hypernetworks can be understood as regular networks where 'external' inputs such as a (continuous) time or depth variable smoothly changes the output of the network for identical inputs. They can thus represent multiple networks or transformations. In [52], hypernetworks are used to represent ODEs and are trained by using ODEsolvers and clever use of the adjoint sensitivity method. A continuous normalizing flow uses neural ODEs as it transformations.
We will now explain the training of a continuous flow. For clarity we will use h to refer to a continuous transformation and g for a discrete one. If θ(t) represents the samples from the distribution at a given time t, when going from t 1 to t 2 , the continuous normalizing flow obeys an ordinary differential equation The change in likelihood associated with this 'step' differs slightly from equations 1 due to the continuous nature of the flow: Assuming a non-stiff ODE the integration can be performed rapidly with state-of-the-art ode-solvers, which in our case is MALI [56]. In addition, we have to solve a trace instead of a determinant, which reduces the complexity, going from O(D 3 ) to at most O(D 2 ) with D being the dimensionality of posterior space, and speed-ups the computation. Moreover, using continuous normalizing flows removes the need to use coupling layers between transformations, instead all parameter dimensions can be dependent on each other throughout the flow. Combining the continuous and conditional flows leads to continuous conditional normalizing flows, where the conditional consists of the GW d and the time t.
We also need a better data representation than the raw strain to train and analyze the data. Therefore, we follow a similar approach as in [44][45][46][47][48], using a single value decomposition (SVD) [57] as summary statistics, reducing the dimension and the noise content of the data. Each of the 256 generated basis vectors is used as a kernel in 1D convolutions used as an initial layer in a residual convolutional neural network (CNN), enabling one to capture the time variance of the signal. Therefore, we do not need to use a Gibbs sampler to estimate the time of the signal as is done in [46][47][48], and can sample over time like any other variable which allows us to retain the likelihood of the samples. Furthermore, we also use a different representation for the angles. Instead of directly using their values, we project them onto a sphere for the sky location and onto a circle for the other angles. This makes for a better-posed domain for these angles, plays on the strong interpolation capacities of the network, and makes the training step easier.
In the end, our framework combines data representation as a hybrid between classical SVD and CNN, followed by the continuous normalizing flow network. It is worth noting that our entire framework is relatively small compared to the ones presented in [46]. It enables the network to run on lower-end GPUs, but also means the network could be limited in his capacity to model the problem.

DATA AND SETUP
To test the capacity of our framework to deal with overlapping signals, we start with a simplified setup. We consider a network made of the two LIGO detectors and the Virgo detector, at design sensitivity [2,58], and with a lower sensitive frequency of 20Hz. We generate Gaussian noise from their power spectral density (PSD) and inject two precessing BBH mergers generated with the IMRPhenomPv2 waveform [59]. Our data frames have an 8 seconds duration. The chirp mass (M c = (M 1 + M 2 ) 3/5 /(M 1 M 2 ) 1/5 ) is sampled from a uniform distribution between 10M and 100M , and the mass ratio (q = M 2 /M 1 ) from a uniform distribution between 0.125 and 1. We constrain the individual component masses between 5M and 100M . During the data generation, the luminosity distance is kept fixed. It is then rescaled to result in a network signal-to-noise ratio value taken randomly between 10 and 50, sampled from a beta distribution with central value 20. The time of coalescence for the two events is set randomly around a time of reference, with t c ∈ [t ref − 0.05, t ref + 0.05], ensuring that the two BBH merge in the high bias regime [21]. The other parameters are drawn from their usual domain. Table I gives an overview of the parameters and the function from which they are sampled. During the training, we continuously generate the data by sampling the prior distributions for the events and generating a new noise realization for each data frame. The training is stopped when convergence is reached and before over-fitting occurs. Our model trained for about 12 days on a single Nvidia GeForce GTX 1080.

RESULTS
First, we show the corner plots recovered for the masses and sky location of the two events in a pair in Fig. 1. These are representative of our results. One can see that the injected values are within the 90% confidence interval. This is the case for most events, regardless of the relative difference in arrival time or the SNR ratio be-tween the two.
To demonstrate the method's reliability, P-P plots for the two signals recovery are shown in Fig. 2. The P-P plot is constructed by sampling the posteriors of 1000 overlapped events. We then compute in which percentile of the distribution the injected value lies. If everything goes as expected the cumulative density functions should align along the diagonal. One can observe that this is the case for our network. Comparing this to the results given in [46] for single signals, there is a broadening of the shell around the diagonal, showing more variability in the signal recovery. This means that our inference is less accurate than for single signals. Possible origins are the degenerate posteriors, increased complexity of the problem, and the reduced size of our network. In addition, an increased variability has been noted when going from single parameter estimation to joint parameter estimation methods in Bayesian approaches [32].
A trick used to alleviate the problem of degeneracies is to time-order the samples. Indeed, for two BBHs, the likelihood is symmetric in the two events. Therefore, the posteriors can get bimodal [32]. While our training method formally labels one event as A and the other event as B, when the characteristics of the events are close, we may get somewhat bimodal posteriors. This probably also contributes to the higher variability of the P-P plot.
Since parameter estimation for overlapping signals is still an active field of research, it is difficult to compare with traditional methods. While methods have been developed in [32], they are not yet fully stable and take a long time to analyze a BBH system. Therefore, making a statistically significant study comparing the two approaches seems a bit premature at this stage. However, to have some sense of the performances of our network compared to traditional methods, we make 15 injections complying with our network's setup and analyze them with the framework presented in [32]. While this is not enough to make a statistical comparison, we can already identify some trends between the two pipelines. The first is that our machine pipeline typically has broader posteriors than the Bayesian approach. As mentioned in Ref. [32], the classical joint parameter estimation approach can sometimes get overconfident, where the recovered injected value lies outside of the 90% confidence interval. Our method is not confronted with this bottleneck as the broader posterior encapsulates the injected value 1 . Fig. 3 represents the two situations: one where the Bayesian approach finds the event correctly, and one where we see that our ML approach covers the injected values while it does not for the classical approach. The origin of the larger posterior, which is not observed in the single parameter estimation machine learning-based methods, is probably due to the increased complexity of the problem combined with our network's small size. One possible avenue is applying importance sampling after the normalizing flow [48,60]. This would increase the computational time, but the time needed to go from events to samples would still be well below the time taken by the traditional methods. However, such methods can be tricky, and additional modifications to our network could be needed. Finally, an important advantage of our method is its speed. After being trained, it can analyze two overlapping BBH signals in about a second, to compare with O(20days) reported in [32]. While it is difficult to estimate the time gain for other types of signals, such as BNSs or NSBHs, we can expect the inference time after training not to be significantly larger than for BBHs. Since the computation time is a crucial aspect for studies in the 3G era, machine learning approaches seem to be more suited to study realistic scenarios for these detectors.

CONCLUSIONS AND PERSPECTIVES
In this work, we have presented a proof-of-concept machine learning-based method to analyze overlapping BBH signals. We focused on a 2G detector scenario with the two LIGO, and the Virgo detectors at design sensitivity, with a lower frequency cutoff of 20Hz. Our approach is based on continuous normalizing flows.
While also using normalizing flows, as in [44][45][46][47][48], we bring extra modifications that seem to help in the inference task. We represent the data through a mixture of SVD and convolutions, enabling us to sample directly over the events' arrival time, removing the need to use additional Gibbs sampling steps over that parameter and retaining the ability to access the likelihood of a sample. We also move to continuous conditional normalizing flows, reducing the computational cost of the method as we need to solve a trace instead of a determinant when going from one step to the other in the transformation. Finally, we also use a particular representation of the angles, projecting them onto circles (for the phase, the polarization, . . . ) and spheres (for the sky location). We believe that these modifications make our network more supple, enabling it to deal with overlapping signals even in a reduced form.
With this simplified setup, we have shown that our approach is reliable, with posteriors consistent with the injected values. Our method takes about one week to train on a single GPU. After that, it only takes about a second to analyze two overlapped BBHs. While, in reality, other types of CBC mergers can happen, their inference after training should not be significantly longer than for BBHs. We also compared our machine learning method with classical Bayesian methods for overlapping signals. While our scheme leads to wider posteriors, it can correctly recover the injected values, even when the Baysian  FIG. 2. P-P plots for the recovered parameters for the two events in the data. In both cases, the lines align along the diagonal, showing that our method can be trusted. The larger spread can be due to bimodalities in the data, the increased complexity of the problem, and the reduced size of our network.
approach gets overconfident and misses the injection. A possibility to correct for the widened posteriors is to use importance sampling.
Our method's combined reliability and speed show that machine learning is a viable approach to analyze CBC mergers in the 3G era. More interestingly, it would even be possible without needing to account for the development of more powerful computational means and could enable some science-case studies for ET and CE soon. For example, once trained for all possible BBH systems, it could help study the BBH mass function in the 3G era.
Still, one should note that extra improvements are needed before using our method in realistic 3G scenarios. One would first need to change our setup to the 3G detectors, a lower frequency and extreme SNRs that could be encountered. In addition, a wider range of objects Comparison between our approach and the approach from [32]. Our posteriors are generally broader but encapsulate the injected value within the 90% confidence interval. This could be corrected by applying importance sampling on the samples given by our network.
should be accounted for. One should account for higherorder modes and eccentricity that could play a crucial role in the 3G era. Other modifications could also be implemented. Additionally, we need to account for the change in noise realization from one event to the other. Some of these steps, like changing the detectors, should be relatively easy. Others are more complex, as it is hard to perform parameter inference for long-lasting mergers due to the computational burden. So, extra developments in parameter estimation using machine learning would be required to get to the realistic 3G scenario. For overlapping signals, one would also benefit from developments in the classical study of the 3G scenario, such as how to deal with the noise characterization or the types of other events that could come into the data. In the end, there is still work to be done before machine learning can be used in realistic 3G scenarios. However, we believe that this work shows it is an interesting avenue and could be practical on a relatively short time scale.