The gravitational waves emitted by an EMRI are among the most complicated compact binary signals due to their long duration and rich harmonic content. This complexity stems from the dynamics of the secondary object inspiralling into the central MBH. Such orbits are characterized by three degrees of freedom associated with the dimensionless fundamental frequencies of Kerr geodesic orbits Ω_{r, θ, ϕ} [42, 43]. These quantities are determined by the dimensionless spin magnitude of the MBH a and the quasi-Keplerian orbital parameters of p (semi-latus rectum or separation in units of M), e (eccentricity), and x_I ≡ cosI (cosine of the angle I which describes the inclination of the orbit from the equatorial plane). At adiabatic order, for a given system with primary mass M, secondary mass μ, initial orbital parameters (p₀, e₀, x_I0), and initial phases Φ_φ0, Φ_θ0, Φ_r0, the orbital evolution can be determined by solving the following system of ordinary differential equations:

The orbital-element fluxes f_{p, e,xI} account for the gravitational wave emission and are interpolated from pre-computed grids. In the FEW package, the system of ordinary differential equations is solved using an explicit embedded Runge-Kutta Prince-Dormand (8, 9) method with adaptive step size [see gsl_odeiv2_step_rk8pd in Galassi [44]], resulting in sparse orbital trajectories of ~ 100 data points in length.

Once the trajectory is fully specified by {p(t), e(t), x_I(t), Φ_φ(t), Φ_θ(t), Φ_r(t)}, we compute the complex gravitational-wave amplitudes A_lmkn(p(t), e(t), x_I(t)). Each amplitude mode is associated with the multipole l of the radiation and the harmonic {m, n, k} of the fundamental orbital frequencies, where m, k, and n label the azimuthal, polar, and radial modes, respectively. In the FEW package, the amplitudes and phases are interpolated using a “not-a-knot” cubic spline to handle the large number of modes using a sparse number of points. Cubic splines are the simplest polynomial basis spline suitable for our purposes because they provide an approximation of the second derivative of frequency with respect to time. Higher order spline functions are more expensive and might overfit data as sparsely evaluated as our trajectory information. Additionally, the root-finding procedure that will be used in constructing the time-frequency correspondence in Eq. 6 is significantly more expensive for higher order polynomials due to the need to perform this task numerically. In the interest of computational efficiency, we control the number of harmonic modes with a further mode-selection step. The mode-selection step is defined by the threshold ϵ of the total power emitted by all modes. Within each time step of the sparse array, we take |Almkn(t)|2 and sort the individual modes in descending order and compute a cumulative summation. From this sorted array, we check whether the additional power contributed to an individual mode falls below the threshold ϵ related to the total power emitted. This gives a set of contributing modes within each time step. The selected modes are given by the union of all contributing modes across all time steps. This ensures to maintain continuity across time. Only the modes that pass this threshold are included in the waveform computation [see Katz et al. [39] for further details].

The gravitational wave signal at a large distance from the source can be represented in terms of the complex time-domain dimensionless strain [40]:

where t is the time of arrival of the gravitational wave at the Solar System barycenter, θ is the source-frame polar viewing angle, ϕ is the source-frame azimuthal viewing angle, and d_L is the luminosity distance of the source. The functions S_lmkn(t, θ) are spin-weighted spheroidal harmonics. For inspiral into a Schwarzschild black hole (the focus of this analysis), these functions are identical to the spin-weighted spherical harmonics. For the general Kerr case, they differ from the spherical harmonics (though they can be very usefully expanded into spherical harmonics [45]), and in particular, they depend on the orbital frequencies. Because the orbital frequencies evolve with time, the spheroidal harmonics evolve as well [40], complicating the description of the waveform. Future studies will describe how we handle this complication. In any case, the usual transverse-traceless gravitational wave polarizations can be obtained from the real and imaginary parts of the complex conjugate of h.

A key element to be noted in Eq. 1 is that the EMRI waveform can be viewed as a sum of simple “voices”, where each voice corresponds to a mode (l, m, k, n). The time-domain waveform can thus be rewritten in terms of a sum over voices V:

where each voice is characterized by a complex amplitude H_V(t) and a phase Φ_V(t). The voice-by-voice decomposition was suggested to one of the present authors by L. S. Finn and first presented in Hughes [46].¹ In this study, we provide the first ready-to-use implementation of this decomposition for frequency-domain EMRI waveforms.

