1 sec:Introduction↩︎

The detection of gravitational waves by the LIGO-Virgo-KAGRA (LVK) detector network has opened a new era of gravitational wave astrophysics [1]–[9]. During the first three observing runs these detectors have already detected a large collection of signals, \({\cal O}(100)\), from binary mergers of compact objects, namely black holes and neutron stars, together with an estimation of their physical parameters. These signals have been mainly detected and studied through the matched-filtering of the experimental data [10]–[18] with theoretical waveform template banks [12], [19], [20].

Even though templates from the coalescence of black holes and neutron stars have dominated the search landscape as candidates for gravitational wave sources, the possibility of alternative compact objects, known generically as exotic compact objects (ECO), has generated an increasing interest in recent years [21]. In particular, much scientific attention has been oriented towards those objects that could generate gravitational-wave signals resembling those of black holes, which could generate confusion in their classification. These objects are typically known as “black hole mimickers”.

A particularly important class of ECOs are those based on gravitationally bound states of ultralight bosonic fields [22], [23]. These fields could appear in particle physics models, such as the string axiverse [24], [25] or as extensions of the Standard Model [26]. Ultralight bosonic fields are capable of forming stationary states, resembling hydrogen orbitals [27], which can also be rotating [28], [29]. In the case of massive complex vector fields — Proca fields — the objects are known as Proca stars [30].

While sharing some similarities, scalar and Proca bosonic stars have distinct behaviors regarding their formation and stability [31]. In particular, some configurations of rotating bosonic stars exhibit a bar-mode instability [32] which can be prevented through the addition of non-linear interaction terms in their action [33], [34] or considering multi-state bosonic configurations [35]. On the other hand, Proca stars are known to have sufficiently generic mechanisms of formation though the so-called “gravitational cooling” mechanism [36], [37] and a stable ground state [38]. A general review on boson stars can be found in [39].

Numerical evolutions of Proca stars have been performed at the fully non-linear regime in 3+1 numerical relativity, both with single stars [40] and head-on collisions [41], [42]. This has allowed the use of the gravitational waveforms extracted from the simulations to compare with actual detections and perform parameter estimations (PE) [43], [44]. In particular, statistical inference on the event GW190521 resulted in a very good agreement, even with a slightly higher Bayes factor than with black hole templates [45].

The standard procedure for Bayesian PE involve sampling over a large number of waveform templates generated with different physical parameters, typically on the order of millions. As numerical relativity simulations are computationally expensive, on the order of thousands of CPU-hours per simulation, generating all the templates by direct simulation is not possible at the practical level. For the specific case of binary black hole (BBH) mergers, the development of waveform templates (or “approximants”) based upon approximations to the equations of motion of the two-body-problem in general relativity, has already a long history spanning over two decades [46]. For the three main model families covering the whole inspiral-merger-ringdown process, namely SEOBNR [47]–[50], TEOBResumS [51]–[53] and IMRPhenom [54]–[57], waveform templates are either calibrated to numerical relativity using analytical or semi-analytical expressions or they are obtained through the “hybridization” of waveforms from numerical relativity simulations and post-Newtonian approximations.

A fourth and more recent line of work is the construction of surrogate waveform models, first proposed in the seminal work of [58]. Such models have been shown to reach an accuracy comparable to numerical relativity and they will be the focus of our work (albeit in the context of Proca-star mergers). For non-spinning quasi-circular BBH mergers up to a mass-ratio \(q=8\), the first attempt to build a surrogate model employing waveforms from numerical relativity was performed by [59]. This work was soon followed by the first surrogate model for precessing systems [60], NRSur4ds2, built from a set of 276 precessing BBH simulations, and by the corresponding extensions to account for unequal masses and generic spins NRSur7dq2 [61] and NRSur7dq4 [62]. The latter was calibrated with 1528 numerical relativity precessing simulations with mass-ratio \(q\leq 4\) and spin magnitudes up to 0.8, covering all harmonics with \(l\leq 4\), retaining all intrinsic degrees of freedom of quasi-circular BBH systems. NRSur7dq4 is the most accurate waveform model for quasi-circular BBH with generic spins. This model has been employed to analyze several exceptional events from the LVK Collaboration [4], [63], to extract measurements of novel observables like orbital precession [64], gravitational-wave recoils [65], [66] or the Chern-Pontryagin pseudo-scalar [67].

