December 19, 2025
The recently discovered gravitational wave event GW231123 was interpreted as the merger of two black holes with a total mass of 190-265 \(M_\odot\), making it the heaviest such merger detected to date. Whilst much of the post-discovery literature has focused on its astrophysical origins, primary analyses have exhibited considerable discrepancies in the measurement of source properties between waveform models, which cannot reliably be reproduced by simulations. Such discrepancies may arise when an unaccounted overlapping signal is present in the data, or from phenomena that produce similar effects, such as gravitational lensing or overlapping noise artifacts. In this work, we analyse GW231123 using a flexible model that allows for two overlapping signals, and find that it is favoured over the isolated signal model with Bayes factors of \(\sim 10^2 - 10^{4}\), depending on the waveform model. These values lie within the top few per cent of the background distribution. Similar effects are not observed in GW190521, another high-mass event. Under the overlapping signals model, discrepancies in the measurement of source properties between waveform models are largely mitigated. We also find that neglecting an additional signal in overlapping-signal data can lead to discrepancies in the estimated source properties resembling those reported in GW231123. Although the overlapping signal model provides a higher Bayesian evidence, the astrophysical prior probability of two short signals overlapping is low. However, we find that the two recovered sources show similar properties. This, taken with the higher evidence of the two signal model, suggests that gravitational lensing may provide an alternative explanation.
The LIGO–Virgo–KAGRA collaboration [1]–[3] recently reported the detection of the gravitational wave (GW) signal, GW231123_135430 (hereafter referred to as GW231123), from the merger of a binary black hole (BBH) system [4]. With an estimated total mass of 190-265 \(M_\odot\), it became the heaviest BBH system detected to date, surpassing GW190521 [5], whose total mass was 126-170 \(M_\odot\). The individual black hole masses of GW231123 lie within or above the pair-instability black hole mass gap predicted by standard stellar evolution theory [6], [7], making the event particularly intriguing from an astrophysical perspective. Several formation scenarios have been proposed, including hierarchical mergers [8]–[10], mergers in active galactic nucleus disks [11], [12], primordial black holes [13], [14], cosmic strings [15], and evolution of massive stars under specific conditions such as low-metallicity [16], [17], moderate magnetic fields [16], and high spins [18], [19]. Constraints on fundamental physics are also made based on the parameter estimation of GW231123 [20]–[22].
Determining the formation channel of a BBH depends on its source properties, i.e., masses, distance, spins, and such. Measuring these properties for an unusual system such as GW231123 is also a challenging task which heavily relies on how accurately we can model the signal. In this scenario, these are known as waveform models which describe the evolution of spacetime [23]–[26]. For GW231123, the measured source properties using different waveform models show significant discrepancies [4], which cannot be reliably reproduced by the simulations. The scale of discrepancies appears to be larger compared to the expected statistical fluctuations and also larger compared to the systematic differences between the waveform models. This behaviour may be attributed to, but not limited to, inaccuracy of the waveform models, missing physical effects, or noise fluctuations in the data [27]–[29].
An additional source for such discrepancies may be the presence of overlapping signals [30]–[38], i.e., signals arriving at the detectors close to each other in time. Since GW detection rate estimates disfavour such scenarios [39], the state-of-the-art parameter estimation pipelines do not account for overlapping signals [40]. However, when not accounted for, they may bias the measurements of source properties such as masses, distance, and spins [32], [33], [41]. The fainter overlapping signal may be fitted differently by different waveforms, leading to waveform systematics which cannot be reproduced in simulations accounting for only one isolated signal. This seems to be the case for GW231123, i.e., the estimated source properties using different waveform models do not agree, raising the question: could there be an additional signal overlapping with the GW231123? To address this question, we analyse GW231123 with a flexible signal model that accounts for two independent overlapping signals, and find that it is supported by the Bayes factor against the single signal model.
In principle, such a model could also fit other additional power in the data that overlaps with the signal, such as faint noise artefacts [27]. More importantly, significant interest in determining the presence of an additional signal near GW231123 is generated by the recent gravitational lensing interpretation of the event [42]–[45]. Specifically, the scenario where a gravitational wave signal is lensed by a massive object in its path [46], and the two images produced by lensing are not well-separated in time, they can imitate the features of two overlapping signals [47]. This lensing scenario may be viewed as a specific sub-case of our overlapping signal case of more generalised overlapping signal analysis, and may affect the significance of its lensing interpretation [48]. To make a confident detection of a lensed gravitational wave signal, it would be crucial to rule out the possibility of coincidental overlap of signals. We calculate the prior odds of astrophysical overlapping signals in the GW231123 scenario, and find it rules out this possibility. In addition, we find the two recovered sources under the overlapping signal model have similar properties and sky locations. These results point to a potential favour towards the lensing interpretation.
We employ the Bayesian model selection technique to determine which model, two overlapping signals or the isolated signal, fits better to the data corresponding to GW231123. Specifically, we compute the Bayesian evidence \(\mathcal{Z}\), \[\label{eq:evidence} \mathcal{Z} = \int \mathrm{d}\vec{\theta}~\pi(\vec{\theta})~\mathcal{L}(\vec{d} | \vec{\theta}),\tag{1}\] for each model. In Eq. 1 , \(\vec{\theta}\) denotes the model parameters, \(\vec{d}\) denotes the detector strain, \(\pi(\vec{\theta})\) denotes the prior on \(\vec{\theta}\) and \(\mathcal{L}(\vec{d} | \vec{\theta})\) denotes the likelihood of \(\vec{d}\) given \(\vec{\theta}\). For stationary Gaussian noise, the log-likelihood function up to a constant for \(N\) signals can be written as [32], [49], \[\label{eq:likelihood} \ln{\mathcal{L}(\vec{d} | \vec{\theta})} = -\frac{1}{2}\left<\vec{d}-\sum_{i=1}^N \vec{h}_i(\vec{\theta}) \right| \left. \vec{d}-\sum_{i=1}^N \vec{h}_i(\vec{\theta})\right>,\tag{2}\] where \(<\cdot|\cdot>\) denotes the noise-weighted inner product [50] and \(\vec{h}(\vec{\theta})\) denotes the signal model.
The detector strains and the noise power spectral densities were obtained from the publicly available data release corresponding to GW231123 [51]. We analyse a 4 seconds long segment centred at the trigger time of the event. The frequency range of our analysis is \([20,448]\) Hz. We adopt the prior distributions from the public data release [51], with the exceptions that the chirp mass and luminosity distance priors are extended to 250 \(M_\odot\) and 12 Gpc respectively. We restrict our analysis to two overlapping signals owing to the computation cost and the lower odds of more than two signals overlapping. When using the overlapping signals model, the prior on each signal is identical to the isolated signal model. Hence, the prior volume for the overlapping signals model is the square of the one for the isolated signal. We do not include detector calibration uncertainties since calibration envelopes are mostly consistent with zero [52]. We use the Bilby codebase (v2.6.0) in combination with the Dynesty sampler (v2.1.5) to compute the evidence [40], [53].
| XPHM | NRSur | TPHM | XO4a | |
|---|---|---|---|---|
| \(\log_{10} \mathcal{Z}_{\mathcal{S}}\) | \(-3774.14\) | \(-3771.68\) | \(-3769.88\) | \(-3767.95\) |
| \(\log_{10} \mathcal{Z}_{\mathcal{OS}}\) | \(-3769.92\) | \(-3769.72\) | \(-3768.21\) | \(-3767.74\) |
| \(\lbf\) | 4.22 | 1.96 | 0.79 | 0.21 |
We refer to the model with one isolated signal by \(\mathcal{S}\) and two overlapping signals by \(\mathcal{OS}\). The \(\log\) Bayes factor comparing the two models is given by \[\label{eq:bfactor} \log_{10} \mathcal{B}^{\mathcal{OS}}_{\mathcal{S}}= \log_{10}{\mathcal{Z}_\mathcal{OS}} - \log_{10}{\mathcal{Z}_\mathcal{S}}.\tag{3}\] A positive \(\log_{10} \mathcal{B}^{\mathcal{OS}}_{\mathcal{S}}\) indicates that the two overlapping signals model is preferred over the isolated signal model.
We compute the \(\mathcal{B}^{\mathcal{OS}}_{\mathcal{S}}\) for four waveform models: IMRPhenomXPHM-SpinTaylor [23], [54], NRSur7dq4 [24], IMRPhenomTPHM [25], and IMRPhenomXO4a [26]. Hereafter, for brevity, we refer to them by XPHM, NRSur, TPHM, and XO4a, respectively. We find that the overlapping signals model is preferred over the isolated signal model (Table 1) for all four waveforms. We note that the overlapping signals model has twice the number of model parameters compared to the isolated signal model. According to Occam’s razor [55], if two models fit the data equally well, the Bayesian model selection will prefer the model with fewer parameters, penalising the model with more parameters. However, the overlapping signals model overcomes this penalty and is preferred over the isolated signal model.
Comparing the XPHM and NRSur analyses, we find that measurements of masses, distance, and spins are significantly different when using the isolated signal model. These discrepancies are resolved when using the overlapping signals model (Figure 1). Specifically, when using the overlapping signals model, the measurement of the masses, spins, and luminosity distance becomes consistent between the XPHM and NRSur analyses for the louder signals. This does not provide direct proof of the overlapping-signal hypothesis, i.e. it is not a sufficient condition. Rather, it may be viewed as a necessary consequence of that hypothesis, and provides an alternative to waveform systematics as the cause of the posterior discrepancies in the single-signal model. For the fainter signals, the low signal-to-noise-ratio (SNR) leads to significantly broader posterior distributions, though they still agree well between the two models. We find that the parameter estimates of the fainter signal, including its masses and sky location, overlap with those of the louder one, suggesting that the two signals could share the same source properties. We also find the difference in the arrival time to be 20 milliseconds between the two signals of the overlapping signals model. These findings are consistent with the lensing analyses of the event, where one typically measures the arrival time difference between two lensed images [42]–[44].
While the component masses are broadly consistent between the overlapping-signal and single-signal models, the inferred spin magnitudes tend to be lower in the overlapping-signal case. Since high spins form the basis of many interpretations of the astrophysical origin of GW231123 [8]–[12], adopting the overlapping-signal results may reduce the statistical significance of those interpretations. In addition, waveform uncertainties are generally smaller in the low-spin region of parameter space, which may help to explain the reduced posterior discrepancies between XPHM and NRSur.
Comparing the TPHM and XO4a analyses, we do not find any significant difference in the measurement of source properties, whether we use the isolated or overlapping signal model. This is consistent with the \(\mathcal{B}^{\mathcal{OS}}_{\mathcal{S}}\) values for XO4a and TPHM analyses, which are much smaller compared to the ones obtained from XPHM and NRSur analyses.
The difference in \(\mathcal{B}^{\mathcal{OS}}_{\mathcal{S}}\) values between different waveforms may be attributed to waveform systematics. Among all four waveforms, the TPHM and XO4a waveforms show stronger evidence in favour of the isolated signal model for GW231123 (first row of Table 1). If we assume that there is an additional signal present in the data overlapping with GW231123, the values in Table 1 imply that the XO4a and TPHM waveforms are more susceptible when fitting to a distorted signal, compared to the other two waveforms.
To verify this, we simulate two overlapping signals in zero noise and recover them with the isolated signal model. We perform four analyses: the simulated data consists of two overlapping signals generated by the NRSur waveform, which are then recovered with the isolated signal model using each of the four waveform models separately. We select the maximum likelihood parameters of the NRSur analysis of GW231123 with the overlapping signal model as injection parameters.
These simulations corroborate the results of the GW231123 analysis (Table 2). Specifically, we find that when we attempt to recover two overlapping signals with the isolated signal model, the XO4a waveform gives the largest evidence, followed by TPHM, NRSur, and XPHM. Besides the evidence values, we show the estimated source properties from these simulations in Figure 2. Specifically, we show the detector frame total mass, mass ratio, and the luminosity distance. We observe discrepancies in the measurements between different waveforms, which are similar to the ones observed for GW231123 [4], indicating that the presence of an additional signal could cause this behaviour.
| XPHM | NRSur | TPHM | XO4a | |
|---|---|---|---|---|
| \(\log_{10} \mathcal{Z}_{\mathcal{S}}\) | 171.08 | 177.81 | 179.60 | 180.21 |


