Effective-one-body waveforms for binary neutron stars using surrogate models

Gravitational-wave observations of binary neutron star systems can provide information about the masses, spins, and structure of neutron stars. However, this requires accurate and computationally efficient waveform models that take $\lesssim 1\mathrm{s}$ to evaluate for use in Bayesian parameter estimation codes that perform $10^7–10^8$ waveform evaluations. We present a surrogate model of a nonspinning effective-one-body waveform model with $\ell =2$, 3, and 4 tidal multipole moments that reproduces waveforms of binary neutron star numerical simulations up to merger. The surrogate is built from compact sets of effective-one-body waveform amplitude and phase data that each form a reduced basis. We find that 12 amplitude and 7 phase basis elements are sufficient to reconstruct any binary neutron star waveform with a starting frequency of 10 Hz. The surrogate has maximum errors of 3.8% in amplitude (0.04% excluding the last $100M$ before merger) and 0.043 rad in phase. This leads to typical mismatches of ${10}^{-5}-{10}^{-4}$ for Advanced LIGO depending on the component masses, with a worst case match of $7\times {10}^{-4}$ when both stars have masses $\ge 2M_{\odot }$. The version implemented in the LIGO Algorithm Library takes $\sim 0.07\mathrm{s}$ to evaluate for a starting frequency of 30 Hz and $\sim 0.8\mathrm{s}$ for a starting frequency of 10 Hz, resulting in a speed-up factor of $\mathcal{O}(10^3)$ relative to the original matlab code. This allows parameter estimation codes to run in days to weeks rather than years, and we demonstrate this with a nested sampling run that recovers the masses and tidal parameters of a simulated binary neutron star system.

Figure 1

Contribution of each tidal multipole moment to the phase evolution of an equal mass BNS system with component masses $(M_A,M_B)=(1.4,1.4)M_{\odot }$ beginning at 30 Hz for the soft EOS SLY (top) and the stiff EOS MS1b (bottom). The phase contribution is given by the difference in phase between waveforms with no tidal interactions ${\mathrm{\Phi }}_{\text{not}}(t)$, only the $\ell =2$ interaction ${\mathrm{\Phi }}_{\ell =2}(t)$, the $\ell =2$, 3 interactions ${\mathrm{\Phi }}_{\ell =2,3}(t)$, and the $\ell =2$, 3, 4 interactions ${\mathrm{\Phi }}_{\ell =2,3,4}(t)$. Also shown by the dashed curve is the error that results from using the fitting functions ${\mathrm{\Lambda }}_3^{\mathrm{fit}}({\mathrm{\Lambda }}_2)$ and ${\mathrm{\Lambda }}_4^{\mathrm{fit}}({\mathrm{\Lambda }}_2)$ instead of the values of ${\mathrm{\Lambda }}_3$ and ${\mathrm{\Lambda }}_4$ calculated from the EOS ${\mathrm{\Phi }}_{\mathrm{fit}}(t)$. Each curve is plotted as a parametric function of the phase difference between the two waveforms $|{\mathrm{\Phi }}_2(t)-{\mathrm{\Phi }}_1(t)|$ versus the frequency of the waveform with no tidal interactions $f_{\text{not}}(t)$. In this way, the phase difference between waveforms is calculated at the same time instead of the same frequency. This can be more directly compared to phase errors in the surrogate model below which are calculated as a function of time.

Figure 2

The training set is constructed from the Chebyshev-Gauss-Lobatto nodes with 16 nodes in each dimension for a total of 4096 waveforms. The waveform parameters have the range $q\in [0.5,1]$, ${\mathrm{\Lambda }}_2^A\in [50,5000]$, and ${\mathrm{\Lambda }}_2^B\in [50,5000]$. The same grid is also used for the Chebyshev interpolation to evaluate the amplitude and phase at the empirical nodes ${\tau }_j$. Red $△\mathrm{s}$ represent the 12 waveforms chosen by the greedy algorithm to generate the amplitude-reduced basis, while blue $○\mathrm{s}$ represent the 7 waveforms chosen for the phase-reduced basis.