EMRI waveforms can be represented in the frequency domain by using the stationary phase approximation. Because amplitude, phase, and frequency in an EMRI evolve slowly (as long as the two-timescale approximation is valid [48]), we expect the stationary phase approximation (SPA) to provide a high-quality approximation to the Fourier transform of the signal [49]. Here, we review the computation of the Fourier transform of an EMRI waveform in the frequency domain presented in the study by Hughes et al. [40]. We begin by taking the Fourier transform of each voice of Eq. 2 as follows:

where we have not written the subscript for ease of notation. To compute the stationary phase approximation to the signal, it is necessary to expand the phase evolution, written as

where

We have introduced the instantaneous frequency of the signal and the instantaneous first and second derivatives of the frequency at the stationary time t = t_S for which F = f. We expanded up to third order in the phase evolution because for some voices, the frequency evolution is not monotonic. If the frequency derivative vanishes (Ḟ = 0), the standard SPA Fourier transform is singular and cannot be used [50]. This is highly relevant to EMRI signals because the frequency associated with many voices rises to a maximum and then decreases. In particular, this occurs for EMRI voices which involve harmonics of the radial frequency. This is because Ω_r becomes very small (approaching |Ω_ϕ|) in the weak-field limit and goes to zero as systems approach the last stable orbit; a maximum of Ω_r exists between these two “small” limits.

We now use Eq. 4 to obtain the time-frequency correspondence for each voice, i.e., we find the time t such that F(t) = f. The frequency evolution of each voice F(t) = (mΩ_φ(t)+kΩ_θ(t)+nΩ_r(t))/(2πM) is directly obtained from the trajectory evolution of each of the fundamental frequencies. During the waveform production process, the fundamental frequency evolution is interpolated with a cubic spline, allowing for a good approximation using a third-order polynomial. Doing so, the condition f = F(t) between two knots can be rewritten as:

where A, B, C, and D are real numbers given by the sum of the spline coefficients of three frequencies multiplied by their mode number, and t_i is the time at the beginning of the spline segment.

The fundamental frequency, fundamental phase, and amplitude spline coefficients are determined along the sparse trajectory prior to the final waveform summation phase [see Katz et al. [39] for more information]. Additionally, using all values of f and each set of frequencies found along the sparse trajectory for each voice, the spline segment, in which each value of f falls, is predetermined before the summation. Once in the summation, the goal is to determine t(f) as efficiently as possible. This is difficult in cases where the voice is not monotonic because t(f) is not a bijective mapping for this case: two values of t can be represented by the same value of f. To address this issue, we solve for the cubic roots of Eq. 6 within each segment, where f is found (one or two segments per voice depending on f). We take the real-valued roots and determine which one or two actually exist between the bounds of the segment, and we are currently examining (it can be two only for the segment in which the turnover occurs). Since the frequency evolution is slowly varying within each segment, the condition Ḟ = 0 cannot occur more than one time per segment. These time values occurring within our segment are then used to determine all necessary quantities for waveform building because all informative splines are fit as a function of t.

We define t_j(f) to be the jth time (first or second) at which F(t) = f, and write Ḟ_j ≡ Ḟ[t_j(f)], F¨j≡F¨[tj(f)]. These higher order derivatives are also computed using the fundamental frequency spline information. With that, for each value of f for each voice, the Fourier transform can be represented as:

where N ∈ {1, 2} is the number of time values associated with each value of f. To perform this integral, we set α = γ + 2πiḞ, with γ real and positive, define β=2πF¨, and use

where K_n(z) is the modified Bessel function of the second type. Taking the limit γ → 0, we find

This result defines our frequency domain waveform. Expanding about “small” F¨J, it can be shown that this yields the standard SPA when F¨j→0 and N = 1 [see Hughes et al. [40]]. For computational efficiency, we approximate the term ixe-iXK1/3(-iX) with two polynomial expansions in X=2πiḞj3/3F¨j2, valid for |X| < 7 and |X| > 7. The largest relative error occurs around the transition region |X| ≈ 7 and is ≈ 10⁻³. We did not explore the potential of applying this method for calculating other oscillatory integrals [51].