In addition, Islam et al. [68] recently used NRSur7dq4 to re-analyze 42 events detected by the LVK Collaboration. Moreover, for aligned-spin binaries, a surrogate model for hybridized numerical relativity waveforms was constructed in [69] using 104 simulations, employing Effective-One-Body for hybridizing the orbital frequency and post-Newtonian results for hybridizing the harmonic amplitudes. This was extended in [70] who developed a non-precessing 2-dimensional surrogate model of hybridized waveforms up to mass-ratio \(q=15\), NRHybSur2dq15, in order to build a template waveform for signal GW190814 [71]. It is also worth mentioning that efforts based on Machine Learning, namely Gaussian Process Regression, have been used to build waveform surrogate models for both non-precessing and precessing BBH systems at points of the parameter space not covered by numerical relativity simulations [72], [73].

The progress in surrogate waveform models for BBH mergers motivates the development of such models for ECOs, in particular if a systematic search for Proca stars is to be pursued. These models should be able to produce accurate waveform estimates, taking as input the physical parameters of the system, using modest computational resources. Some attempts have been made in this direction by application of generative adversarial networks (GAN) [74], already producing promising results. In this paper, we perform a first approach to the use of surrogate models for PE in Proca star binaries by restricting to head-on collisions of Proca stars with equal masses. The choice of head-on collisions instead of quasi-circular inspirals is entirely motivated by the much higher technical difficulty and computational cost of the Proca inspirals, which makes it challenging to perform numerical-relativity simulations of such systems at the present time.

Therefore, in our case, the system depends on a single physical parameter, the normalised frequency of the system \(\omega/\mu_V\), where \(\mu_V\) denotes the mass of the ultralight vector boson forming the stars. For the PE procedure to confidently resolve the physical parameters of the system, the error introduced by the surrogate model must be much lower than the typical differences between templates in the sampling range. For this reason, we analyze a variety of possible surrogate architectures that can be used to fit the simulation data. We compare their performances on a testing set and we finally select the most reliable one for the actual Bayesian analysis. The models in this study are data-driven, based entirely on a catalog of simulations [42], with no built-in physical information.

In the first class of models, which we call two-phase architectures, we encode the waveform data into a feature space of lower dimension. This drastically reduces the amount of numerical parameters that need to be fitted, thus reducing the model size. The parameters in the feature space (sometimes called latent space) can then be interpolated in terms of \(\omega/\mu_V\). As we have only one independent variable, we use standard cubic spline interpolation. The algorithms that we analyze for the dimensional reduction phase are singular value decomposition, empirical interpolant representation [75] and a deep convolutional autoencoder. The first two algorithms decompose the data as a linear combination of a suitable reduced basis, while the last one uses a nonlinear mapping based on neural networks. The second type of model is a more straightforward approach, which we call monolithic, as a single neural network is used to map the frequency values \(\omega/\mu_V\) to the temporal values of the waveform.

We use our model to perform Bayesian parameter inference on both simulated and real data. In the first case, we perform parameter inference on numerical relativity waveforms injected in zero noise, demonstrating that our surrogate can accurately estimate the parameters of such injections. In the second case, we perform parameter inference on the data containing the signal GW190521, showing that our surrogates faithfully match those obtained through numerical relativity waveforms, presented in [43], [45], while using significantly less computational resources. This indicates that this type of data-driven surrogate models are a very promising tool to encapsulate the catalogs of numerical boson-star waveforms in order to be efficiently used in detection data analysis.

The paper is structured as follows. In Section 2 we present the generation of the training data and its preparation before it can be fed to the surrogate models. We also discuss the criteria used for the evaluation and comparison of the different models. Section 3 describes the different architectures of two-stage models, i.e., those who can be divided into encoder and translator modules. A different type of models, the monolithic surrogates, are exposed in Section 4. The results of model evaluations, as well as the Bayesian parameter estimations performed with the models, are discussed in Section 5. We conclude in Section 6.

2 The Dataset↩︎

2.1 The Proca Merger Catalog↩︎

The models are trained on a catalog of 59 equal-mass simulations from [42], of Proca stars in Einstein-Proca theory \[\mathcal{S} = \int d^4 x \sqrt{-g} \left(\frac{R}{16\pi} - \frac{1}{4} F_{\mu\nu} \bar F^{\mu\nu} - \frac{1}{2} \mu^2 A_\mu \bar A^\mu \right) \; ,\] where \(A_{\mu}\) and \(F_{\mu\nu}\) are the Proca potential and field strength, respectively, \(\mu\) is the Proca field mass parameter, and the bar denotes complex conjugation. We also take \(c = G = 1\). The interested reader is addressed to [42] and references therein for details on the formalism. The simulations were performed in the Einstein Toolkit [76], which uses the Cactus framework [77], under the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation using the McLachlan [78], [79] and Lean [80] thorns. Details on the simulations and construction of initial data can be found in [30], [40], [41]. The thorn containing the Proca equations was implemented in [81] and is available in [82].