Figure 3

The seven orthonormal reduced basis functions $\{e_i^{\mathrm{\Phi }}(t){\}}_{i=1}^7$ for the waveform phase. These functions accurately capture the features in the waveform phases within the range of parameters considered. The $x$ axis is linear in the range $[-10^1,10^1]$ and logarithmic elsewhere.

Figure 4

The seven elements of the empirical interpolation operator $\{B_j^{\mathrm{\Phi }}(t){\}}_{j=1}^7$ for the waveform phase. The empirical interpolation nodes $\{{\tau }_j^{\mathrm{\Phi }}{\}}_{j=1}^7$ are shown as black dots. The $x$ axis is linear in the range $[-10^1,10^1]$ and logarithmic elsewhere.

Figure 5

Error between the surrogate and the $16^3$ training-set waveforms used to construct the reduced basis. Larger points represent larger errors. The fractional amplitude and phase errors, maximized over time and waveform parameters, are $\mathrm{\Delta }A/A=7.7\times {10}^{-5}$ and $\mathrm{\Delta }\mathrm{\Phi }=0.014$, respectively. These errors are due to the finite number of reduced bases.

Figure 6

Error between the surrogate and $10^4$ waveforms with randomly sampled parameters not in the training set. Larger points represent larger errors. Top: Fractional amplitude error $\mathrm{\Delta }A/A$ maximized over all times except the last $100M$. Middle: Fractional amplitude error $\mathrm{\Delta }A/A$ maximized over all times. Bottom: Phase error $\mathrm{\Delta }\mathrm{\Phi }$ maximized over all times.

Figure 7

Final $\sim 10$ cycles of the waveform with the largest phase error. The parameters are $\{q,{\mathrm{\Lambda }}_2^A,{\mathrm{\Lambda }}_2^B\}=\{0.618,420,81\}$. The amplitude errors are always largest during the last gravitational-wave cycle, and are about 2 orders of magnitude smaller before the last gravitational-wave cycle.

Figure 8

Histogram of the mismatch between the surrogate and the $10^4$ randomly sampled EOB waveforms. The smaller mass $M_B$ is set to either $1M_{\odot }$ or $2M_{\odot }$. The PSD corresponds to the aLIGO design sensitivity and the lower and upper frequency integration bounds are $f_{\mathrm{low}}=30\mathrm{Hz}$ and $f_{\mathrm{high}}=2048\mathrm{Hz}$.

Figure 9

Performance of the surrogate implemented in LAL for a binary with component masses $(1.4,1.4)M_{\odot }$. We used sampling frequencies of 4096 Hz and 16384 Hz. The evaluation time at each starting frequency is averaged over 16 evaluations on a 2.7 GHz Intel Xeon CPU using one core.

Figure 10

Posterior densities for the tidal parameters $\tilde{\mathrm{\Lambda }}$ (top) and $\delta \tilde{\mathrm{\Lambda }}$ (bottom). The simulated signal had component masses $(1.4,1.4)M_{\odot }$ and EOS MS1b, so that $\tilde{\mathrm{\Lambda }}=1286$ and $\delta \tilde{\mathrm{\Lambda }}=0$, as indicated by the vertical red lines.

Figure 11

Relative errors $\mathrm{\Delta }{\mathrm{\Lambda }}_3/{\mathrm{\Lambda }}_3$ and $\mathrm{\Delta }{\mathrm{\Lambda }}_4/{\mathrm{\Lambda }}_4$, with $\mathrm{\Delta }{\mathrm{\Lambda }}_{\ell }={\mathrm{\Lambda }}_{\ell }-{\mathrm{\Lambda }}_{\ell }^{\mathrm{fit}}$ (only 300 points shown).