As previously mentioned, this waveform is built into the framework of FEW. All parts of the waveform generation prior to the “waveform summation module” are identical to the original time-domain implementation [39]. The waveform summation described above is implemented for both CPUs and GPUs. For the purposes of GPU parallelization, the summation is grouped by harmonic voice and segment of the splines defined in time. This creates a rectangular computation grid because all spline information occurs on the same sparse time evolution returned from the trajectory integrator yielding a constant number of segments across all spline constituents. On the GPU, we place each segment within each voice on a separate GPU block (larger parallelization unit). Within each block, we use the parallel threads (smaller parallelization unit) to evaluate all values of interest in frequency for that specific voice and segment. The waveform is then constructed by “atomically” adding from each of these blocks and threads to an overall waveform vector stored in GPU global memory.

An EMRI waveform for non-rotating black holes is uniquely determined in the reference frame of the Solar System barycenter by setting the following parameters:

where θ_S and ϕ_S are the polar and azimuthal sky location angles given in the Solar System barycenter reference frame, and θ_K and ϕ_K are the azimuthal and polar angles describing the orientation of the orbital angular momentum. The relation to the source frame angles θ, ϕ can be found in the study by Katz et al. [39]. We do not include the effect of the detector response function, but we leave this to future studies [52].

We check the performance and accuracy of the FD waveform over the EMRI parameter space confined to log(M/M⊙)∈[log105,log107), log(μ/M) ∈ [log10⁻⁶, log10⁻⁴), e₀ ∈ [0.001, 0.7), and Φ_φ0, Φ_r0 ∈ [0, 2π) by uniformly drawing 5,000 parameter realizations. The initial semi-latus rectum is fixed to obtain an inspiral with a time duration of 0.99T, where T is the duration of the observation. We consider h₊ and h_× in the Solar System barycenter as our two channels, and we fix {d_L = 1Gpc, θ_S = π/3, ϕ_S = π/3, θ_K = π/3, ϕ_K = π/3} throughout this study. To check our implementation, we quantify how similar two waveforms are using the mismatch²

where we have introduced the inner product ·· of each channel α,

where f_i = i/T, with i from 1 up to the total number of frequency bins N_f = N/2 with N number of time points. Notably, we will also consider frequency arrays with spacing different from the one previously mentioned. The tilde indicates the Fourier transform, and S_n(f) is the one-sided noise spectral density for the LISA data streams, which we take from the study by Babak et al. [53]. We use the mismatch M(hFD,hTD) to compare the frequency domain implementation with the DFT of the time-domain implementation sampled with a time interval Δt = T/N. The conversion from continuous to discrete samples changes the underlying Fourier transform into the discrete Fourier transform. This leads to a type of distortion called aliasing. The choice of an appropriate time sampling Δt is the key to minimizing such distortion. The Nyquist rate is the maximum resolvable frequency for the DFT of a discrete-time signal, and its value equals 1/(2Δt). When the signal is digitized, we need to make sure that the Nyquist frequency is larger than the largest harmonic frequency F_V(t) of the signal, i.e. max_VF(t) < 1/(2Δt), where V = (l, m, n, k). By ensuring this condition, the resulting discrete-time sequence is free of aliasing. Similarly, the conversion from a very long (or infinite) sequence to a manageable size entails a type of distortion called leakage, which is manifested as a loss of resolution in the DFT. We mitigate this effect by applying the Hann window to the time domain waveform [54]. To account for the effect of windowing in the FD waveform, we perform a convolution of the FD signal with the DFT of the window function. We do not maximize over the phase since the initial conditions of the frequency and phase evolution of the TD and FD waveforms are the same by construction. Global relative time and phase shifts are not observable parameters, so by not optimizing over them, we are being conservative with the model and folding data processing into the model comparison.

We show the impact of windowing and spectral leakage in Figure 1 for a system with μ = 10M_⊙, M=106M⊙, e₀ = 0.6, p₀ = 9.56 for an observation of T = 1 year, a mode content threshold ϵ = 10⁻², and a sampling interval of Δt = 10 s. We show the amplitude squared, |h~+(f)|2, of the frequency domain waveform (dashed orange line) and DFT of the time domain waveform (solid blue line) for the plus polarization in the upper panel of Figure 1. For reference, we also show the LISA sensitivity curve used to calculate the inner products. In the upper panel of Figure 1, we do not apply any windowing, and therefore, the FD waveform is zero only where there is no mode content. Conversely, the DFT of the TD waveform is affected by spectral leakage, especially for frequencies larger than 10⁻²Hz. If we apply the Hann window consistently to both the FD and TD signals as previously described, we obtain the spectrum shown in the lower panel of Figure 1. The effect of windowing reduces the total power but improves the agreement between the waveforms such that the mismatch decreases from 8 × 10⁻⁴ to 4 × 10⁻⁶.