Each simulation produces gravitational waveforms for the Newman-Penrose scalar \(\psi_4\), which are decomposed in spherical harmonic modes. Our study will be centered on the fundamental \((l,m) = (2,2)\) and \((l,m) = (2,0)\) modes. An example of the waveforms, together with their prediction by one of the surrogates presented below, is shown in Figure 1.

Even though the ADM mass of the initial data may vary due to the construction of the Proca star binaries, the numerical relativity waveforms are rescaled so that \(M_{\rm{ADM}} = 1\) before they are treated in the surrogate pipelines.

2.2 Pre-Processing of the Data↩︎

The complex \(\psi_4\) waveforms of the original dataset are subjected to a pre-processing pipeline before they are submitted to the actual model. There are a number of reasons for this. Firstly, the algorithms will always tend to be more effective when the samples are similar to each other. For instance, having large time shifts between otherwise very similar waveforms can unnecessarily increase the cost of the fitting. Secondly, many machine learning algorithms work best with data of numerical value of order one. As the values of \(\psi_4\) vary usually on the order of \(10^{-4}\), it is highly advisable to rescale the data before making any learning attempt on it. Additionally, numerical relativity simulations usually include a burst of unphysical junk radiation at the beginning of the waveform, which has to be removed. And last but not least, it is highly convenient to have all waveforms on the same time grid, which makes it possible to treat them on equal footing. To this end, a Gaussian curve is fitted to the modulus \(|\psi_4|\) of the waves, as \[f(t) = A e^{- (t - t_0)^2/s^2}\; .\] Then, the waveforms are time-shifted by \(t_0\) so that all the waves are centered at \(t = 0\). Then, they are divided by \(A\) so their range of oscillation becomes of order one. In order to perform the transformations, the original waveform is interpolated by cubic splines and re-evaluated on a fixed grid of \(N_\psi\) equally spaced points between -150 and 150. \(\psi_4\) is set to zero for \(t < 100\), to remove the junk radiation. Figure 2 shows an example of the Gaussian curve fit to the amplitude of a wave after the transformations.

In order for the surrogate to be used for parameter estimation, the relative time positions of the different modes have to be preserved. The values of the time shift and amplitude \((t_0, A)\) are stored in the dataset, as they also need to be interpolated in order to recover the information later. Once the model is evaluated, we multiply the waveforms by their respective amplitudes \(A\) and we time-shift them so that the peak of the reference mode \((l,m) = (2,2)\) lies at \(t = 0\). This way, the original relative time positions of all the modes are preserved.

The method of the Gaussian fit for the identification of the peak and amplitude of the waveform, as opposed to just taking the absolute maximum of \(|\psi_4|\), guarantees that the paramters \(A\) and \(t_0\) will vary continuously with \(\omega/\mu_V\). The absolute maximum of \(|\psi_4|\), on the other hand, can suffer discontinuous jumps between relative maxima, which would create discontinuous jumps on the data, thus damaging the performance of the fitting algorithms.

2.3 Model Testing↩︎

As in any data-driven learning algorithm, the surrogate models are prone to suffer from overfitting. In other words, the models will tend to perform better on the data they have been trained than on new data that they have not seen before. An example of this is shown in Fig. 3. Therefore, before the surrogate model can be used on parameter estimation samplers, an estimate must be done of its reliability when generating waveforms outside of its own training dataset.

The full waveform dataset from the original Proca catalog is divided into training and testing datasets with 80% and 20% of the samples, respectively. All the model fittings will be performed on the training dataset, and the test set will be reserved as a final evaluation of the model. As a metric for the accuracy of the model we employ the mismatch [83]–[85] between the predicted and the actual waveform, defined as \(1 - M(\psi_1, \psi_2)\). Here, the match \(M(\psi_1, \psi_2)\) is given by \[M(\psi_1,\psi_2) = \max_{t_s,\phi}\frac{\langle \psi_1 | \psi_2 \rangle}{\sqrt{\langle \psi_1 | \psi_1 \rangle \langle \psi_2 | \psi_2 \rangle}},\] with \(\langle \psi_1 | \psi_2 \rangle\) the inner product defined as \[\langle \psi_1 | \psi_2 \rangle = 4 \times {\rm Re}\int_{\nu_0}^{\nu_{\rm 1}} \frac{\tilde{\psi_1}(\nu) \, \tilde{\psi_2}^{*}(\nu)}{S_n(\nu)}\,d\nu\,,\] where \(\tilde{\psi}_i(\nu)\) are the Fourier transforms, and \(\tilde{\psi}_i^*(\nu)\) their complex conjugates. The match is maximized over relative time shifts \(t_s\) and phases \(\phi\) between \(\psi_1\) and \(\psi_2\). \(S_n(f)\) is the one-sided noise power spectral density (PSD), which we take to be flat.