Figure 2: Posterior distributions when we recover two simulated overlapping signals using the isolated signal model with four different waveforms. The top panel shows the measurements of detector frame total mass and mass ratio; with the contours showing 90% confidence interval. The bottom panel show the measurements of luminosity distance. With these simulations, we are able to reproduce the discrepancies in the estimated source properties between different waveforms reported in the GW231123 data release (c.f. Figures 7 and 8 in Ref. [4])..
To estimate the significance of \(\log_{10} \mathcal{B}^{\mathcal{OS}}_{\mathcal{S}}\) comparing the overlapping and the isolated signal models, we estimate the corresponding background statistic. Specifically, we want to find out how often an isolated BBH signal could imitate the features of two overlapping signals due to the presence of noise fluctuations. To that end, we simulate a set of 100 BBH signals; 50 of which are chosen from the up-to-date BBH population distribution [39], and the remaining 50 are chosen from the posterior distribution of GW231123 [51]. The latter subset is to ensure we cover the parameter space of high total mass and high spin BBHs. To compute the background statistic, we perform parameter estimation on these isolated, simulated BBH signals using the isolated signal and the overlapping signals models.
When estimating background using GW231123-like signals, we rely on NRSur waveform for injection and recovery. We consider GW231123-like signals in simulated Gaussian noise as well as detector noise around GW231123 but excluding the event. The latter choice is to estimate how the noise artifact may falsely indicate the presence of overlapping signals. When estimating the background using astrophysical BBH signals, we rely on XPHM waveform for injection and recovery due to limited parameter-space coverage of NRSur [24], [56]. For the astrophysical background, we only consider simulated Gaussian noise. The case of astrophysical BBH in detector noise should give a conservative estimate of the background compared to the GW231123-like signals in detector noise so we do not consider it here. We present the collective results of the background estimation in Figure 3.
In Figure 3, the green (brown) line shows the foreground \(\log_{10} \mathcal{B}^{\mathcal{OS}}_{\mathcal{S}}\) for NRSur (XPHM) analyses. The black lines show the background statistic for events drawn from an astrophysical population (solid) and GW231123-like events (dashed and dotted lines for real detector and Gaussian noise, respectively). We find that all of the astrophysical background is below the foreground statistic of the XPHM analysis (brown), which indicates strong support for the overlapping signals model when using the XPHM waveform. Further, the large difference between the brown and green lines may be caused by waveform systematics between NRSur and XPHM. The background estimated using the detector noise is consistently above that of simulated Gaussian noise, indicating that the noise artefacts play an important role when comparing the two models, as seen for our GW231123-like injections. Additionally, we find that the GW231123-like signals are more likely to falsely imitate the signs of overlapping signals than the broader BBH population.
Although background studies usually focus on false positives caused by noise fluctuations, we note that waveform systematics can also produce false positives in overlapping-signal analyses, as waveform errors may be absorbed by an additional signal. This may be particularly relevant for the high-mass scenario of GW231123, since the signal is dominated by the merger stage, where waveform uncertainties are largest. As pointed out by [4], a more accurate waveform model would likely be required before drawing confident conclusions for GW231123.
Similarly, we analyse GW190521 [5], the heaviest BBH system prior to the detection of GW231123, using the isolated signal and the two overlapping signals model. We obtain \(\log_{10} \mathcal{B}^{\mathcal{OS}}_{\mathcal{S}}\) of 0.27 for XPHM, \(-0.02\) for NRSur, \(0.56\) for TPHM, and \(0.54\) for XO4a waveform. The values of \(\mathcal{B}^{\mathcal{OS}}_{\mathcal{S}}\) across different waveforms are more comparable, indicating the effect of waveform systematics may be negligible for GW190521. We also note that the inferred luminosity distance of the fainter signal consistently shifts toward higher values, which, in combination with the Bayes factors, disfavour the presence of an additional signal for this event.
Although the Bayesian model selection favours the overlapping signals model, the prior odds of such an occurrence are expected to be low [33]. In addition, both signals in the overlapping model lie in the high-mass tail of the astrophysical compact binary population [39], further reducing the prior odds. To get an estimate of the prior odds, we adopt the high-mass BBH merger rate inferred from the observation of GW190521, \(\mathcal{R} = 0.13^{+0.30}_{-0.11} \mathrm{Gpc}^{-3} \mathrm{yr}^{-1}\), and assume this merger rate is constant across the comoving volume. Given our prior bounds (\(z_\mathrm{max}=1.6\) and \(\Delta T = 0.4\)s), the probability of observing another high-mass BBH event within a time window \(\Delta T\) when one such signal is already present can be expressed using Poisson statistics as [33] \[p(\mathcal{OS}|\mathcal{S}) = 1- \exp \left( -\Delta T\int_0^{z_\mathrm{max}} \frac{\mathcal{R}dV_c}{(1+z)dz}dz \right),\] which yields \(p(\mathcal{OS}|\mathcal{S}) \sim 10^{-7}\). Multiplying by the Bayes factor, we obtain the final posterior odds ratio of observing two overlapping BBH signals to be \(\lesssim 10^{-3}\), an exceedingly small number. Since the overlapping signal model accommodates a number of phenomena, the low prior odds of two overlapping signals implies the other ones, such as gravitational lensing or overlapping glitches, would be sensible alternatives.
We perform Bayesian model selection analyses on GW231123 to determine which model, between the overlapping signals and the isolated signal, is favoured by the data. For all waveform models used for the analysis, we find that the flexible signal model allowing for two overlapping gravitational wave signals is favoured over the isolated signal model when analysing GW231123. The overlapping signals model largely mitigates the discrepancies in measurements of mass, distance, and spins between different waveforms. In an additional set of simulations, we are able to reproduce the discrepancies between different waveform models similar to GW231123, when neglecting the presence of the overlapping signals.
However, interpreting high-mass binary black hole systems is challenging because only the last few cycles of the signal are observed, and a variety of models may fit the data, thereby affecting the inferred source properties. We cannot confidently establish an overlapping-signal interpretation for GW231123, partly because noise features may mimic signatures of overlapping signals, as indicated by the background estimates, and partly because of waveform systematics, as indicated by the comparison between different waveform analyses. That said, our findings on GW231123 remain valuable, as this is the first GW event for which the overlapping-signal model provides a better fit than the single-signal model; no comparable behaviour is seen in another similarly high-mass BBH event, GW190521.
While the detection of two overlapping gravitational wave signals is disfavoured by the current detection rate estimates, the overlapping signal model could be viewed as a superset of other potential phenomena which could be detected with current detector sensitivities. One possibility is the presence of faint noise artifacts near the signal [27], although similar noise transients occurring coherently in both detectors are not commonly observed. The overlapping-signal model can also describe the case of gravitationally lensed signals in which two images overlap in time due to a small arrival-time delay. This degeneracy between lensing and overlapping signals makes it important to confidently rule out the overlapping-signal scenario. In this work, we find that the overlapping signal interpretation is disfavoured by astrophysical rate estimates, while the inferred source properties remain consistent with a possible lensing interpretation. These results therefore motivate further investigation of the gravitational-lensing scenario for GW231123.
As the detector sensitivity increases and as we transition towards the era of third-generation gravitational wave detectors [57], [58], it will become possible to observe subtle effects from both lensing and overlapping signals. In such scenarios, it becomes crucial to develop a methodology that can distinguish between the two effects. Our work takes a first step in that direction.
The authors would like to thank Jonathan Gair, Yifan Wang, Zhen Pan, and the LIGO–Virgo–KAGRA parameter estimation group for helpful discussions and suggestions. The authors are grateful for the computing resources provided by Cardiff University funded by STFC grant ST/I006285/1 and ST/V005618/1, and by the Consortium des Équipements de Calcul Intensif (CÉCI) funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11 and by the Walloon Region. QH and JV are supported by STFC grant ST/Y004256/1. JH is a FRIA grantee of the Fonds de la Recherche Scientifique - FNRS. HN, MW, and CVDB are supported by the research programme of the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (Netherlands Organisation for Scientific Research, NWO). HN and CVDB further acknowledge the support of NWO through the grant number OCENW.XL21.XL21.038. We also acknowledge the use of the computing infrastructure of Nikhef, which is co-supported by the NWO. This material is based upon work supported by NSF’s LIGO Laboratory which is a major facility fully funded by the National Science Foundation.