It is often quoted that the average mismatch for a waveform needs to be 1 − 0.9^1/3 ≈ 1 − 0.965 to recover 90% of the signals in an experiment [55, 56]. However, the requirements on the mismatch for parameter estimation (PE) studies are not as straightforward to characterize, primarily as there is no one-to-one relationship between the level of mismatch and the level of bias that is incurred during inference. A sensible requirement for a PE study is to ask that the systematic (or mismodeling) error coming from an approximate waveform is lower than the statistical error coming from the noise in the data. While the former is independent of the signal-to-noise ratio (SNR=〈h+h+〉+〈h×h×〉), the latter depends on the SNR, and therefore the mismatch requirement on a waveform for a PE study is SNR-dependent. The mismatch times the SNR² is approximately half the model waveform error: M(h1,h2)×SNR2≈〈h1-h2h1-h2〉/2 [57]. Therefore, we perform PE on the system that yields the largest value of M×SNR2.

In particular, we check that the posterior distribution obtained with an FD waveform is equivalent to the one obtained with a TD waveform. To do this, we use the likelihood given by [58]:

where we ignored the normalization factor that needs to be included when fitting for the shape of the noise spectral density S_n. We do not inject any noise in the datastream since we want to check whether any waveform systematic arises due to the difference between the FD and TD implementations. Notably, if we apply a window to our data, we should take into account its effects in the likelihood. In fact, the expectation value of the inner product of the noise n(t) = s(t)−h(t) is different if a window is applied

We do not take into account the effect of windowing in the likelihood because this would lead to a non-diagonal noise covariance matrix and an increase in the computational cost of the inner product. This is not a problem for the purpose of our study because we are consistently neglecting the effect of windowing in the likelihood when analyzing the data with different waveform models (FD or TD) [59]. We remind the reader that we apply the windowing solely because the two waveforms (FD and TD) are in two different domains and we perform the parameter estimation only to check the presence of waveform systematics. However, we remark that the effect of windowing (and gaps) in LISA data is an important issue that remains to be addressed, but it is beyond the scope of this study. To deal with the data gaps expected in LISA, it would be worth investigating the usage of non-uniform DFT.

We calculate the mismatch between the frequency and time domain waveforms as described in Section 2.3. The mismatches obtained over the parameter space are shown in the upper part of Figure 2. The median of the mismatch distributions are between 4 × 10⁻⁶ and 10⁻⁵, whereas the largest 95% quantile of the mismatch is 9.7 × 10⁻⁵.

The distributions obtained with a different mode content ϵ = 10⁻⁵ (solid blue histogram) and ϵ = 10⁻² (dashed green histogram) are almost identical. This is expected since these systems are affected in the same way by spectral leakage. When comparing the mismatches obtained with different sampling intervals Δt = 5 s (solid blue histogram) and Δt = 10 s (dotted orange histogram), we obtain lower mismatches for larger sampling intervals. We attribute this difference to the fact that the total number of points T/Δt contributing to the mismatch is different. For the systems with lower T/Δt, the spectral leakage is weaker since the total power is lower. To confirm this, we take the ratio of the mismatches MΔt=10s/MΔt=5s, and we find it to be approximately constant as a function of the mass. For the same reason, the mismatches obtained for T = 2 years (dashed-dotted red histogram) are slightly lower than the ones obtain for the fiducial value T = 4 years (solid blue histogram). However, the mismatches obtained for (T, Δt) = (2yrs, 5s) (dashed-dotted red histogram) are slightly larger than the ones obtained for (T, Δt) = (4yrs, 10s)(dotted orange histogram). This is probably due to the fact that we fix p₀ for a given inspiral duration, and the maximum p₀ allowed by the current implementation is p₀≈16.