2.4 Model Selection↩︎

Before the final evaluation on the test set, we have to choose among a number of different architectures of models, as well as the hyperparameters intrinsic to each architecture. It is not convenient, however, to check the performance of each configuration on the final test set. If we did so, we could again be selecting the architecture that best suits our particular test set, without it generalizing to arbitrary data.

While this could be easily solved by splitting the training set into two subsets —one for yet again training and the other for validation— it would further reduce the number of samples in each subset. Given the limited number of simulations in our catalog, it would be convenient to use a scheme where we could extract information from all the data in the original training set. A good solution for this is \(k\)-Fold cross-validation on the training set. In this scheme, the training set is itself divided in \(k\) (approximately) equal parts, called folds. We then train \(k\) versions of the model. For each version, we train the model on \(k-1\) of the folds, and we validate it on the remaining fold. With this method, all of the training data samples have been used as a validation of the model architecture, without ever being used both in the training and validation sets. It is important to keep in mind that the final testing set does not intervene at all during the model selection phase. In this paper we will always use \(k=5\).

In Fig. 4 we depict the mismatches of all the \((l,m) = (2,2)\) waveforms in the dataset, evaluated on their corresponding model, in this case based on the empirical interpolant representation (see Sections 3.1 and 3.2). The waveforms of the training set, which we depict as circular dots, are divided in the 5 cross-validation folds, corresponding to the 5 different colors of the dots. For the validation on each of the folds, the model has been trained on the remaining 4 folds. On the other hand, for the waveforms in the final test set, the model has been trained on the full training set (all 5 folds together). The distribution of mismatches is quite similar for all the subsets, which indicates that our statistical estimates of the surrogate performance are consistent with each other. The points at the edge of the domain, close to \(\omega = 0.93\), have significantly higher mismatches as extrapolation is being used.

For the sake of consistency, in all of the dataset splittings, we keep the different modes of each simulation together. In other words, if the mode \((l,m) = (2,2)\) is in one particular subset, then so are the other values of \((l,m)\).

3 Two-Stage Models↩︎

The two-stage surrogate models are divided in a pipeline of two processes, that we chose to name as encoder-decoder and translator.

• The encoder-decoder phase is an invertible process of dimensional reduction from the \(N_\psi\) complex values of \(\psi_4^i\) to a feature space of lower dimension \(N_Z\), which can be complex or real depending on the algorithm.

• The translator phase interpolates the features \(Z^j\), as well as the amplitude and peak position, as a function of the physical parameters of the model. In this case, the only physical parameter is the frequency \(\omega/\mu_V\) of the merging Proca stars.

Schematically, \[\begin{align} \mathbb{R} &\rightarrow \mathbb{C}^{N_Z} \rightarrow \mathbb{C}^{N_{\psi}} \\ \omega &\mapsto Z^j \mapsto \psi_4^i \; . \end{align}\] When the model is evaluated for parameter estimation, the process is then reversed. First, the physical frequencies are used in the translator to obtain the coordinates \(Z^j\) in feature space, which then are passed through the decoder to recover the values of the Newman-Penrose scalar \(\psi_4\).

3.1 Singular Value Decomposition↩︎

Our first approach for dimensional reduction is singular value decomposition (SVD). This produces an orthonormal complex basis for our space of waveforms, in such a way that successive basis vectors are less and less important in the accurate description of the data. The notion of importance is quantitatively encoded in the singular values of the basis. Namely, if we construct an \(N_\psi \times N_t\) matrix \(M\), whose \(N_t\) columns are the training waveforms, we can decompose it as \[M = U \Sigma V^\dagger\; .\] Here the columns of \(U\) are the adapted orthonormal basis \(u_i\), and the diagonal elements of \(\Sigma\) are their corresponding singular values \(\sigma_i\). By taking the first \(N_Z\) of such vectors, we can express a compressed version of any waveform as a linear combination of them. Their \(N_Z\) coordinates define the feature space of the model. The fact that \(U\) is a unitary matrix greatly simplifies its inversion, as \(U^{-1} = U^\dagger\), so the encoding-decoding of waveforms becomes trivial.

