Fast prediction and evaluation of eccentric inspirals using reduced-order models
FFast prediction and evaluation of eccentric inspirals using reduced-order models
D´aniel Barta ∗ and M´aty´as Vas´uth † Institute for Particle and Nuclear Physics, Wigner Research Centre for Physics,Hungarian Academy of Sciences, Konkoly-Thege Mikl´os ´ut 29-33., H-1121 Budapest, Hungary (Dated: March 2, 2018)A large number of theoretically predicted waveforms are required by matched-filtering searchesfor the gravitational-wave signals produced by compact binary coalescence. In order to substantiallyalleviate the computational burden in gravitational-wave searches and parameter estimation withoutdegrading the signal detectability, we propose a novel reduced-order-model (ROM) approach withapplications to adiabatic 3PN-accurate inspiral waveforms of sources that evolve on either highlyor slightly eccentric orbits. We provide a singular-value decomposition-based reduced-basis methodin the frequency domain to generate reduced-order approximations of any gravitational waves withacceptable accuracy and precision within the parameter range of the model. We construct effi-cient reduced bases comprised of a relatively small number of the most relevant waveforms over3-dimensional parameter-space covered by the template bank (total mass 2 . M (cid:12) ≤ M ≤ M (cid:12) ,mass ratio 0 . ≤ q ≤
1, and initial orbital eccentricity 0 ≤ e ≤ . PACS numbers: 04.30.-w, 04.30.Nk, 95.30.Sf,
I. INTRODUCTION
This present work is a response to the growing demand for computationally-efficient generation of eccentricwaveform families in gravitational-wave (GW) searches. To the extent of our knowledge, surrogate model buildingfor this particular family of waveforms has not yet been tested. As ROM techniques have proved exceedingly efficientfor other models (such as for aligned-spin BBHs), we thus anticipate similar benefits of extremely large speedups inthe time-consuming process of generating eccentric waveform. Our aim is to demonstrate its exceptional expedienceby design and to offer a novel and practical way to dramatically accelerate parameter estimations.Compact binary coalescences (CBCs) in binary compact objects, such as stellar-mass binary black holes (BBHs)and/or binary neutron stars (BNSs), are among the most promising GWs sources for ground-based GW detectors. [1]Binaries that evolved through typical main sequence evolution [2] are expected to shed their formation eccentricitiesover time due to gravitational radiation reaction. For this reason, isolated compact binaries are commonly assumed tomove on quasicircular orbits by the time they spiral into the sensitive frequency band of terrestrial GW observatories.[3, 4] Some relatively young sources, nevertheless, which had too short time for the gravitational radiation reaction tocompletely circularize their orbits retain some residual eccentricity. [5] Therefore, CBC inspirals with non-negligibleorbital eccentricities are plausible sources. [3] Some results [6, 7] support the qualitative conclusion that neglectingresidual orbital eccentricities (even small ones) in CBCs may seriously deteriorate matched-filter detection perfor-mance. A number of possible astrophysical scenarios and mechanisms allows the formation of observationally relevanteccentric ultracompact binaries (see in [8–10]). Short-period CBCs may form by dynamical capture in dense stellarenviroments, present in both galactic central regions and globular clusters, or by tidal capture of compact object byNSs which is described in great detail in [11–13]. Stable hierarchical triple star-systems may form in globular clusterswhere multi-body interactions are involved. It has been estimated that ∼
30% of binaries formed in systems where theKozai resonance increased the eccentricity of the inner binary will have initial eccentricities e > . ∼ e > .
9, where e denotes the initial eccentricity of the binaryby the time it enters the lower part of the frequency band of detectors. [12] Roughly 0 . −
10 eccentric inspiral events ∗ Electronic address: [email protected] † Electronic address: [email protected] a r X i v : . [ g r- q c ] M a r per year up to redshift z ∼ . e (cid:38) .
002 causes systematic errors that exceed statistical errors in aLIGOmeasurements. [15] Since the phasing of the GW signal is significantly more important for parameter estimation,and eccentricity modify the phasing beginning at 1.5 and 0PN orders, eccentricity corrections to the SPA (stationaryphase approximation) phase have to be included at leading order.Putative frequency modulated GW signals (also known as ‘ chirps ’) from CBC inspirals will be buried in thenoisy data streams of the advanced detectors. Data analysis of targeted search techniques operate by matched-filtering to extract any possible signal from the white Gaussian noise by cross-correlating the discrete-time sequencesof the detector data against a large set of theoretical waveform templates (or filters ) which approximate potentialastrophysical signals. [16] The utility of this technique rests partly on how accurately the applied template waveformsmodel the signal being sought. Stellar-mass BBHs and BNSs in the inspiral regime are adequately described by high-order
Post-Newtonian (PN) waveform templates. PN approximations [17] are expansions of Einstein field equationsto any specified order in a small parameter v/c which provide a powerful formalism for modelling CBCs during theinspiral phase, when the orbital speed of the binary v is much smaller than the speed of light c . [14] A PN extensionof order ( v/c ) n to the Newtonian expression of gravity is said to be of ( n/ CBwaves software, developed by the VirgoGroup at Wigner RCP. The 3PN-accurate equations of motion (and the spin precession equations if needed) of theorbiting bodies are integrated by a fourth-order Runge–Kutta method while the far-zone radiation field is determinedby a simultaneous evaluation of analytic waveforms. The waveforms involve all high-order relativistic contributionsfor generic eccentric orbits up to 2PN-order accuracy. [18] The choice of PN template families used in this paper isalso motivated by the very fact that such templates are available in the
LSC Algorithms Library (LAL) which appliessome of them in targeted searches, although current searches do not exceed 2PN order. [20]GWs are parametrized by a set of intrinsic and extrinsic parameters λ = λ extrinsic + λ intrinsic , associated with theastrophysical model of their respective sources. The earlier are intrinsic to the source (such as the masses, spins andeccentricity of the compact objects) while the extrinsic parameters are those which depend on the relative location ofthe source with respect to the detector (such as the time of arrival t of the signal at the detector and the phase of thesignal φ at a reference time t ). Each template has a specific set of values for its parameters which are hereinaftercollectively referred to as model parameters . A collection of points in a p -dimensional parameter space , provided that p is the number of model parameters, is called a template bank (or template grid ). [21] A template bank generatedwith minimal match M M could contain a large number of templates that scales as L ∼ (1 − M M ) − p/ . The numberof templates L required for correlations grows rapidly with p and the number of GW cycles L cyc . [22] A fully coherentGW search for a CBC with p = 8 parameters lasting for L cyc = 10 cycles would require as much as L = 10 waveform model evaluations. [23]Over the last three decades methods have been developed for setting up template banks which minimize thecomputational cost in GW searches without degrading the signal detectability, measured by the signal-to-noise ratio (SNR). [24–26] Since the 1990s a method most feasible for small-dimensional parameter space ( p = 2 ,
3, or 4 atmost) has been popular to address the problem of template placement by associating the parameter space with apositive-definite metric space. In this geometric framework, the metric measures the fractional loss in squared SNR ofa predicted signal (at one point in the parameter space) filtered through the optimal waveform template correspondingto a nearby point in the parameter space. [27] In 2009 a template placement algorithm was developed that is suitablefor any number of dimensions, provided that the metric distance between two points in the parameter space is largeor well-defined. [28]Beside the issue of ensembling sufficiently large template banks, parameter estimation (PE) carries a number ofchallenges unique to large data sets. The exploration of the parameter space of BBHs relies on numerical relativity(NR) simulations of the field equations to discover how such mergers evolve. [29] Even a very coarse survey of theparameter space would require an enormous number, typically L = 10 − [30], of expensive NR simulations whichimpose a computationally insuperable obstacle. The required number is in fact subtantially greater than the combinednumber of all simulations ever performed by each and every NR group [31, 32]. Consequently, techniques which canestimate the astrophysical parameters fast and accurately are needed to overcome this computational bottleneck. [33] Reduced-order modeling or model order reduction is a practical mathematical tool to extract the fundamentalfeatures of a computationally demanding high-order model through exploiting only a reduced set of information.Investigations [34–37] over the last few years have revealed that GW templates exhibit significant redundancy inthe parameter space, suggesting that the amount of information required to represent a fiducial waveform model isappreciably smaller than commonly anticipated. The reduction of information content is achieved through expressingthe essential information by means of only a remarkably few, reduced number of representative waveforms r (cid:28) L toconstruct a reduced-order model (ROM) also known as a surrogate model . ROMs provide compressed approximationsof any selected waveforms within the same physical model. They are projection-based techniques that aim to lowerthe computational complexity in the simulations by mapping the original full-order model (FOM) onto an appropriatesubspace of much lower dimension spanned by a reduced-order basis (RB). To find these representative waveformsthat constitutes the RB several methods, including singular value decompsition (SVD) and greedy methods havebeen proposed, usually combined with the empirical interpolation method (EIM). [34, 38] SVD-based methodshave been applied in Ref. [37, 39, 40] to interpolate time-domain inspiral waveforms. We are going to provide aneffcient (fast and accurate) representations of approximated waveforms for any desired parameter values withinthe model by using the information provided by only r RB waveforms instead of the total number L . [36, 41]The SVD-based approach to significantly accerelate PE process used in Ref. [16] is to directly interpolate thelikelihood function over a significant portion of the parameter space. Moreover there is yet another method, pre-sented in Ref. [29, 42], that defines special reduced-order quadrature (ROQ) rules to assist in fast likelihood evaluation.The rest of the paper is organized as follows: Sec. II deals with the procedure for generating fiducial PN waveformsby CBwaves, with respect to the statistics of the cost of computing individual waveforms to estimate the total costof building template banks. Sec. III proposes the simplest strategy (regular spacing) for template placement inthe intrinsic parameter space, followed by the representation of the fiducial waveform templates on a common, finelysampled and regularly spaced frequency grid. Sec. IV gives a general description of our approach to construct efficientROM assembled from the reduced bases and of its characteristic features, particularly the truncation error. Sec. Vis dedicated to assess the overall performance of the ROM, including the accuracy of the surrogate model and itscomputational cost relative to that of the fiducial model. Conclusions, remarks, limitations and an outlook for futureresearch will be given in Sec. VI. II. FIDUCIAL WAVEFORM MODELS
Current searches for GWs from NS and stellar-mass BH binaries use restricted stationary-phase approximations tothe Fourier transform of 3.5PN-accurate circular inspiral-only waveforms, such as spin-aligned
TaylorF2 or SpinTay-lorT4 . [14] The first part of this section describes a procedure for constructing PN eccentric waveforms by
CBwaves model in the time domain. The second part deals with the statistics of the cost of computing individual time-domain(TD) waveforms, drawn from a relatively large number of sample points in a finite-sample distribution.
A. Construction of eccentric post-Newtonian waveform templates
The
CBwaves open-source software was developed by the Virgo Group at Wigner Research Centre for Physics withthe intent of providing efficient computational tool capable of generating gravitational waveforms produced by genericspinning binary configurations moving on eccentric closed or open orbit within the applied PN framework. A detailedexamination of the software’s performance is given in Ref. [18]. The source release and binary packages supportedboth on x86 and x86 64 platforms are available at the group’s website [19].In the PN formalism the spacetime is assumed to be split into the near and wave zones. The field equations for theperturbed Minkowski metric are solved in both regions. A fourth-order Runge–Kutta (RK4) method with adaptivestep-size control is carred out to numerically solve for the 3PN-accurate near-field radiative dynamics at each time t > t , where t is the time of arrival of the signal at the detector, while the far-zone radiation field [43] decomposedas h ij = 2 Gµc D ( Q ij + P . Q ij + P Q ij + P Q SO ij + P . Q ij + P . Q SO ij + P Q ij + P Q SS ij ) (1)is determined in harmonic coordinates by a simultaneous evaluation of orbital elements ( φ, r, n ) where D is thedistance (typically a few Mpc) to the GW source of that consists of two point particles of masses m and m and µ = m m /M is its reduced mass. The term Q ij is the Newtonian mass quadrupole moment, P . Q ij , P . Q ij , P Q ij are higher-order relativistic corrections up to 2PN order beyond the Newtonian term while P Q SO ij , P . Q SO ij , P Q SS ij are corrections arising from spin–orbit and spin–spin effects, respectively. Here, for brevity, we will not repeat lengthyPN coefficients. They are written out explicitly in the appendix of Ref. [18].Radiative orbital dynamics involving all possible correction terms up the 3PN order beyond the Newtonian termare written out explicitly in terms of mean motion n and orbital eccentricity e in Ref. [10]. The secular evolution istreated adiabatically, assuming that the timescales of the shrinkage of orbits (due to gravitational energy radiation ˙ E )and the precession (due to angular momentum flux ˙ J ) are much longer than that of the orbital period. Consequently,the functions ( ˙ x , ˙ e . . . ) in the equations derived from ˙ E and ˙ J depend only on the eccentricity e , and not oneccentric anomaly u . Hence, the adiabatic evolution equations for x ≡ ( M ω ) / and e form a closed system, and canbe solved independently of the Kepler’s equation. Given initial conditions x (0) and e (0), we can solve the systemof ordinary differential equations numerically to obtain x ( t ) and e ( t ). The integration of the equations of motion isterminated at the innermost stable circular orbit (ISCO), which is located at r ISCO = 6
GM/c (2)in Schwarzschild spacetime (for a non-spinning NS/BH). The orbital angular frequency at the ISCO is f ISCO = c / (6 √ πGM ) , (3)which marks the end of the inspiral phase. Fig. 1 demonstrates that the integration run-time t int depends sensitivelyboth on the initial eccentricity and on the disparity of components’ masses ( m , m ) in a binary system. The t int increases exponentially with decreasing total mass M . The mass disparity, defined by ¯ q ≡ − q , allows bettercomparability with e than q itself, considering that t int asymptotically increases – faster than with decreasing M –towards infinity as either e (left panel) or ¯ q (right panel) tends to 1. The physical interpretation of these competingtrends is very simple:1. The lighter the components of the binary are, the longer it takes for them to gradually descend onto their ISCOthrough a sequence of increasingly circular orbits. [44]2. The more eccentric the orbit was initially, the longer it takes to shed its residual eccentricity over many orbitalperiods. [3]3. Among different configurations of equal total mass, the one with the largest mass disparity has the longestinspiral time for harbouring the lightest component. [44]At the high total-mass region on Fig. 1, the influence of first trend grows comparable to that of the last two to reversethe trend of decreasing integration run-time. Fig. 3 shows the influence of M and q on the length of integrationrun-time t int from a different aspect. Excluding the red and yellow dots, each point in the coloured triangular regionis assigned to a hue level running from dark to light as the value of t int increases on a logarithmic scale. The dark blue‘basin’ represents the region where M and q simultaneously lower the value of t int to its minimum. Isoclines runningin parallel are connecting points at which t int has the same value, therefore they are associated with horizontal linesin Fig. 1 (b). The influence of growing q becoming comparable to and gradually greater than that of M accountsalso for the drift from the linear rising trend in the curvature of isoclines that occurs at the high- q region on Fig. 3.Although Fig. 3 suggests that over 85% of the waveform templates of initial eccentricity e = 0 are computed up to10 seconds, in fact, only 4 .
6% of all waveform templates require less then 10 seconds to integrate, as demonstratedin Fig. 2. Out of a total of 1800, only those 120 templates are shown in Fig. 3 that are located in the e = 0 plane.Still, the figure illustrates well that in the same e -plane the frequency of templates with little t int is extremely highcompared to that of templates with large t int , regardless of the value e .In the next section we shall give a quantitative description of the summary statistics computed from the relativefrequency of occurrence (or empirical probability) of the integration time-runs. B. Probability distribution of integration run-times
Let T = { t int1 , t int2 , . . . , t int L } be a univariate independent and identically distributed (IID) finite data sample drawnfrom the probability (or relative frequency) distribution of the discrete random variable t ∈ T while a discrete set of L time-domain input waveforms h ( t ) ≡ { h ( t ; λ l ) } Ll =1 (4)is computed at each parameter point λ l (see Sec. III A) by evaluating eq. (1) at a distance D = µ simultaneouslywith the integration of the equations of motions at 3PN order that requires integration run-times t int l .Since we do not make any prior assumption about the probability distribution, we shall use a non-parametric modelwhere the statistical measures are determined by the finite data sample T . In statistics, kernel density estimation(KDE) is a fundamental data-smoothing technique that provides a non-parametric estimate, based on observed data T , (a) Integration run-times forequal-mass systems (¯ q = 0). (b) Integration run-times forsystems on circular orbit( e = 0). FIG. 1: The integration run-time t int increases exponentially with decreasing total mass M . With increasing initial eccentricity e (left panel) or mass disparity ¯ q (right panel), t int grows asymptotically at a significantly faster rate than with decreasing M . The integration time of those template waveforms that enter a detector’s sensitivity band at a frequency of 10 Hz hasbeen measured 20 times, each at 11 distinct values of M ∈ [2 . M (cid:12) , M (cid:12) ] for three distinct values of initial eccentricity; e = { , . , . } and mass disparity; ¯ q = { , . , . } represented by blue, orange and green dots, respectively. The templatewaveforms were generated at a uniform sampling frequency 16 .
384 kHz. Around each median curve of corresponding t int values,the shaded bands represent their respective 95% point-wise confidence band. of an unobservable underlying probability density function (PDF) of the continuous random variable inf T ≤ t ≤ sup T .A PDF, denoted by (cid:102) t and illustrated in Fig. 2, is a non-negative Lebesgue-integrable function that defines thecumulative distribution function (CDF) of a real-valued random variable t , evaluated at a value t (cid:48) as (cid:70) t [ t (cid:48) ] ≡ Pr[ t ≤ t (cid:48) ] = (cid:90) t (cid:48) −∞ (cid:102) t [ τ ] dτ. (5)It represents the probability that the random variable t , with the expected value given by E [ t ] = (cid:90) ∞−∞ t (cid:48) d (cid:70) t [ t (cid:48) ] , (6)takes on a value less than or equal to t (cid:48) and its kernel density estimator isˆ (cid:102) h [ t ] = 1 Lh L (cid:88) l =1 K (cid:20) t − t int l h (cid:21) (7)where K ≥ h > h as small as the data sample T allows; however, there is always a trade-offbetween the bias of the estimator and its variance. Another option is the use of adaptive bandwidth kernel estimatorsin which the bandwidth changes as a function of t .A specific quantitative measure of the probability distribution is the n -th moment µ n ≡ E [( t − c ) n ] (8)of the continuous random variable t about some central value c (e.g. the mean, denoted by µ ) where E is the expectedvalue of t defined by eq. (6). The graphical representation of the most common measures of central tendency (mean,median, mode) is depicted on Fig. 2 with solid, dashed and dotted red lines, respectively. The PDF rapidly increaseswith the random variable t up to a point at t = 0 . t = 4 .
438 sec, which marks the mode , i.e. the most frequent value in the distribution.The median which represents the value separating the higher half of the probability distribution from the lower halfis located at t = 20 .
615 sec. The mean which represents the first moment of the PDF ( µ ≡ µ in eq. (8)) is situatedat t = 77 .
499 sec.The central tendency of distributions is typically contrasted with its dispersion that measures the extent to which adistribution stretched or squeezed. Common measures of statistical dispersion are the variance and standard deviation:The variance of t is the second central moment, given by (8) as Var[ t ] ≡ E [( t − µ ) ] and the standard deviation isits square root, denoted by σ . For the given distribution σ = 266 .
885 sec. Finally, the shape (or asymmetry) ofprobability distributions is quantitatively measured by the third and fourth central moments, called skewness and kurtosis and denoted by Skew[ t ] ≡ E [( t − µ ) /σ ] = 18 .
56 and Kurt[ t ] ≡ E [( t − µ ) /σ ] = 490 .
04, respectively. - - - FIG. 2: PDF denoted by (cid:102) t ( t ) (blue line) and CDF by (cid:70) t ( t ) (orange line) are displayed as functions of the random variable t ∈ [inf T , sup T ], corresponding to t int -values, which is measured in seconds on the lower horizontal axis and in standarddeviation ( σ = 266 .
885 sec) round the mean value of t on the upper horizontal axis. The smooth KDE with adaptive bandwidthis based on the data sample T collected from the integration run-times of L = 1800 waveform templates that were generated inthe parameter space Ω, described in Sec. III A. The location of the mode , the median and the arithmetic mean are illustratedby dotted, dashed, and solid red lines, respectively in ascending order of their locations. This order of the measures of centraltendency is a characteristic feature of right-skewed (positive skewness) distributions. III. TEMPLATE PLACEMENT AND THE FREQUENCY GRID
We may now discuss the placement of a grid of TD waveform templates in a compact parameter space, followedby the generation of a sequence of frequency-domain (FD) templates on a common finely sampled uniform frequencygrid. The TD waveforms are Fourier transformed and split into their amplitude and phase parts. These functions areaccuretely represented on a sparse frequency grid with only O (10 ) nodes, with a sampling frequency recorded wellabove the Nyquist frequency of the shortest time-series in the template bank. A. Template placement in an associated 3-dimensional parameter space
The set of input waveforms (4) is computed by
CBwaves , described in Sec. II A, at corresponding parameter points λ ≡ { λ l ∈ Ω } Ll =1 (9)in a compact p -dimensional parameter space Ω ⊂ R p where p is the number of model parameters. We restrictourselves to a feasible 3-dimensional parameter space consisting of totally ordered one-dimenstional sets of values ofcorresponding model parameters ( m , m , e ) that define a sparse grid of points λ ≡ m ⊗ m ⊗ e = { ( m i , m j , e k ) : i ∈ [0 , i max ]; j ∈ [ i, i max ]; k ∈ [0 , k max ] } (10)covering the desired parameter range in the particular model involved. Owing to the invariance of input waveformsunder exchange of the components’ masses ( m , m ), the values of the 2-dimensional index pair ( i, j ) are constrainedto a triangular sub-region in the positive quadrants where i ≤ j . Considering that the waveform templates span a3-dimensional parameter space, each template is successively placed into a single vector (4) as indexed by l ≡ (cid:20)(cid:18) i max − i − (cid:19) i + j (cid:21) k max + k + 1 (11)in the range of values 1 ≤ l ≤ L . This flat index corresponds to the position of templates in the parameter space.The total number of templates in the set is then expressed as L = (cid:104) (2 i max + 1) − (cid:105) k max / . (12)It is desirable to work with a dense grid of short waveforms encompassing the late inspiral phase to make a bettercoverage of the selected region of the parameter space. For the sake of simplicity, we sample at equidistant parametercombinations within the region. Nevertheless, using a template placement algorithm that is based on a template-space metric over the parameter space makes a far more efficient coverage. [45, 46] Generally, the algorithms that usegeometrical techniques concentrate more points near the boundaries of the region and at lower mass-ratios. t int ( sec ) FIG. 3: The template bank h ( t ) of L = 1800 waveform templates was set up over a domain { ( M, q, e ) | M ∈ [2 . M (cid:12) , M (cid:12) ] , q ∈ [0 . , , e ∈ [0 , . } ⊂ Ω by computing eq. (4) with uniform grid spacings { ∆ M = 14 . M (cid:12) , ∆ q =0 . , ∆ e = 0 . } . Red points, confined within a triangular region with a boundary ∂ Ω (thick gray line), represet the pa-rameter points of those 120 input waveforms that are situated in the k = 0 plane section of the parameter space Ω. In orderto measure the accuracy of the ROM of waveforms, eq. (36) is evaluated at equidistant parameter points (yellow points) fromtheir respective nearest basis-waveform neighbors. Each background point in the coloured triangular region is assigned to ahue level running from dark to light as the value of the integration run-time t int increases on a logarithmic scale. t int increasesexponentially with decreasing total mass M and grows asymptotically at a significantly faster rate with increasing mass dis-parity ¯ q . The dark blue ‘basin’, where the great majority of template waveforms are concentrated, represents the region where M and q simultaneously lower the value of t int to its minimum. Isoclines (thin gray curves) running in parallel are connectingpoints at which t int has the same value. A drift from the linear rise in the curvature of isoclines occurs at the high- q region,where the influence of growing q becomes comparable to and gradually greater than that of M . The set of initial eccentricity { e k : 1 ≤ k ≤ k max } is chosen to cover the entire interval [0 , .
95] and the mass ratio q ≡ m /m ≤ q = 1 and relatively extreme systems at q ≈ .
01 withtotal mass
M/M (cid:12) ∈ [2 . , η = µ/M the model covers approximately η ∈ [0 . , . L = 120 templates (red dots) that are situated in the k = 0plane section of the parameter space Ω, out of a total of 1800 templates and are collected in h ( t ). These templatesare confined within a triangular region with a boundary ∂ Ω (thick gray line).
B. Production of frequency-domain waveforms
For optimal orientation all time-domain waveforms in eq. (4), composed of their two fundamental polarizations h + and h × in the dominant (cid:108) = (cid:109) = 2 mode are represented by complex-valued GW strain amplitudes h n ( λ l ) ≡ h + ( t n ; λ l ) − ih × ( t n ; λ l ) (13)at N equidistant grid points { t n = n ∆ t } n ∈ [0 ,N − (14)as elements of a finite sequence of N regularly spaced samples of the complex-valued TD waveforms { h ( λ l ) , h ( λ l ) , . . . , h N − ( λ l ) } . The sequence is converted by a fast Fourier transform (FFT), denoted by a lin-ear operator (cid:70) : h → ˜ h , into an other equivalent-length sequence of regularly spaced samples { ˜ h k ( λ l ) } k ∈ [ − N/ ,N/ − = (cid:70) { h n ( λ l ) } n ∈ [0 ,N − (15)evaluated at the same N equidistant frequency grid points { f − N/ , . . . , f , . . . , f N/ − } considering that ROM con-struction, to be discussed in Sec. IV, will require a set of values that reside in the same grid points over all thewaveforms in the template bank.1. This is achieved by having the length of all frequency series truncated to that of the shortest waveform in time,denoted by T = t N − − t . (16)This particular waveform is associated with the highest mass, lowest eccentricity configuration ( i = i max , j = i max , k = 0) in the template bank and its position in the parameter space, given by eq. (11), is l short =( i max + 3) i max k max / h k ( λ l ) ≡ N − (cid:88) n =0 h n ( λ l ) e − πikn/N , k ∈ [0 , N − , (17)are complex-valued functions of the frequency f k which encodes both the amplitude and the phase,˜ h ( A ) kl = (cid:113) Re[˜ h k ( λ l )] + Im[˜ h k ( λ l )] /N, ˜ h ( φ ) kl = − i ln (cid:16) ˜ h k ( λ l ) / | ˜ h k ( λ l ) | (cid:17) , (18)respectively. In this interpretation, ˜ h k ( λ l ) corresponds to the cross-correlation of the time sequence h n ( λ l ) and an N -periodic complex sinusoid e πikn/N at a frequency point f k ≡ k/N that represents k cycles of the sinusoid. Therefore,eq. (17) acts in place of a matched filter for that frequency. Now, the sequence of frequency-domain waveforms (15)can be re-expressed as ‘chirps’ in a simple form { ˜ h k ( λ l ) } = { ˜ h ( A ) kl exp( i Λ˜ h ( φ ) kl ) } (19)where the oscillation degree Λ is a large number. The behavior of GWs in the late inspiral phase is highly oscillatory,but the amplitude and the phase themselves are smoothly varying functions of frequency. [49] It will thus be moreexpedient to perform high-accuracy parametric fits of the phase and amplitude given by (18) rather than of thecomplex waveform (13) itself. The preprocessed amplitudes and phases are collected in the columns of separatetemplate matrices {H ( A ) , H ( φ ) } ∈ R N × L , H = (˜ h ) kl ∈ C N × L (20)where we have droped the amplitude or phase labels for brevity and where L is the total number of templates, andeach template ˜ h l ( f k ) is given on a common freqency grid of length N . We may choose to represent the waveforms ata large number of frequency points so that N (cid:38) L . C. Definition of a regularly spaced high-resolution frequency grid
Provided that T in (16) is the longest time length, the time spacing is defined as∆ t = T / ( N −
1) (21)by eqs. (14) and (16). The time spacing and the number of time steps N in the grid (14) are chosen such that theFD waveforms (19) are sampled at a rate of f s and cover a suitable and well-resolved frequency range [ f low , f high ].1. The lower limit of the frequency range f low is specified by the low-frequency cutoff of the detector noise spectrumwhich is close f cutoff = 10 Hz for advanced detectors design.2. The upper limit f high is determined to be at f ISCO = 2 .
045 kHz by the ISCO frequency (3) for the lowest totalmass configuration of interest M = 2 . M (cid:12) . - - - - FIG. 4: The panels illustrate the inspiral evolution of equal-mass BBH/BNS systems (¯ q = 0) starting at a Keplerian meanorbital frequency of 5 Hz at a distance D = µ . The four most distinct template waveforms were generated by CBwaves at auniform sampling frequency of 16 .
384 kHz with extreme values of total mass M = { . M (cid:12) , M (cid:12) } and initial eccentricity e = { . , . } in the investigated parameter space Ω. (See large red points on Fig. 3.) The top inset panel presents the last N = 15 ,
000 points of the longest waveform (blue) projected onto an equal number of points of the shortest waveform (red).
The Nyquist criterion requires the sampling frequency to be at least twice the highest frequency contained in thesignal to avoid aliasing. Thus, the smallest sufficient sampling frequency is f s = 4096 samples per second for beingthe first power of 2 to meet the criterion. Note that the typical sampling rate being used by aLIGO and aVirgoobservatories in ongoing searches for GWs is at 2048 Hz. [50] Instead, an equidistant grid with N = 4000 grid pointsis sampled at f s = 16 .
384 kHz in the frequency band
M f ∈ [0 . , . G = c = 1). The totalmass M is expressed in units of geometrized solar mass by M (cid:12) [s] = G/c × M (cid:12) [kg] ≈ . × − sec. At the timeresolution ∆ t = 1 /f s ≈ . M which corresponds to a Nyquist frequency f Ny = f s / ≈ . × − M − , (22)a waveform long enough for the BNS system of total mass M = 2 . M (cid:12) down to f low = 2 . × − M − is givenand is about T = ( N − t ≈ . × M long in time. The spacing in frequency domain is∆ f = 2 f Ny /N, (23)so the power will be either in positive or negative frequencies, depending on conventions and we need to consider onlyhalf of the FFT. Combining this with the relations (21–22), one has∆ f = N − N T ≈ . × − M − . (24)Only half of the points in the FFT spectrum are unique, the rest are symmetrically redundant. Thus, the points ofnegative frequencies contain no new information on the periodicity of the random number sequence. Which amountsto swapping the left and right half of the result of the transform. IV. SVD-BASED REDUCED-ORDER SURROGATE MODEL BUILDING
In this section we summarize some of the characteristic features of SVD that are especially useful for reduced-ordermodeling and discuss our approach to construct a compressed approximate representation of a collection of fiducialwaveforms at the cost of truncation error. Next, projection coefficients of the waveforms are determined in terms ofthe reduced basis. In conclusion, the ROM is assembled from the reduced basis and projection coefficients interpolatedover the parameter space. Our procedure follows the well-established strategy that has been pursued by P¨urrer andby Cannon for building frequency-domain ROMs. [36, 37, 47, 48]0
A. Singular values and truncation error
Formally, the decomposition of the template matrix
H ∈ C N × L in eq. (20) is expressed by a factorization of theform H = V Σ U † , (25)where the complex unitary matrices V = [ v | . . . | v L ] ∈ R L × L ,U = [ u | . . . | u N ] ∈ R N × N (26)are orthogonal sets of non-zero eigenvectors of the non-negative self-adjoint operators H † H and HH † so that U † U = I and V † V = I . The rank–nullity theorem states that the SVD (25) provides a decomposition of the range of H . [47]Accordingly, the left-singular vectors (or eigensamples ) { v i ∈ V } provide an orthonormal basisrange( H ) = span { v , . . . , v R } (27)for the range of H (columnn space) where the maximal number of linearly independent columns of H is R ≡ rank( H ) ≤ L . In a qualitative sense, each v i represents a typical waveform pattern. The right-singular vectors (or eigenfeatures ) { u i ∈ U } provide a basis for the domain of H (row space) and represent the evolution of the magnitude of eachwaveform along the frequency gridpoints. The diagonal entries of the rectangular matrix Σ ∈ R N × L correspond tothe non-negative real singular values (SVs) σ ≥ . . . ≥ σ s ≥ s = min( N, L ). SVs are roots of eigenvalues of H † H (and of HH † ) describing the spectrum of the template matrix H , arranged in monotonically decreasing order(see Fig. 5). If the number of frequency points is significantly larger than the number of waveforms (i.e. L (cid:28) N )then a thin SVD is a more compact and ‘economical’ factorization of eq. (25) than the full-rank
SVD that comprizesall R eigensamples. In practice, low-rank matrices are often contaminated by errors, and for that reason they featurean effective rank R eff smaller than its exact rank R . The reduced-rank approximation of the template matrix H isexpressed by H r = r (cid:88) i =1 σ i v i ⊗ u Ti (28)which comprises only those r < R singular vectors which correspond to singular values of a significant magnitude.The approximated representation (28) of the fiducial template bank H is the r -th partial sum of the outer-productexpansion of the expression (25) where r denotes the desired target rank . The Eckart–Young theorem [51] impliesthat the low-rank SVD in eq. (28) provides the optimal rank- r reconstruction of the template matrix H r ≡ argmin rank( H (cid:48) r )= r (cid:107) H − H (cid:48) r (cid:107) (29)in the least-square sense where the truncation error of approximated representation (28) in both the spectral andFrobenius norm is given by (cid:107) H − H r (cid:107) = σ r +1 ( H ) , (cid:107) H − H r (cid:107) F = (cid:113)(cid:80) min( N,L ) i = r +1 σ i ( H ) , (30)respectively.Fig. 5 shows ˆ σ i ≡ σ i /σ on logarithmic scale as a function of the number of SVD components i = { , . . . , R } involved in the approximated representation. Each ˆ σ i , which describes the relative magnitudes of the correspondingeigenfeatures, is computed from the truncated SVD (28) of template matrices with three distinct full-ranks R = { , , } (i.e. total number of templates). The truncation error in the approximation, in accordance with eq.(30), decreases with the number of SVD components retained. The ultimate accuracy (or minimal error) achievableis limited by the total number of templates L that the original template matrix H contains. The growing rate ofdecay in the SV spectrum demonstrates that the individual SVD components gradually lose their relevance for beingincluded in the approximation. In this respect, the spectrum has three clearly distinctive regions characterized by therate at which SVs decrease:1. Overreduced
SVDs ( k (cid:46) − − − . The initial steep exponential fall attests1that the information contained in the corresponding eigenfeatures is predominantly relevant. In fact, the firstfew components shown on Fig. 6 contain roughly 90% of all the information on the input waveforms, regardlessof rank( H ). Then, SVs decrease at a much lower, yet a slowly increasing rate, practically indistinguishable fordifferent values of full rank R .2. Sufficiently reduced
SVDs (400 (cid:46) k (cid:46) − R is, the moreSVD components have to be kept to achieve the same accuracy of representation.3. Underreduced
SVDs ( k (cid:38) − r is highly dependent on the objective. One either desires a highly accurate recon-struction of the fiducial waveform templates, or a very low dimensional representation of the fundamental features inthe templates. In the former case r should be chosen close to the effective rank , while in the latter case r might bechosen to be much smaller. Fig. 5 demonstrates that choosing a target rank r = 456 for the smallest among fiducialtemplate matrices will result in a truncation error related to ˆ σ = 2 . × − at r = 456. - - - - Overred.Suf. red.
Underred. - - - FIG. 5: Normalized singular-value spectra of the template matrix for full ranks R = { , , } are illustrated by blue,orange and green colour, respectively. The horizontal axis represents the index of SVs, while the vertical axis represents therelative variance of SVs. The main panel displays the relative variance of ˆ σ i of the matrix H ( A ) which encodes the amplitudepart of waveform templates while the corresponding relative variance of ˆ σ i for the phase is shown in the top inset. At r = R − σ r , falls onto a dotted black line given by log ˆ σ r − log σ ≈ − . − . R . The rate at which the ratiodecreases is significantly lower under the dashed black line given by − . − . R . Excluding waveforms in the lowersection causes less errors by a magnitude much smaller than in the upper section. B. Assembly of the surrogate model
The basis for the amplitude or phase space is given in the columns B i of the matrix B ≡ (cid:40) V L ∈ R N × L , if N > LV ∈ R N × N , if N ≤ L (31)and a full-rank basis is desired. If N < L , then the information from L waveforms at N grid points is containedin a basis of dimension N . The reduced basis waveforms only resemble the physical behavior of frequency domain2 - - (a) Basis functions for the amplitude modes. - - (b) Basis functions for the phase modes. FIG. 6: Reduced basis functions for the first 5 amplitude and phase SVD modes are represented at N = 4000 grid points inthe frequency domain. The basis functions become increasingly oscillatory as their index i increases. amplitudes and phases for the first basis function, the higher basis functions are oscillatory (see Fig. 6). To compressthe model, a reduced basis of rank r is selected from the full-rank basis (27) in the form B r = V r = [ v | . . . | v r ] ∈ R N × r for r < R ≤ N. (32)For any r the columns of V r are an optimal orthonormal bases for the starting waveforms. Notice that B r ⊂ B r +1 ,which demonstrates the underlying hierarchical nature of the generated template banks. [47] Fig. 7 may serve as anillustration of the underlying sparsity of the selected basis in the parameter space. The identification of parametervalues associated with the basis waveforms selected by SVD from the full template bank is not that straightforwardas a greedy algorithm would pick values that parameters take. Nevertheless, it may safely be said that a verysmall part of the parameter space volume is covered; the parameter points are heavily concentrated at low-mass andlow-eccentricity values.Hereafter the label r on the rank- r reduced basis will be dropped for brevity. Given the reduced bases B ( A ) and B ( φ ) , we compute projection coefficient vectors (cid:126)µ for any given input waveform ˜ h ∈ R N as follows (cid:126)µ (˜ h ) ≡ B T ˜ h ∈ R r , (33)where the labels referring to amplitude or phase were dropped for brevity. The projection coefficient vectors for allwaveform templates are packed in the matrices M ( A ) and M ( φ ) with entries M kl = µ k (˜ h l ) = ( B T H ) kl ∈ R r × L . (34)Comparing with eq. (25) we see that M = B T H = − Σ U T for a full-rank basis B = V . It follows that the projectioncoefficient matrices are ordered in the same way as the individual waveforms in H . To undo the packing of thewaveforms in the matrices M we just partition the linear index l that enumerates the waveforms in H and obtain atensor M k,l q ,l e = µ k (˜ h ( l q ,l e ) ) ∈ R r × L q × L e . (35)To complete the model we define the projection coefficient vectors at any position in the chosen parameter spaceby suitable interpolants I [ M ]( λ ) ∈ R r for the amplitude and phase coefficient tensors M ( A ) , M ( φ ) . For each inputwaveform we have two corresponding r -vectors of projection coefficients (for amplitude and phase) that are interpolatedover the parameter space. The frequency-domain ROM representation of waveform templates is then constructed inthe form ˜ h S ( λ ; f ) ≡ A ( λ ) I f [ B ( A ) · I [ M ( A ) ]( λ )] × exp { iI f [ B (Φ) · I [ M (Φ) ]( λ )] } , (36)3
50 100 150 200 250 m / M ⊙
50 100150200250 m / M ⊙ e e FIG. 7: The SVD-based reduced-basis parameter choices in the 3-dimensional parameter space ( m , m , e ). Comparing thepositions of the retained r = 600 templates to the placement of the original R = 1800 template shown in Fig. 3, it becomesclear that primarily those parameters are selected that are associated with low-mass and low-eccentricity systems. Only a smallfraction of the whole volume of the parameter space is covered. where · denotes matrix multiplication, I f [ · ] interpolates vectors in frequency on a suitable grid, and A ( λ ) is anamplitude prefactor which is stored before the SVD takes place and an interpolant is computed over the parameterspace. V. ACCURACY AND SPEEDUP FOR SURROGATE MODEL PREDICTIONS
Once a ROM is built, any surrogate waveform can be evaluated as a sum of reduced basis elements with incrementalerrors within the parameter range covered in the particular model. The main criteria for a successful ROM are thatit facilitates data analysis applications that were infeasible with the fiducial waveform model and that it representswaveforms accurately. [48] This section is dedicated to appraise the overall performance of the ROM building discussedin Sec. IV. The first part of this section assesses the accuracy of surrogate model predictions in terms of the matchbetween the surrogate model and the fiducial model. In the second part we provide an overview of the computationalefficiency of the ROM with respect to computational complexity and cost relative to the cost of the fiducial model. - - - - - - - - FIG. 8: The linear trend in the change of surrogate error (39) as a function of the resolution of the frequency grid. Higherresolution of sampling times (i.e. lower resolution for sampling frequencies) result lesser uncertainty in estimating the amplitudesand phases. Surrogate error ∆˜ h = 1 . × − is marked with an orange point for a frequency-grid spacing ∆ f = 5 . × − M − which was obtained in eq. (24). The value of surrogate error corresponds to the mean relative error of the amplitude∆˜ h ( A ) ≈ × − shown in Fig. 9. A. Reconstruction errors
The overlap integral of two normalized waveforms, say, of a fiducial
CBwaves waveform h CB and its surrogatemodel prediction h S , is given by the mismatch (or unfaithfulness ) between the two waveforms and is defined as thenormalized inner product (38) maximized over time and phase shifts (cid:77) ≡ − max t ,φ (cid:104) h CB , h S (cid:105)(cid:107) h CB (cid:107)(cid:107) h S (cid:107) (37)with an inherited norm given by (cid:107) h (cid:107) ≡ (cid:104) h, h (cid:105) . A natural inner product between the two waveforms is given by thecomplex scalar product (cid:104) ˜ h CB , ˜ h S (cid:105) ≡ (cid:90) f high f low ˜ h CB ( f )˜ h ∗ S ( f ) S ˜ h ( f ) df (38)where the tilde denotes Fourier transformation given in eq. (15), ˜ h ∗ S ( f ) is the complex conjugate of ˜ h S ( f ), S ˜ h ( f ) isthe one-sided power spectral density (PSD) of the detector noise and f low , f high are suitable cutoff frequencies fordetector sensitivity. The low-frequency cutoff depends on the PSD and is at 10 Hz for advanced detectors design.The high-frequency cutoff is at 2 .
045 kHz, which is the ISCO frequency of the lowest total-mass configuration in ourfiducial model, discussed in Sec. III C. - - × - × - × - × - × - × - FIG. 9:
Top panel:
The amplitude and the phase part of the waveform associated with l = 1. There is visual agreementamong the fiducial CBwaves waveform and its surrogate prediction throughout the entire frequency range.
Bottom panel:
Therelative errors (41) with moving average of 50 points, defined by eq. (41), in the amplitude and the phase difference betweenthe fiducial waveform and its surrogate model prediction. The differences are smaller than the errors intrinsic to the surrogatemodel itself, as well as those of state-of-the-art numerical relativity simulations.
A discrete version of the normed difference between a fiducial waveform and its surrogate is what we can actuallymeasure: ∆˜ h ( λ ) = f s N − (cid:88) k =0 (cid:12)(cid:12)(cid:12) | ˜ h CB ( f k ; λ ) − ˜ h S ( f k ; λ ) (cid:12)(cid:12)(cid:12) (39)where f s is the sampling frequency discussed in Sec. III C. The square of the normed difference between two waveforms,refered to as the surrogate error , is directly related to their overlap (38). It is the dominant source of error in thesurrogate model that translates directly into errors in the fits of the parameters for building the surrogate. [33] Fig.8 shows the linear correlation of the surrogate error in eq. (39) with the time spacing ∆ f in the regularly spacedgrids (14–15). The surrogate model gradually converges to the fiducial one at finer time scales (i.e. larger sampling5frequencies). Other errors of interest are the pointwise ones (separately for the amplitude and phase). They areencoded in the l th surrogate model prediction (36) as˜ h ( A )S ( f k ; λ l ) ≡ I f [ B ( A ) · I [ M ( A ) ]( λ l )] , ˜ h ( φ )S ( f k ; λ l ) ≡ I f [ B (Φ) · I [ M (Φ) ]( λ l )] , (40)respectively. The relative errors in approximating the amplitude and phase of a fiducial waveform by its surrogatemodel prediction is then expressed by∆˜ h ( A ) ( f k ; λ l ) = (cid:12)(cid:12)(cid:12) − ˜ h ( A )S ( f k ; λ l ) / ˜ h ( A )CB ( f k ; λ l ) (cid:12)(cid:12)(cid:12) , ∆˜ h ( φ ) ( f k ; λ l ) = (cid:12)(cid:12)(cid:12) − ˜ h ( φ )S ( f k ; λ l ) / ˜ h ( φ )CB ( f k ; λ l ) (cid:12)(cid:12)(cid:12) (41)where the amplitude and phase parts of fiducial waveforms, ˜ h ( A )CB ( f k ; λ l ) and ˜ h ( φ )CB ( f k ; λ l ), respectively, are given byeq. (18) on N discrete frequency points f k .Fig. 9 shows a comparison between the surrogate and fiducial model, using the template assigned to l = 1. The toppanel shows that the fiducial and surrogate waveforms are visually indistinguishable. The bottom panel demonstratesthat both amplitude and phase pointwise errors (41) increase with frequency. Nevertheless, the errors are indeed assmall as predicted on Fig. 5. A moving average of 50 points was used to smooth out short-term fluctuations in theerror and highlight longer-term trends. B. Computational cost and speedup for surrogate model predictions
Apart from the requirements for accuracy or reliability, a ROM building is considered efficient if it generatescost-efficient surrogate models. The major advantage of using surrogate model predictions in lieu of actual waveformevaluations is their significantly reduced resource consumption. Now we discuss the computational cost, in termsof operation counts and run-time , of ROM building and present the desired speedup that can be achieved whenevaluating surrogate models.As described in Sec. IV B, the complete surrogate model (36) is assembled with the evaluation of r projectioncoefficients µ l ( f ) given in (33) and 2 r fitting functions { ˜ h ( A ) l ( λ ) } rl =1 and { ˜ h ( φ ) l ( λ ) } rl =1 given in (40). In order toconstruct a surrogate model for some parameter λ , one only needs to evaluate each of those 2 r fitting functions at λ ,recover the r complex values { ˜ h ( A ) l ( λ ) exp[ − i ˜ h ( φ ) l ( λ )] } rl =1 , and perform the summation over the index l . Each µ l ( f )is a complex-valued frequency series with N samples. Therefore, the total operation count to evaluate the surrogatemodel at each λ is (2 r − N plus the cost to evaluate the fitting functions. [33] The entire process of constructinga small, efficient ROM which is comprized of only r = 550 waveform templates sampled at N = 4000 grid pointsrequires the execution of approximately 4 . × operations (excluding the cost of evaluating the fitting functions).The notion of ‘speedup’, in our terminology, is the number that evaluates the relative performance of generatingthe same waveforms on the same processor by the execution of CBwaves code and of the surrogate model. Morespecifically, we test the acceleration of waveform generation by measures on the length of time required to performeach computational process. Let as note that the time which was denoted by t int and was referred to as ‘integration run-time’ in Sec. II A is actually the execution time during which the processor is actively working on our computations.It is referred to as CPU time (or run-time) and will be denoted by t CPU . In contrast, the actual elapsed real timeaccounts for the whole duration from when the computational process was started until the time it terminated. Thedifference between the two can arise from architecture and run-time dependent factors such as waiting for input/outputoperations (e.g. saving waveform templates). Consequently, the elapsed real time is greater than or equal to the CPUtime.Fig. 10 shows (on top) the computation time or CPU time for CBwaves waveforms (solid lines) against correspondingsurrogate waveforms (dashed lines) as a function of total mass of the binary system. The total mass M is measured inthe same 11 points as in Fig. 1 for three different initial eccentricities ( e = { . , . , } ) of equal-mass configurations,each associated with different colours. The computation time t CPU for surrogates is multiplied by a factor of 300 inorder to shift the curves close to their respective CBwaves counterparts and enable visual comparison. The bottompanel demonstrates that surrogates are several thousand times faster around 10 − M (cid:12) to evaluate as compared6 FIG. 10:
Top panel:
Computational time t CPU to generate fiducial waveforms by CBwaves code (dots; connected by solid lines)against the cost of evaluating corresponding surrogates by ROM (rectangles; connected by dashed lines). The computationaltime was measured for three different initial eccentricities of equal-mass configurations, each associated with different colours.
Bottom panel:
The speedup in evaluating the surrogate model is several thousand times faster around 10 − M (cid:12) thangenerating CBwaves waveforms. For high total mass the speedup falls off to several hundreds. The speedup is roughly twice asgreat for configurations having extremely high initial eccentricity at e = 0 .
98 (blue line) as for circular ones at e = 0 (greenline). to the cost of generating CBwaves waveforms. The speedup falls off to several hundreds as the total mass increases.Moreover, the speedup grows when the initial eccentricity e is increased in much the same way as with the massdisparity ¯ q (cf. Fig. 6 in [48]). The speedup is roughly twice as great for configurations having extremely high initialeccentricity ( e = 0 .
98) as for circular ones ( e = 0). The resemblance of the influence of e and ¯ q exercize on thespeedup can be attributed to their asymptotical nature as it had been pointed out earlier in Sec. II A. It is also evidentthat the speedup culminates when waveforms for configurations of very low total mass and very high eccentricity aregenerated. Such waveforms are prohibitively expensive to generate with CBwaves in contrast to surrogates that aregenerated at the same cost, regardless of the parameters of the configuration.Let us note that successive versions of SEOBNR ROMs have been developed and put to use within LAL (LSCAlgorithms Library) to shorten data analysis applications carried out since the first observation runs have begun.[52] It has been shown in [33, 34] that the cost of evaluating the surrogate model is linear in the number of samples N (cf. our Fig. 8 where the surrogate error depends on the sampling rate). Depending on the sampling rate, thespeedup is between 2 and almost 4 orders of magnitude. The speedup in evaluating surrogate models comparedto generating NR waveforms with the LAL analysis routines is crucial for searches and theoretical parameter esti-mations. SEOBNR (aligned-spin Effective-One-Body), IMRPhenomD (Inspiral-Merger-Ringdown PhenomenologicalModel ‘D’) and PhenSpinTaylorRD waveform approximants are among the best available GW models for genericspinning, compact binaries. In comparison with our results, the speedup achieved at the typical rate of 2 .
048 kHzused by aLIGO and aVirgo observatories is roughly 2300. [50]
VI. CONCLUDING REMARKS AND OUTLOOK
The primary goals of the present paper have been to propose a potential extension of the ROM techniques toalleviate the computational burden of constructing waveform templates for coalescing compact binaries with anyresidual orbital eccentricity and to validate the applicability of ROMs to this particular family of waveforms. ROMshave been applied to several waveform families (SEOBNR, IMRPhenomP and PhenSpinTaylorRD) in LAL routinesfor gravitational-wave data analysis. [34–41] The aforementioned waveform families provide efficient descriptions ofgravitational waves emitted during the late IMR (inspiral, merger, and ringdown) stages of compact binary systems,but only in the zero-eccentricity limit. The major motivation for extending the scope of application beyond the zero-eccentricity limit is based on the ground, referred to in Sec. I, that the great majority of compact objects formed indense stellar environments retain some non-negligible eccentricity when entering the frequency band of ground-basedGW detectors [3, 4], as well as the impact of eccentricity on the accuracy of parameter estimation for BNSs [15].Our approach to construct frequency-domain ROMs has been predominantly based on the method outlined inRefs. [47, 48] (see Sec. IV). Input waveforms comprised in the ROM are Fourier transformed and split into their7amplitude and phase parts (see Sec. III B). These functions are accuretely represented on a common, finely sampledand regularly spaced frequency grid defined in Sec. III C with only N = 4000 equidistant nodes, with a samplingfrequency recorded well above the required Nyquist frequency, at f s = 16 .
384 kHz. Fig. 8 demonstrates that,beside the degree of model order reduction, the accuracy of surrogate-waveform representation relies on the samplingfrequency. The upper and lower limits of frequency contained in the grid are determined from the ISCO frequencyfor the lowest total-mass configuration of interest (which is roughly 2 kHz in the present work) and the low-frequencycutoff of the detector noise spectrum (which is close to 10 Hz for aLIGO design). The ROM is designed to be capableof producing surrogates for GWs from coalescing compact binaries of total mass between 2 . M (cid:12) and 215 M (cid:12) , therebycovering the entire total-mass range of stellar-mass BBH/BNS systems of interest for ground-based GW detectors.The mass ratio is allowed to range between equal mass at q = 1 and relatively high mass-ratio at q ≈ .
01 while theinitial orbital eccentricity changes over a relatively wide range of values from e = 0 (circular orbits) up to e = 0 . M (cid:12) , which, of course, are excluded as inconceivable astrophysical sources. Despite the fact thatthe investigation has been restricted to a feasible 3-dimensional subset of the full 8-dimensional parameter space ofGW signals (see Fig. 3), the conclusions of Sec. IV, in agreement with that of Refs. [33–35, 47, 48], suggest that afull representation of the 8-parameter space might actually be achievable with a relatively compact reduced basis (cf.Ref. [33]). Template placement algorithms based on template-space metric (such as in Ref. [45, 46]) make admittedlyfar more effective coverage of the parameter space than the uniform spacings we used in this preliminary study. Asa matter of fact, Fig. 7 illustrates that the large majority of parameters of the selected templates constituting thereduced basis are concentrated along the axes of the parameter space.The reduced bases were built separatelly for the input amplitude and phase (see Fig. 6) by the decompositionof template matrices that comprise 550, 936, and 1800 input waveforms, respectively. The projection coefficientsfor corresponding input waveforms projected onto their reduced bases, were calculated as functions of the modelparameters ( M, q, e ) and were interpolated by tensor product cubic lines over the parameter space. Finally, theROM which preserves fundamental features of the original full-order model is assembled from its constituent parts.Fig. 5 demonstrates the underlying hierarchical nature of the generated template banks and indicates that thetruncation error in the approximated representation of surrogates decreases with the number of SVD componentsretained, characterized by a rate at which SVs decrease. Extremely little ( r (cid:46) r (cid:38) − r = 456). The first part of Sec. V assess the error of surrogate model predictions for waveforms thatwere originally not present in the original template bank, with special regard to the impact of frequency on thereconstruction error. To that end, reference waveforms were generated by CBwaves in all the intersection pointsright between the grid poinst of the original template bank (see the yellow in Fig. 3). Finally, the surrogates wereevaluated in the corresponding parameter-space points for comparison and the relative error was measured along allthe N = 4000 frequency points. The bottom panel of Fig. 9 attests that the relative error of the approximatedrepresentation is consistent with the error estimates derived from the singular values (∆˜ h ( A ) ≈ − , ∆˜ h ( φ ) ≈ − )over a large portion of the frequency range, but larger than expected at around the starting frequency (∆˜ h ( A ) ≈ − , ∆˜ h ( φ ) ≈ − ). The figure indicates that the relative error of the amplitude and phase increases with thefrequency. Our results provide clear examples of the construction and use of ROMs for eccentric inspiral waveforms.Our results also provide strong evidence that large increases in the speed of computation are obtained throughthe use of ROMs. Fig. 1 has exposed that the cost of computating input waveforms increases exponentially as thetotal mass decreases, but rises asymptotically at an even faster rate than the initial eccentricity or mass disparityincrease. In contrast to the cost of EOB waveform (full IMR) generation that rises steeply as the starting frequencyis decreased (see Ref. [48]), the cost of CBwaves waveform (inspiral-only) generation rises more gradually. The costof input waveform generation varies considerably in the region of parameter space ( M, q, e ) explored and depictedin Fig. 3, but Fig. 2 has revealed that only a surprisingly small fraction of waveforms of high-eccentricity andhigh-mass-disparity configurations are actually responsible for the prohibitively large time-consumption of integratinga large number of 3PN-accurate equation of motion over the investigated range of parameters. As discussed in thesecond part of Sec. V (based on Ref. [47]), the cost of generating surrogate waveforms (shown in the top panel ofFig. 10) comprises a constant cost of the spline interpolation at each frequency point and a cost of performing theinterpolations of coefficients over the parameter space. The speedup in evaluating the surrogate model, shown in thebottom panel of Fig. 10, is 2–3 orders of magnitude faster than generating corresponding CBwaves waveforms overall,reaching a factor of several thousand around 10–50 M (cid:12) .Finally, the method presented in this paper is limited to building surrogate models of inspiral-only PN inputwaveforms for the reason that eccentric binaries circularize in the last few cycles before the merger. Nevertheless,8composite waveforms that fully cover all the IMR stages can be constructed as prescribed in Ref. [47, 48] by matchingthe inspiral and NR waveforms of merger stages in either the time or frequency domain and then fitting this ‘hybrid’waveform to the ring-down part, described by damped exponentials. The gap between the initial part of the waveformand its final ring-down part, described by damped exponentials, is bridged by a phenomenological phase. The practicalimplementations of ‘hybrid’ waveforms that comprise eccentric inspirals of will be left for future work. We anticipatesubstantial speedup factors to come for predicting NR waveforms with a surrogate model compared to the expensivenumerical simulations for the same parameters. Developing an efficient template placement technique (such as in Ref.[45, 46]) for better coverage of the parameter space and an adaptive sampling technique in the frequency domain arecritical factors in the operational efficiency of ROMs and have been left for future work. All these ultimately leadingto computationally feasible and successful exploration of the full 8-dimensional parameter space of GW signals. Acknowledgments
This work was supported by the NKFIH 124366 and NKFIH 124508 grants. Partial support comes from NewComp-Star, COST Action Program MP1304. [1] B. Sathyaprakash and B. F. Schutz, Living Rev. Relat. , 2 (2009).[2] V. Kalogera et al., Astrophys. J. Lett. , L137 (2004).[3] P. C. Peters, Phys. Rev. , B1224 (1964).[4] F. Antonini and H. Perets, ApJ , 104 (2012).[5] S. L. Shapiro and S. A. Teukolsky, Astrophys. J. , L41 (1985).[6] V. Pierro and I. Pinto, Nuovo Cimento B , 1517 (1996).[7] K. Martel and E. Poisson, Phys. Rev. D , 124008 (1999).[8] E. A. Huerta, P. Kumar, S. T. McWilliams, R. O’Shaughnessy, and N. Yunes, Phys. Rev. D , 084016 (2014).[9] N. Yunes, K. G. Arun, E. Berti, and C. M. Will, Phys. Rev. D. , 084001 (2009).[10] C. K¨onigsd¨orffer and A. Gopakumar, Phys. Rev. D , 124012 (2006).[11] J. Samsing, M. MacLeod, and E. Ramirez-Ruiz, Astrophys. J. , 71 (2014).[12] R. M. O’Leary, B. Kocsis, and A. Loeb, Mon. Not. R. Astron. Soc. , 2127 (2009).[13] W. H. Lee, E. Ramirez-Ruiz, and G. van de Ven, Astrophys. J. , 953 (2010).[14] D. A. Brown and P. J. Zimmerman, Phys. Rev. D , 024007 (2010).[15] M. Favata, Phys. Rev. Lett. , 101101 (2014).[16] R. J. E. Smith, K. Cannon, C. Hanna, D. Keppel, and I. Mandel, Phys. Rev. D , 122002 (2013).[17] L. Blanchet, Living Rev. Relat. , 2 (2014).[18] P. Csizmadia, G. Debreczeni, I. R´acz, and M. Vas´uth, Class. Quant. Grav. , 245002 (2012).[19] P. Csizmadia, G. Debreczeni, I. R´acz, and M. Vas´uth, CBwaves simulation software (2011), URL http://grid.kfki.hu/project/virgo/cbwaves/ .[20] A. Buonanno, B. R. Iyer, Y. P. E. Ochsner, and B. S. Sathyaprakash, Phys. Rev. D , 084043 (2009).[21] B. Allen, W. G. Anderson, P. R. Brady, D. A. Brown, and J. D. Creighton, Phys. Rev. D , 122006 (2012).[22] B. J. Owen, Class. Quant. Grav. , 6749 (1996).[23] S. Babak, J. R. Gair, and E. K. Porter, Class. Quant. Grav. , 135004 (2009).[24] B. S. Sathyaprakash and S. V. Dhurandhar, Phys. Rev. D , 3819 (1991).[25] B. S. Sathyaprakash, Phys. Rev. D , 7111 (1994).[26] S. V. Dhurandhar and B. S. Sathyaprakash, Phys. Rev. D , 1707 (1994).[27] B. J. Owen and B. S. Sathyaprakash, Phys. Rev. D , 022002 (1999).[28] I. W. Harry, B. Allen, and B. S. Sathyaprakash, Phys. Rev. D , 104014 (2009).[29] P. Ca˜nizares, S. E. Field, J. R. Gair, and M. Tiglio, Phys. Rev. D , 124005 (2013).[30] J. Aasi and others (LIGO-Virgo Scientific Collaboration), Phys. Rev. D , 062001 (2013).[31] A. H. Mrou´e, M. A. Scheel, B. Szil´agyi, H. P. Pfeiffer, M. Boyle, et al., Phys. Rev. Lett. , 241104 (2013).[32] L. Pekowsky, R. O’Shaughnessy, J. Healy, and D. Shoemaker, Phys. Rev. D , 024040 (2013).[33] S. E. Field, C. R. Galley, J. S. Hesthaven, J. Kaye, and M. Tiglio, Phys. Rev. X , 031006 (2014).[34] S. E. Field, C. R. Galley, F. Herrmann, J. S. Hesthaven, E. Ochsner, and M. Tiglio, Phys. Rev. Lett. , 221102 (2011).[35] S. E. Field, C. R. Galley, and E. Ochsner, Phys. Rev. D , 084046 (2012).[36] K. Cannon, A. Chapman, C. Hanna, D. Keppel, A. C. Searle, and A. J. Weinstein, Phys. Rev. D , 044025 (2010).[37] K. Cannon, C. Hanna, and D. Keppel, Phys. Rev. D , 084003 (2011).[38] Y. Maday and O. Mula, NA , 221 (2013).[39] K. Cannon, C. Hanna, and D. Keppel, Phys. Rev. D , 081504 (2012).[40] K. Cannon, A. Chapman, C. Hanna, D. Keppel, A. C. Searle, and A. J. Weinstein, Phys. Rev. D , 044025 (2010). [41] S. Privitera, S. R. P. Mohapatra, P. Ajith, K. Cannon, N. Fotopoulos, M. A. Frei, C. Hanna, A. J. Weinstein, and J.Whelan, Phys. Rev. D , 024003 (2014).[42] P. Ca˜nizares, S. E. Field, J. R. Gair, V. Raymond, R. Smith, and M. Tiglio, Phys. Rev. Lett. , 071104 (2015).[43] L. E. Kidder, Phys. Rev. D , 821 (1995).[44] H. P. Pfeiffer, Class. Quant. Grav. , 124004 (2012).[45] C. Kalaghatgi, P. Ajith, and K. G. Arun, Phys. Rev. D , 124042 (2015).[46] T. Cokelaer, Phys. Rev. D , 102004 (2007).[47] M. P¨urrer, Class. Quant. Grav. , 195010 (2014).[48] M. P¨urrer, Phys. Rev. D , 064041 (2016).[49] E. Candes, P. Charlton, and H. Helgason, Appl. Comput. Harmon. Anal. , 14 (2008).[50] J. A. Aasi et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. D , 022002 (2013).[51] C. Eckart and G. Young, Psychometrika1