In the lower panel of Figure 2, we show the behavior of the mismatch as a function of the central black hole mass M. The correlation coefficients between the mass M and the mismatch M are approximately ~0.6. This is caused by the shape of the LISA sensitivity curve and the fact that for larger M systems, the spectrum shifts to lower frequencies, and the spectral leakage present at the high frequencies becomes more dominant. For large M, the fiducial distribution of mismatches (blue crosses) is wider than the one obtained with Δt = 10 s, but both distributions reach approximately the same mismatches of ~10⁻⁵ for M~10⁷. This is no longer the case when M becomes smaller. For M~10⁵, the mismatches of Δt = 10 s are generally smaller than the ones obtained from Δt = 5 s.³ This confirms that for larger masses, the spectral leakage dominates, whereas for smaller masses, the mismatch is more affected by the total power in the signal. This behavior is also confirmed for systems with T = 2 years (red crosses). In fact, if we compare the systems with (T, Δt) = (4yrs, 10s) and (T, Δt) = (2yrs, 5s) that have the same power, the distribution for the latter system reaches larger values than the distribution for the former system at all masses because the maximum resolvable frequency, 1/(2Δt), is higher. We note that for lower masses, the distribution for (T, Δt) = (2yrs, 5s) has lower mismatches than the distribution for (T, Δt) = (4yrs, 5s) because the total power T/Δt is smaller.

We use Bayesian parameter estimation to assess whether the FD waveforms are affected by systematic errors. We use the Eryn package [60] to perform a Markov chain Monte Carlo analysis to estimate the posterior distribution of the EMRI intrinsic parameters. Since this analysis is computationally expensive and cannot be performed over the 5000 realizations considered in the mismatch analysis, we decided to focus on the EMRI parameters that yield the largest M×SNR2. For the values of T = 4 yrs, Δt = 5 s or Δt = 10 s, ϵ = 10⁻⁵ (blue solid and dotted orange histograms), we find that the largest value of M×SNR2 is reached by the system with parameters shown in Table 1. The posterior distribution of this system is shown in Figure 3 for the three scenarios in which the FD and TD are used as injections or templates. The Kullback–Leibler divergences [61] between each pair of posteriors are of order 10⁻⁴ for all the considered parameters. Therefore, we conclude that we do not find any significant loss of accuracy or precision for the worst point in parameter space.

The computational cost of the time and frequency domain implementations differ only in the final waveform summation and output domain. Here, we compare the waveform speeds without including the DFT cost. We make this choice to highlight the difference in the waveform generation cost, and we provide the speed of the DFT for a reference system. The speed of the waveform is evaluated using the NVIDIA A100 GPU and the computer cluster Hypatia.⁴ for GPU and CPU timings reported below. We warn the reader that the computational cost of the waveform can vary depending on the computing resources used.

By drawing 5,000 realizations from the parameter space as discussed in the previous section, we show the speed-up factor given by the ratio of the TD and FD waveform evaluation speeds per realization. The results for the GPU speed-up factor are shown in Figure 4. The upper panel shows the speed-up as a function of sampling interval, observation time, and mode content, while the lower panel shows the dependence on eccentricity and central mass for the fiducial system. Notably the median speeds of all considered configurations are 0.044 and 0.055 s for the FD and TD waveforms, respectively. The distributions shown in the upper panel of Figure 4 depend on the priors chosen for the parameter space.

Since the mode threshold ϵ controls the mode content, we can deduce that the FD waveform generation is faster than the TD one when the harmonic content is larger, i.e., ϵ is smaller (see solid blue and dashed green histograms). This is a key result because future expansion of the FEW package will implement EMRI systems with a central rotating black hole, where the number of harmonics is expected to increase by an order of magnitude.

As shown by the solid blue and red dashed-dotted histograms in the upper panel of Figure 4, the FD waveform creation is faster than the TD one for longer signals. This is expected since the TD generation increases with the number of time points. Instead, the FD generation slows down with an increasing frequency resolution. This is confirmed by comparing the orange and blue histograms, which differ only by the time sampling interval Δt. The FD generation for our fiducial case is two times faster for systems with M > 2 × 10⁶ and e₀ > 0.2. The performance of the FD waveform compared with the TD one improves with increasing eccentricity and central mass M. However, the speed-up trend shown in the lower panel of Figure 4 moves to higher eccentricities and higher masses when T = 2 years, resulting in a smaller speed-up region. In the limit of short signal durations, we expect the TD generation to be faster than the FD one, and contrary in the long duration regime.