The main hyperparameter on the SVD models is the number of vectors in the reduced basis \(N_Z\). We have thus performed an analysis of such basis size on the performance of the model, as shown in Fig. 5. The averaged mismatch after the \(k\)-Fold cross validation process decays approximately as an exponential law for \(1 \leq N_Z \leq 7\), and remains fairly constant after that. In other words, having more than 7 vectors in the reduced basis does not seem to increase the precision of the model. We have performed our parameter estimation with a basis of 10 vectors, to keep a safe margin.

3.2 Empirical Interpolant Representation↩︎

The empirical interpolant representation (EIR) is an additional step we take after SVD to obtain a more convenient basis for the space of waveforms in our model. The idea behind EIR algorithms is to choose a reduced basis \(B_j\) of the waveform space and a subset of indices \[\mathcal{I} \subset \mathbb{Z}_{N_\psi} \, , \quad |\mathcal{I}| = N_Z\] such that \[B^i_j = \delta^i_j \quad \forall j \in \mathcal{I} \; .\] In other words, the coordinates of \(\psi_4(t)\) in the basis \(B_j\) are its values at certain selected times \(\psi_4(t = T^j) = \psi_4^j\), for \(j \subset \mathcal{I}\), as \[\psi_4^i = \sum_{j \in \mathcal{I}} \psi_4^j \, B^i_j\; .\] This process is performed by the Python library rompy [75], [86], which takes as an input the SVD basis \(u_i\, ,(i = 1, \dots, N_Z)\) and computes the basis vectors \(B_j\) as well as the relevant indices \(\mathcal{I}\).

The encoding procedure of the waveforms into the reduced feature space is then straightforward, as it just involves taking the elements of \(\psi_4\) in the selected positions \(\mathcal{I}\). The decoding is also simple, just multiplying the feature space coordinates by the basis vectors \(B_j\).

In a similar way as in the case of SVD, we have studied the mismatch dependence on the number of EIR basis vectors \(N_Z\), as shown in Fig. 6. The plot looks very similar to that of SVD, with an exponential decay which stabilizes for \(N_z \geq 7\).

3.3 Deep Convolutional Autoencoder↩︎

Autoencoders provide a nonlinear alternative to EIR dimensional reduction. In this case, we set up two 1D convolutional neural networks, which will act as an encoder and a decoder, respectively. The encoder will learn a map from the real and imaginary parts of the waveform \(\psi_4^i\) to the feature space \(Z^j\), which is now taken to be real. In other words, the encoder will become a nonlinear mapping from \(\mathbb{C}^{N_\psi}\) to \(\mathbb{R}^{N_Z}\). On the other hand, the decoder will learn the inverse mapping, trying to recover \(\psi_4^i\) from the features \(Z^j\). In the context of autoencoders, the reduced feature space \(Z^j\) is typically known as latent space.

The training process is performed on a concatenation of the encoder and decoder networks, in such a way that the goal of the full system is to learn the identity function, with the output being as close as possible to the input. This process is nontrivial as the number of parameters has to be drastically reduced from \(N_\psi\) to \(N_Z\) as the information is transmitted through the bottleneck between the networks. A schematic representation of this configuration is shown in Fig. 7. We define the loss function of the system as \[\mathcal{L} = \frac{1}{2N_\psi} \sum_{i=1}^{N_\psi} |\psi_4^i - \tilde{\psi}_4^i|^2\, ,\] with \(\psi_4^i\) and \(\tilde{\psi}_4^i\) the original and recovered waveforms, respectively. The loss is also averaged over the instances in the training batch. After the training, the encoder and decoder modules can be used separately to map the waveforms to their coordinates in the latent space and vice-versa. Figure 8 shows the embedding of the waveforms in a latent space of \(N_Z=2\), where the model has spontaneously ordered the waveforms by their value of \(\omega/\mu_V\). In this case the mapping is nonlinear, so the number of coordinates in the latent space is usually lower than the number of basis vectors that are needed in SVD or EIR encoding.