The CPU waveform generation can take up to hundreds of seconds depending on (T, ϵ, Δt). Due to the significantly higher computational cost, we do not present the full scan of the parameter space, but from 10 parameter realizations, we conclude that for (T, ϵ, Δt) = (4yrs, 10⁻², 10s), the FD waveform generation is on average five times faster than the TD one on CPUs and takes 5 s in median. For reference, we also report that in Table 2, the main findings after timing the FD and TD speeds for the system of Table 1 with (T, ϵ, Δt) = (4yrs, 10⁻⁵, 10s). The FD waveform is seven times faster than the TD waveform on CPUs but only 1.3 times faster on GPUs. The DFT operation is particularly expensive on CPUs taking 7 s compared with the 0.4 milliseconds on a GPU.

One of the main advantages of the FD generation is that we can evaluate the waveform on a sparser frequency array than the one imposed by the frequency resolution 1/T. We use an evenly-spaced frequency array that spans only the region where the injected waveform is non-zero f ∈ [0, f_max]. We change the frequency spacing to downsample the array, and we define the total number of frequency bins to be N_f = f_max/Δf. We show that in Table 2, the computational cost of the waveform generation in frequency and time domain on GPU and CPU obtained for an EMRI system with the parameters specified in Table 1 and (T, ϵ, Δt) = (4yrs, 10⁻⁵, 10s). As shown in Table 2, the frequency downsampling strongly affects the CPU speed-up factors, reducing the cost by a factor of 34. For GPUs, instead, the frequency downsampling mildly affects the speed-up factor between FD and TD generation since the GPU timing is already of the order of tens of milliseconds.

One of the downsides of reducing the frequency resolution is the loss of accuracy in the inner product calculation, with a possible impact on parameter inference. Here, we study the loss of accuracy as a function of the number of frequency bins. As a proof of concept, we consider the source with the parameters reported in Table 1 and obtain the posterior distributions using the FD waveform with different input frequencies and with (T, ϵ, Δt) = (4yrs, 10⁻², 10s). In this analysis, we did not include windowing since the injected and template waveforms are in the same domain. The absence of windowing also affects the total SNR, so we rescaled the distance to keep the SNR reported in Table 1 for the finest frequency resolution. The input array is given by evenly spaced frequencies between zero and a maximum frequency f_max. The spacing is set by selecting the total number of frequency bins N_f.

In Figure 5, we show the posteriors obtained with different levels of downsampling and report the CPU and GPU likelihood speeds. The solid blue posterior is obtained using a frequency array defined by the DFT resolution f_max = 1/(2Δt) and N_f = 6311631, whereas the maximum frequency of the dashed orange and dash-dotted green posteriors is set by the maximum harmonic frequency present in the signal f_max = 2.8 mHz. The only difference between the three cases is due to the approximation of the innner product in calculating the likelihood. Even using only 0.05% of the frequencies, the downsampled analysis with N_f = 3554 (dashed orange histograms) provides a posterior indistinguishable from the one of the full frequency array (solid blue). When further decreasing N_f to 35, the downsampling starts to affect the obtained posteriors. The frequency downsampling leads to a faster likelihood evaluation which drops from 4.26 s for the full array to 0.34 for CPUs. Further reducing the number of frequency bins does not significantly improve the speed, highlighting that the bottleneck of the likelihood computation might be due to other waveform computation steps, such as mode selection [39]. The GPU likelihood timing is reduced by a factor of two when using downsampling.

The reason for the deviation in the posteriors can be explained by looking at the absolute square of the characteristic strain |fh~(f)|2, as shown in Figure 6. The FD waveform spectrum obtained with N_f = 3554 (dashed orange line) well represents the spectrum obtained with the full array N_f = 6311631 (solid blue line) in the frequency range, where the signal is above the LISA sensitivity. Instead, the spectrum obtained with N_f = 35 (dotted green line) has a very sparse frequency array that cannot describe all the features present in the EMRI signal. We conclude that the usage of frequency downsampling can be an important tool for exploratory studies of EMRI parameter inference on CPUs. This motivates the FD implementation as an alternative to the TD one for CPU facilities. However, we stress that it is essential to check the validity of the downsampling before using it.