The structure of the networks is fairly symmetric between encoder-decoder. The encoder receives the training minibatch of waveforms as a tensor of dimensions \((N_B, 2, N_\psi)\). \(N_B = 16\) is the number of instances in the training minibatch, 2 is the number of channels (corresponding to the real and imaginary parts of the wave), and \(N_\psi = 512\) is the length of the time series containing the wave. The batch is first passed through a 1-dimensional convolutional layer of 16 filters of kernel size 4, stride 2 and padding 1. After that a leaky ReLu activation function, with a negative slope of 0.2, is applied. Then the batch goes through a series of 6 more convolutions with identical kernel configurations, but each time duplicating the number of filters. After each convolution, we apply a batch normalization and another leaky ReLu with a negative slope of 0.2. This configuration of filters reduces the length of the data samples by 2 at each convolution, while doubling the number of channels. Finally a last convolution is applied, this time with a stride of 1 and padding reduced to 0, with \(N_Z\) filters. This reduces the minibatch to a tensor of size \((N_B, N_Z, 1)\), where the \(N_Z\) channel values are taken to be the coordinates in the latent space.

The architecture of the decoder is quite exactly inverse from the encoder. The input minibatch of dimensions \((N_B, N_Z, 1)\) is passed through a transposed convolution with 64 filters of kernel size 4, stride 1 and padding 0. Then, a series of 7 more (transposed) convolutions are applied, now with a stride of 2 and padding of 1, each one duplicating the size of the data instances and halving the number of filters. After each convolution (except the very last one), we apply batch normalization and a ReLU activation. The last convolution actually reduces the number of filters down to 2 (real and imaginary parts of \(\psi_4\)), resulting in an output of dimensions \((N_B, 2, N_\psi)\), as expected.

The training is performed on a Nvidia RTX 4090 24GB GPU over minibatches of 16 waveforms going through 500 iteration epochs with a learning rate of \(10^{-3}\), using an Adam optimizer with parameters \((\beta_1, \beta_2) = (0.5, 0.999)\), in the framework PyTorch [87].

Figure 9 shows the average mismatch of the surrogate as a function of the latent space dimension \(N_Z\), computed over the \(k\)-Fold cross validation process. Even though the waveforms depend on a single physical parameter, namely the frequency \(\omega/\mu_V\), the model has proven to perform much better if we allow at least 2 dimensions in the latent space. A latent dimension higher than that does not have a large influence on the mismatch. For this reason, the final selected model has ben chosen to have \(N_Z = 2\).

3.4 The Translator Phase↩︎

The task of the translator phase of the model is to interpolate the coordinates in feature space \(F^j\) as a function of the Proca star frequency \(\omega/\mu_V\). As there is only one independent variable, order 3 spline interpolation is an efficient and fast method. These are provided by the scipy.interpolate function interp1d, which is applied to every coordinate in the latent space, as well as the peak position and amplitudes \((t_0, A)\), as a function of the variable \(\omega/\mu_V\).

4 Monolithic Model↩︎

We have, for completeness, trained a monolithic model based on fully connected neural networks. In this case, instead of performing a preliminary dimensional reduction before the final interpolation in \(\omega/\mu_V\), we allow a single neural network to receive \(\omega/\mu_V\) and produce \(\psi_4^i\) directly. In this case, the real and imaginary part of the waveform are concatenated in a single vector of dimension \(2N_\psi = 1024\), as complex values are not trivially handled by common activation functions. Cubic splines are used to learn the dependence of the peak position and amplitudes \((t_0, A)\), as defined in Section 2.2.

The wave-fitting network has an input layer with 1 neuron (receiving \(\omega/\mu_V\)) and successive layers of sizes 64, 256, 512 and 1024 (the output layer). We use GELU activation functions after each layer except the last one.

As in the case of the autoencoder, training is performed on a Nvidia RTX 4090 24GB GPU over 500 iteration epochs divided in minibatches of 16 waveforms, with an Adam optimizer of learning rate of \(10^{-3}\), and \((\beta_1, \beta_2) = (0.5, 0.999)\).

In order to gauge the dependence of the model on hyperparameters, we introduce a number of replicas of the layer of size 64, exclusively on the wave-fitting network. As an example, the network with 2 extra hidden layers would have sizes \((1, 64, 64, 64, 256, 512, 1024)\). In Fig. 10 we observe that the results improve significantly when we add one extra layer, but plateau after that. For this reason, we used one extra layer for the final model.

5 Results↩︎

5.1 Model Evaluation↩︎

Figure 11 shows histograms of the mismatches obtained on the final testing set by the four categories of models. We consider the two dominant modes, namely \((l,m) = (2,2)\) and \((l,m) = (2,0)\). In both cases we can observe that the SVD model seems to work best, followed by EIR, the autoencoder, and finally the monolithic fully-connected network. There is, however, a significant overlap between the distributions.

There are several possible reasons for the better performance of SVD. A lower number of parameters makes the model less prone to overfitting, and therefore more robust. Also, SVD is based on standard linear algebra algorithms, which do not require training by gradient descent, which removes the possible accuracy losses due to early or late stopping of the training process.

The significantly better performance of the SVD model on the testing set, together with its simplicity and the fact that its training is much faster than those based on neural networks, clearly indicate that this type of model is the most convenient for our physical system. Indeed, the parameter estimation Bayesian parameter inferences in Sections 5.2 and 5.3 have been performed with SVD models with 10 basis vectors and \(N_\psi = 3001\), in this case trained with the complete dataset of 59 waveforms.

5.2 Injection Recovery↩︎

We test the performance of our model in parameter inference and model selection tasks. To this end, we inject in zero noise a family of 30 signals from head-on Proca-star mergers of varying \(\omega/\mu_V\), directly obtained through numerical-relativity simulations. We recover these injections using our SVD surrogate model together with the parameter inference code Bilby [15] in its parallelizable version Parallel Bilby [88]. We set uniform priors in the detector-frame total mass and \(\omega/\mu_V\). In addition, we set isotropic priors both on sky-location and source orientation, together with a uniform prior in the polarisation angle and a prior in distance uniform in co-moving volume. These priors match those used for the analyses presented in [45] and [43] restricted equal-mass systems. We choose a three-detector network formed by the two Advanced LIGO detectors and Advanced Virgo equipped with their PSD at the time of the GW190521 event. We perform our analysis directly on the Newman-Penrose scalar \(\psi_4\), as described in [43]. We sample the parameter space using the sampler Dynesty [89] with \(N=1024\) live points.
Figure 12 shows the posterior distributions, together with \(90\%\) credible intervals, for the field frequency \(\omega/\mu_V\), as a function of the true injected value \(\omega_{\rm true}/\mu_V\). The true values are denoted by red circles while maximum likelihood (ML) values (i.e. those producing the best fits) are denoted by yellow stars. As somewhat expected from our match analyses, we find that the posterior distribution always peaks very close to the true value, with the latter being always contained within the \(90\%\) credible intervals. The same is found in most cases for the ML values, which always lay very close to the true parameters found by the sampler. We note that the sampler is not designed to find the actual best-fit parameters but to reconstruct the corresponding posterior probability distributions. Due to this, in some cases, the recovered ML parameters may not match the true ones, even if the waveform model is infinitely accurate. This is the case, for instance, when there are large degeneracies in the parameter space or when the sampling is not extremely aggressive. Nevertheless, since in certain situations parameter recovery can be strongly influenced by Bayesian priors causing shifts in the posteriors, we report the location of the ML points to show that such shifts are indeed sourced by priors and not by systematic errors due to the waveform model.

We find that certain values of \(\omega_{\rm true}\) lead to varying features in the posteriors. First, some cases lead to wide distributions while others lead to significantly sharper ones. This reflects the fact that, in certain regions of the parameter space, small changes in \(\omega_{\rm true}\) cause significant changes in the waveform morphology while smoother changes are induced in other regions. Second, for some cases we find somewhat bi-modal distributions. Far from denoting any issues in the analysis or systematic errors in the waveform model, these features come as a consequence of the physics governing Proca-star waveforms. In particular, mergers corresponding to the lower and upper end of our \(\omega/\mu_V\) range, are significantly louder than those corresponding to intermediate frequencies. Therefore, such configurations are highly favoured by our distance prior, producing secondary peaks of the posterior distribution either at large frequencies (for the lowest frequency injections), both large and low frequencies (for the mid-frequency injections), and almost imperceptible peaks at low frequencies for the highest-frequency injections. For the latter, the posterior also somewhat extends artificially to lower frequencies due to the upper frequency cutoff imposed by the model, which prevents the posterior distribution to extend beyond such cutoff.

Figure 13 shows the posterior distributions for the total mass \(M_{\rm total}\) of the binary and the corresponding boson mass \(\mu_V\). First, we note that the true values are again always contained within the \(90\%\) credible intervals. Second, we note that the posteriors for \(M_{\rm total}\) are systematically shifted towards lower values. Again, far from reflecting systematic errors in our model, this is due to the impact of the prior for the luminosity distance [43], [45]. Recall that the source-frame total mass is \(M_{\rm total}=M_{\rm det}/(1+z)\), where \(M_{\rm det}\) and \(z\) are, respectively, the detector-frame total mass and the redshift. Since the distance prior strongly favours large distances, we find that such estimates, and therefore those for \(z\), are always shifted towards large values, causing the opposite effect on \(M_{\rm total}\). The same effect is present but less pronounced for \(\mu_V\), owing to its strong dependency on the – better captured – value of \(\omega/\mu_V\). Finally, we note that ML values are always closer (when not directly on top) to the true injected values than the bulk of the posterior distributions.

5.3 GW190521 Parameter Estimation↩︎

We perform parameter inference on the GW event GW190521, which has been previously analysed directly using numerical relativity simulations for head-on Proca-star waveforms, both restricted to equal-mass cases [43], [45] and extended to unequal-mass ones [44]. We set the exact same priors described in the previous section, which match those used in [43], [45]. We will compare our results to those reported in the second column of Table III in [43], which unlike those in [45], make direct use of the Newman-Penrose scalar. The two main differences between these two analyses are: a) The usage of a surrogate model instead of numerical relativity waveforms, and b) the continuous sampling of the parameter space as opposed to the exploration of a discrete set of values of \(\omega/\mu_V\).

Table 1 reports the parameter estimates obtained from our analysis together with those obtained in [43], expressed as median values together with symmetric \(90\%\) credible intervals. In addition, Figure 14 shows the corresponding posterior distributions for the field frequency \(\omega/\mu_V\) and the ultralight-boson mass \(\mu_V\). We note that the results are in very good agreement, indicating again the reliability of our surrogate model. Table 1 also shows the corresponding natural log-Bayes’ Factors with respect to the noise hypothesis, which match within the sampling uncertainty of 0.1. Finally, we find that the numerical relativity model recovers a slightly larger log likelihood, indicating that the best-fitting waveform fits the data marginally better than in the Surrogate case. We attribute this to the following fact. While for the case of the surrogate model we perform a single inference run over all waveform parameters, for the NR model we combine the result of many runs that sample over the extrinsic waveform parameters and total mass for fixed \(\omega/\mu_{V}\) [45], which leads to a more aggressive sampling of the parameter space.

6 Discussion↩︎

We have constructed several accurate surrogate models that generate waveforms of head-on collisions of equal-mass Proca stars. The models are purely data-driven, without any predefined physical input, and are trained on a catalog of 59 numerical relativity waveforms with several multipoles each. Different model architectures have been analyzed. On the one hand, a collection of two-stage models, which implement a dimensional reduction algorithm followed by cubic spline interpolation of the governing parameters. The dimensional reduction algorithms discussed comprise singular value decomposition, empirical interpolant representation, and finally a deep convolutional autoencoder. On the other hand, we have also tested a monolithic model where a fully-connected neural network interpolates the full waveform time series as a function of the physical parameter \(\omega/\mu_V\).

While all four models have proven to be accurate, with all mismatches below \(10^{-3}\) for the mode \((l,m) = (2,2)\), the model based on singular value decomposition has proven to be the most precise. In terms of precision, SVD is followed by the empirical interpolant representation, the autoencoder and finally the monolithic fully-connected network. The two first models are not only the most precise but also much faster than those based on neural networks, as they do not require gradient descent optimization. However, the neural networks are able to encode waveforms nonlinearly from very few input real parameters (2 for the autoencoder, 1 for the fully-connected), while SVD and EIR require at least 7 complex parameters for our dataset.

We have tested the SVD surrogate model within the task of Bayesian parameter inference through the recovery of numerical relativity injections in zero noise, obtaining excellent results. Finally, we have also performed parameter estimation on GW190521 data, finding strong agreement with the results obtained directly using numerical relativity templates. As a general summary, we conclude that the type of surrogate models investigated in this work can be regarded as a convenient tool to efficiently perform Bayesian analysis on gravitational waveforms from bosonic star mergers.

A natural extension of this work is the inclusion of waveforms from Proca star mergers with unequal masses. Our early attempts at this extension, however, have encountered serious difficulties due to the interference effects between the stars at merger [42]. When the Proca fields from the two colliding stars oscillate at different frequencies \(\omega_1\) and \(\omega_2\), the relative phase at the point of encounter depends on the particular values of these frequencies. This causes an interference pattern, where the amplitude of the emitted gravitational waves presents strong variations. Additionally, at the regions of lowest amplitude, the shape of the waveform changes very rapidly with \(\omega_1\) and \(\omega_2\). At present, our catalog of simulations is not sufficiently dense in these difficult regions for the model to reach the desired precision.

References↩︎