Search

In a previous work we introduced, in the context of gravitational wave science, an initial study on an automated domain-decomposition approach for reduced basis through hp-greedy refinement. The approach constructs local reduced bases of lower dimensionality than global ones, with the same or higher accuracy. These ``light'' local bases should imply both faster evaluations when predicting new waveforms and faster data analysis, in particular faster statistical inference (the forward and inverse problems, respectively). In this approach, however, we have previously found important dependence on several hyperparameters, which do not appear in global reduced basis. This naturally leads to the problem of hyperparameter optimization (HPO), which is the subject of this paper. We tackle the problem through a Bayesian optimization, and show its superiority when compared to grid or random searches. We find that for gravitational waves from the collision of two spinning but non-precessing black holes, for the same accuracy, local hp-greedy reduced bases with HPO have a lower dimensionality of up to $4 \times$ for the cases here studied, depending on the desired accuracy. This factor should directly translate in a parameter estimation speedup, for instance. Such acceleration might help in the near real-time requirements for electromagnetic counterparts of gravitational waves from compact binary coalescences. In addition, we find that the Bayesian approach used in this paper for HPO is two orders of magnitude faster than, for example, a grid search, with about a $100 \times$ acceleration. The code developed for this project is available as open source from public repositories., Comment: This paper is an invited contribution to the Special Issue "Recent Advances in Gravity: A Themed Issue in Honor of Prof. Jorge Pullin on his 60th Anniversary''

We introduce hp-greedy, a refinement approach for building gravitational wave surrogates as an extension of the standard reduced basis framework. Our proposal is data-driven, with a domain decomposition of the parameter space, local reduced basis, and a binary tree as the resulting structure, which are obtained in an automated way. When compared to the standard global reduced basis approach, the numerical simulations of our proposal show three salient features: i) representations of lower dimension with no loss of accuracy, ii) a significantly higher accuracy for a fixed maximum dimensionality of the basis, in some cases by orders of magnitude, and iii) results that depend on the reduced basis seed choice used by the refinement algorithm. We first illustrate the key parts of our approach with a toy model and then present a more realistic use case of gravitational waves emitted by the collision of two spinning, non-precessing black holes. We discuss performance aspects of hp-greedy, such as overfitting with respect to the depth of the tree structure, and other hyperparameter dependences. As two direct applications of the proposed hp-greedy refinement, we envision: i) a further acceleration of statistical inference, which might be complementary to focused reduced-order quadratures, and ii) the search of gravitational waves through clustering and nearest neighbors., Comment: 16 pages, code available from authors upon request

We analyze a prospect for predicting gravitational waveforms from compact binaries based on automated machine learning (AutoML) from around a hundred different possible regression models, without having to resort to tedious and manual case-by-case analyses and fine-tuning. The particular study of this article is within the context of the gravitational waves emitted by the collision of two spinless black holes in initial quasi-circular orbit. We find, for example, that approaches such as Gaussian process regression with radial bases as kernels do provide a sufficiently accurate solution, an approach which is generalizable to multiple dimensions with low computational evaluation cost. The results here presented suggest that AutoML might provide a framework for regression in the field of surrogates for gravitational waveforms. Our study is within the context of surrogates of numerical relativity simulations based on Reduced Basis and the Empirical Interpolation Method, where we find that for the particular case analyzed AutoML can produce surrogates which are essentially indistinguishable from the NR simulations themselves., Comment: 15 pages, 7 figures

We present an introduction to some of the state of the art in reduced order and surrogate modeling in gravitational wave (GW) science. Approaches that we cover include Principal Component Analysis, Proper Orthogonal Decomposition, the Reduced Basis approach, the Empirical Interpolation Method, Reduced Order Quadratures, and Compressed Likelihood evaluations. We divide the review into three parts: representation/compression of known data, predictive models, and data analysis. The targeted audience is that one of practitioners in GW science, a field in which building predictive models and data analysis tools that are both accurate and fast to evaluate, especially when dealing with large amounts of data and intensive computations, are necessary yet can be challenging. As such, practical presentations and, sometimes, heuristic approaches are here preferred over rigor when the latter is not available. This review aims to be self-contained, within reasonable page limits, with little previous knowledge (at the undergraduate level) requirements in mathematics, scientific computing, and other disciplines. Emphasis is placed on optimality, as well as the curse of dimensionality and approaches that might have the promise of beating it. We also review most of the state of the art of GW surrogates. Some numerical algorithms, conditioning details, scalability, parallelization and other practical points are discussed. The approaches presented are to large extent non-intrusive and data-driven and can therefore be applicable to other disciplines. We close with open challenges in high dimension surrogates, which are not unique to GW science., Comment: Invited article for Living Reviews in Relativity. 93 pages

The combination of the Reduced Basis and the Empirical Interpolation Method (EIM) approaches have produced outstanding results in many disciplines. In particular, in gravitational wave (GW) science these results range from building non-intrusive surrogate models for GWs to fast parameter estimation adding the use of Reduced Order Quadratures. These surrogates have the salient feature of being essentially indistinguishable from or very close to supercomputer simulations of the Einstein equations, but can be evaluated in the order of a milisecond per multipole mode on a standard laptop. In this article we clarify a common misperception of the EIM as originally introduced and used in practice in GW science. Namely, we prove that the EIM at each iteration chooses the interpolation nodes so as to make the related Vandermonde-type matrix as invertible as possible; not necessarily optimizing its conditioning or accuracy of the interpolant as is sometimes thought. In fact, we introduce two new variations of the EIM, nested as well, which do optimize with respect to conditioning and the Lebesgue constant, respectively, and compare them through numerical experiments with the original EIM using GWs. Our analyses and numerical results suggest a subtle relationship between solving for the original EIM, conditioning, and the Lebesgue constant, in consonance with active research in rigorous approximation theory and related fields., Comment: Accepted in Classical and Quantum Gravity. 16 pages, 7 figures. In this version the section on Fekete nodes was left for future research

We introduce a new approach for finding high accuracy, free and closed-form expressions for the gravitational waves emitted by binary black hole collisions from ab initio models. More precisely, our expressions are built from numerical surrogate models based on supercomputer simulations of the Einstein equations, which have been shown to be essentially indistinguishable from each other. Distinct aspects of our approach are that: i) representations of the gravitational waves can be explicitly written in a few lines, ii) these representations are free-form yet still fast to search for and validate and iii) there are no underlying physical approximations in the underlying model. The key strategy is combining techniques from Artificial Intelligence and Reduced Order Modeling for parameterized systems. Namely, symbolic regression through genetic programming combined with sparse representations in parameter space and the time domain using Reduced Basis and the Empirical Interpolation Method enabling fast free-form symbolic searches and large-scale a posteriori validations. As a proof of concept we present our results for the collision of two black holes, initially without spin, and with an initial separation corresponding to $25-31$ gravitational wave cycles before merger. The minimum overlap, compared to ground truth solutions, is $99\%$. That is, $1\%$ difference between our closed-form expressions and supercomputer simulations; this is considered for gravitational (GW) science more than the minimum required due to experimental numerical errors which otherwise dominate. This paper aims to contribute to the field of GWs in particular and Artificial Intelligence in general., Comment: 12 pages, 6 figures. Adding code for GWs and fix minor typos based on journal reports. Accepted in Nature Scientific Reports

On September 14, 2015, the newly upgraded Laser Interferometer Gravitational-wave Observatory (LIGO) recorded a loud gravitational-wave (GW) signal, emitted a billion light-years away by a coalescing binary of two stellar-mass black holes. The detection was announced in February 2016, in time for the hundredth anniversary of Einstein's prediction of GWs within the theory of general relativity (GR). The signal represents the first direct detection of GWs, the first observation of a black-hole binary, and the first test of GR in its strong-field, high-velocity, nonlinear regime. In the remainder of its first observing run, LIGO observed two more signals from black-hole binaries, one moderately loud, another at the boundary of statistical significance. The detections mark the end of a decades-long quest, and the beginning of GW astronomy: finally, we are able to probe the unseen, electromagnetically dark Universe by listening to it. In this article, we present a short historical overview of GW science: this young discipline combines GR, arguably the crowning achievement of classical physics, with record-setting, ultra-low-noise laser interferometry, and with some of the most powerful developments in the theory of differential geometry, partial differential equations, high-performance computation, numerical analysis, signal processing, statistical inference, and data science. Our emphasis is on the synergy between these disciplines, and how mathematics, broadly understood, has historically played, and continues to play, a crucial role in the development of GW science. We focus on black holes, which are very pure mathematical solutions of Einstein's gravitational-field equations that are nevertheless realized in Nature, and that provided the first observed signals., Comment: 41 pages, 5 figures. To appear in Bulletin of the American Mathematical Society

Simulating a binary black hole (BBH) coalescence by solving Einstein's equations is computationally expensive, requiring days to months of supercomputing time. Using reduced order modeling techniques, we construct an accurate surrogate model, which is evaluated in a millisecond to a second, for numerical relativity (NR) waveforms from non-spinning BBH coalescences with mass ratios in $[1, 10]$ and durations corresponding to about $15$ orbits before merger. We assess the model's uncertainty and show that our modeling strategy predicts NR waveforms {\em not} used for the surrogate's training with errors nearly as small as the numerical error of the NR code. Our model includes all spherical-harmonic ${}_{-2}Y_{\ell m}$ waveform modes resolved by the NR code up to $\ell=8.$ We compare our surrogate model to Effective One Body waveforms from $50$-$300 M_\odot$ for advanced LIGO detectors and find that the surrogate is always more faithful (by at least an order of magnitude in most cases)., Comment: Updated to published version, which includes a section comparing the surrogate and effective-one-body models. The surrogate is publicly available for download at http://www.black-holes.org/surrogates/ . 6 pages, 6 figures

Inferring the astrophysical parameters of coalescing compact binaries is a key science goal of the upcoming advanced LIGO-Virgo gravitational-wave detector network and, more generally, gravitational-wave astronomy. However, current parameter estimation approaches for such scenarios can lead to computationally intractable problems in practice. Therefore there is a pressing need for new, fast and accurate Bayesian inference techniques. In this letter we demonstrate that a reduced order modeling approach enables rapid parameter estimation studies. By implementing a reduced order quadrature scheme within the LIGO Algorithm Library, we show that Bayesian inference on the 9-dimensional parameter space of non-spinning binary neutron star inspirals can be sped up by a factor of 30 for the early advanced detectors' configurations. This speed-up will increase to about $150$ as the detectors improve their low-frequency limit to 10Hz, reducing to hours analyses which would otherwise take months to complete. Although these results focus on gravitational detectors, the techniques are broadly applicable to any experiment where fast Bayesian analysis is desirable., Comment: To appear in Physical Review Letters

Many relevant applications in gravitational wave physics share a significant common problem: the seven-dimensional parameter space of gravitational waveforms from precessing compact binary inspirals and coalescences is large enough to prohibit covering the space of waveforms with sufficient density. We find that by using the reduced basis method together with a parametrization of waveforms based on their phase and precession, we can construct ultra-compact yet high-accuracy representations of this large space. As a demonstration, we show that less than $100$ judiciously chosen precessing inspiral waveforms are needed for $200$ cycles, mass ratios from $1$ to $10$ and spin magnitudes $\le 0.9$. In fact, using only the first $10$ reduced basis waveforms yields a maximum mismatch of $0.016$ over the whole range of considered parameters. We test whether the parameters selected from the inspiral regime result in an accurate reduced basis when including merger and ringdown; we find that this is indeed the case in the context of a non-precessing effective-one-body model. This evidence suggests that as few as $\sim 100$ numerical simulations of binary black hole coalescences may accurately represent the seven-dimensional parameter space of precession waveforms for the considered ranges., Comment: 5 pages, 3 figures. The parameters selected for the basis of precessing waveforms can be found in the source files

[Abridged] We propose a solution to the problem of quickly and accurately predicting gravitational waveforms within any given physical model. The method is relevant for both real-time applications and in more traditional scenarios where the generation of waveforms using standard methods can be prohibitively expensive. Our approach is based on three offline steps resulting in an accurate reduced-order model that can be used as a surrogate for the true/fiducial waveform family. First, a set of m parameter values is determined using a greedy algorithm from which a reduced basis representation is constructed. Second, these m parameters induce the selection of m time values for interpolating a waveform time series using an empirical interpolant. Third, a fit in the parameter dimension is performed for the waveform's value at each of these m times. The cost of predicting L waveform time samples for a generic parameter choice is of order m L + m c_f online operations where c_f denotes the fitting function operation count and, typically, m << L. We generate accurate surrogate models for Effective One Body (EOB) waveforms of non-spinning binary black hole coalescences with durations as long as 10^5 M, mass ratios from 1 to 10, and for multiple harmonic modes. We find that these surrogates are three orders of magnitude faster to evaluate as compared to the cost of generating EOB waveforms in standard ways. Surrogate model building for other waveform models follow the same steps and have the same low online scaling cost. For expensive numerical simulations of binary black hole coalescences we thus anticipate large speedups in generating new waveforms with a surrogate. As waveform generation is one of the dominant costs in parameter estimation algorithms and parameter space exploration, surrogate models offer a new and practical way to dramatically accelerate such studies without impacting accuracy., Comment: 20 pages, 17 figures, uses revtex 4.1. Version 2 includes new numerical experiments for longer waveform durations, larger regions of parameter space and multi-mode models

One of the main bottlenecks in gravitational wave (GW) astronomy is the high cost of performing parameter estimation and GW searches on the fly. We propose a novel technique based on Reduced Order Quadratures (ROQs), an application and data-specific quadrature rule, to perform fast and accurate likelihood evaluations. These are the dominant cost in Markov chain Monte Carlo (MCMC) algorithms, which are widely employed in parameter estimation studies, and so ROQs offer a new way to accelerate GW parameter estimation. We illustrate our approach using a four dimensional GW burst model embedded in noise. We build an ROQ for this model, and perform four dimensional MCMC searches with both the standard and ROQs quadrature rules, showing that, for this model, the ROQ approach is around 25 times faster than the standard approach with essentially no loss of accuracy. The speed-up from using ROQs is expected to increase for more complex GW signal models and therefore has significant potential to accelerate parameter estimation of GW sources such as compact binary coalescences., Comment: 18 pages, 12 figures (2 updated) and 2 tables. Minor edits based on referee report

We present an algorithm to generate application-specific, global reduced order quadratures (ROQ) for multiple fast evaluations of weighted inner products between parameterized functions. If a reduced basis (RB) or any other projection-based model reduction technique is applied, the dimensionality of integrands is reduced dramatically; however, the cost of approximating the integrands by projection still scales as the size of the original problem. In contrast, using discrete empirical interpolation (DEIM) points as ROQ nodes leads to a computational cost which depends linearly on the dimension of the reduced space. Generation of a reduced basis via a greedy procedure requires a training set, which for products of functions can be very large. Since this direct approach can be impractical in many applications, we propose instead a two-step greedy targeted towards approximation of such products. We present numerical experiments demonstrating the accuracy and the efficiency of the two-step approach. The presented ROQ are expected to display very fast convergence whenever there is regularity with respect to parameter variation. We find that for the particular application here considered, one driven by gravitational wave physics, the two-step approach speeds up the offline computations to build the ROQ by more than two orders of magnitude. Furthermore, the resulting ROQ rule is found to converge exponentially with the number of nodes, and a factor of ~50 savings, without loss of accuracy, is observed in evaluations of inner products when ROQ are used as a downsampling strategy for equidistant samples using the trapezoidal rule. While the primary focus of this paper is on quadrature rules for inner products of parameterized functions, our method can be easily adapted to integrations of single parameterized functions, and some examples of this type are considered., Comment: 27 pages, 9 figures, uses svjour3

Using the Reduced Basis approach, we efficiently compress and accurately represent the space of waveforms for non-precessing binary black hole inspirals, which constitutes a four dimensional parameter space (two masses, two spin magnitudes). Compared to the non-spinning case, we find that only a {\it marginal} increase in the (already relatively small) number of reduced basis elements is required to represent any non-precessing waveform to nearly numerical round-off precision. Most parameters selected by the algorithm are near the boundary of the parameter space, leaving the bulk of its volume sparse. Our results suggest that the full eight dimensional space (two masses, two spin magnitudes, four spin orientation angles on the unit sphere) may be highly compressible and represented with very high accuracy by a remarkably small number of waveforms, thus providing some hope that the number of numerical relativity simulations of binary black hole coalescences needed to represent the entire space of configurations is not intractable. Finally, we find that the {\it distribution} of selected parameters is robust to different choices of seed values starting the algorithm, a property which should be useful for indicating parameters for numerical relativity simulations of binary black holes. In particular, we find that the mass ratios $m_1/m_2$ of non-spinning binaries selected by the algorithm are mostly in the interval $[1,3]$ and that the median of the distribution follows a power-law behavior $\sim (m_1/m_2)^{-5.25}$.

Many evolution problems in physics are described by partial differential equations on an infinite domain; therefore, one is interested in the solutions to such problems for a given initial dataset. A prominent example is the binary black hole problem within Einstein's theory of gravitation, in which one computes the gravitational radiation emitted from the inspiral of the two black holes, merger and ringdown. Powerful mathematical tools can be used to establish qualitative statements about the solutions, such as their existence, uniqueness, continuous dependence on the initial data, or their asymptotic behavior over large time scales. However, one is often interested in computing the solution itself, and unless the partial differential equation is very simple, or the initial data possesses a high degree of symmetry, this computation requires approximation by numerical discretization. When solving such discrete problems on a machine, one is faced with a finite limit to computational resources, which leads to the replacement of the infinite continuum domain with a finite computer grid. This, in turn, leads to a discrete initial-boundary value problem. The hope is to recover, with high accuracy, the exact solution in the limit where the grid spacing converges to zero with the boundary being pushed to infinity. The goal of this article is to review some of the theory necessary to understand the continuum and discrete initial-boundary value problems arising from hyperbolic partial differential equations and to discuss its applications to numerical relativity; in particular, we present well-posed initial and initial-boundary value formulations of Einstein's equations, and we discuss multi-domain high-order finite difference and spectral methods to solve them., Comment: To appear in Living Reviews in Relativity, 191 pages

We present a new fully first order strongly hyperbolic representation of the BSSN formulation of Einstein's equations with optional constraint damping terms. We describe the characteristic fields of the system, discuss its hyperbolicity properties, and present two numerical implementations and simulations: one using finite differences, adaptive mesh refinement and in particular binary black holes, and another one using the discontinuous Galerkin method in spherical symmetry. The results of this paper constitute a first step in an effort to combine the robustness of BSSN evolutions with very high accuracy numerical techniques, such as spectral collocation multi-domain or discontinuous Galerkin methods., Comment: To appear in Physical Review D

We construct compact and high accuracy Reduced Basis (RB) representations of single and multiple quasinormal modes (QNMs). The RB method determines a hierarchical and relatively small set of the most relevant waveforms. We find that the exponential convergence of the method allows for a dramatic compression of template banks used for ringdown searches. Compressing a catalog with a minimal match $\MMm=0.99$, we find that the selected RB waveforms are able to represent {\em any} QNM, including those not in the original bank, with extremely high accuracy, typically less than $10^{-13}$. We then extend our studies to two-mode QNMs. Inclusion of a second mode is expected to help with detection, and might make it possible to infer details of the progenitor of the final black hole. We find that the number of RB waveforms needed to represent any two-mode ringdown waveform with the above high accuracy is {\em smaller} than the number of metric-based, one-mode templates with $\MMm=0.99$. For unconstrained two-modes, which would allow for consistency tests of General Relativity, our high accuracy RB has around $10^4$ {\em fewer} waveforms than the number of metric-based templates for $\MMm=0.99$. The number of RB elements grows only linearly with the number of multipole modes versus exponentially with the standard approach, resulting in very compact representations even for many multiple modes. The results of this paper open the possibility of searches of multi-mode ringdown gravitational waves., Comment: Edits to match the final version to appear in Classical and Quantum Gravity

We introduce a reduced basis approach as a new paradigm for modeling, representing and searching for gravitational waves. We construct waveform catalogs for non-spinning compact binary coalescences, and we find that for accuracies of 99% and 99.999% the method generates a factor of about $10-10^5$ fewer templates than standard placement methods. The continuum of gravitational waves can be represented by a finite and comparatively compact basis. The method is robust under variations in the noise of detectors, implying that only a single catalog needs to be generated., Comment: Minor changes in some of the phrasing to match the version as published in PRL

Within linearized perturbation theory, black holes decay to their final stationary state through the well-known spectrum of quasinormal modes. Here we numerically study whether nonlinearities change this picture. For that purpose we study the ringdown frequencies of gauge-invariant second-order gravitational perturbations induced by self-coupling of linearized perturbations of Schwarzschild black holes. We do so through high-accuracy simulations in the time domain of first and second-order Regge-Wheeler-Zerilli type equations, for a variety of initial data sets. We consider first-order even-parity $(\ell=2,m=\pm 2)$ perturbations and odd-parity $(\ell=2,m=0)$ ones, and all the multipoles that they generate through self-coupling. For all of them and all the initial data sets considered we find that ---in contrast to previous predictions in the literature--- the numerical decay frequencies of second-order perturbations are the same ones of linearized theory, and we explain the observed behavior. This would indicate, in particular, that when modeling or searching for ringdown gravitational waves, appropriately including the standard quasinormal modes already takes into account nonlinear effects.

We perform a statistical analysis of the binary black hole problem in the post-Newtonian approximation by systematically sampling and evolving the parameter space of initial configurations for quasi-circular inspirals. Through a principal component analysis of spin and orbital angular momentum variables we systematically look for uncorrelated quantities and find three of them which are highly conserved in a statistical sense, both as functions of time and with respect to variations in initial spin orientations. We also look for and find the variables that account for the largest variations in the problem. We present binary black hole simulations of the full Einstein equations analyzing to what extent these results might carry over to the full theory in the inspiral and merger regimes. Among other applications these results should be useful both in semi-analytical and numerical building of templates of gravitational waves for gravitational wave detectors.

We report on early results of a numerical and statistical study of binary black hole inspirals. The two black holes are evolved using post-Newtonian approximations starting with initially randomly distributed spin vectors. We characterize certain aspects of the distribution shortly before merger. In particular we note the uniform distribution of black hole spin vector dot products shortly before merger and a high correlation between the initial and final black hole spin vector dot products in the equal-mass, maximally spinning case. These simulations were performed on Graphics Processing Units, and we demonstrate a speed-up of a factor 50 over a more conventional CPU implementation., Comment: Added one reference

We present two methods to include the asymptotic domain of a background spacetime in null directions for numerical solutions of evolution equations so that both the radiation extraction problem and the outer boundary problem are solved. The first method is based on the geometric conformal approach, the second is a coordinate based approach. We apply these methods to the case of a massless scalar wave equation on a Kerr spacetime. Our methods are designed to allow existing codes to reach the radiative zone by including future null infinity in the computational domain with minor modifications. We demonstrate the flexibility of the methods by considering both Boyer-Lindquist and ingoing Kerr coordinates near the black hole. We also confirm numerically predictions concerning tail decay rates for scalar fields at null infinity in Kerr spacetime due to Hod for the first time., Comment: 12 pages, 9 figures; matches the published version

We present numerical simulations of orbiting black holes for around twelve cycles, using a high-order multipatch approach. Unlike some other approaches, the computational speed scales almost perfectly for thousands of processors. Multipatch methods are an alternative to AMR (adaptive mesh refinement), with benefits of simplicity and better scaling for improving the resolution in the wave zone. The results presented here pave the way for multipatch evolutions of black hole-neutron star and neutron star-neutron star binaries, where high resolution grids are needed to resolve details of the matter flow.

Using recently developed efficient symbolic manipulations tools, we present a general gauge-invariant formalism to study arbitrary radiative $(l\geq 2)$ second-order perturbations of a Schwarzschild black hole. In particular, we construct the second order Zerilli and Regge-Wheeler equations under the presence of any two first-order modes, reconstruct the perturbed metric in terms of the master scalars, and compute the radiated energy at null infinity. The results of this paper enable systematic studies of generic second order perturbations of the Schwarzschild spacetime. In particular, studies of mode-mode coupling and non-linear effects in gravitational radiation, the second-order stability of the Schwarzschild spacetime, or the geometry of the black hole horizon., Comment: 14 pages

We compute the contribution to the Lagrangian from the leading order (2.5 post-Newtonian) radiation reaction and the quadrupolar gravitational waves emitted from a binary system using the effective field theory (EFT) approach of Goldberger and Rothstein. We use an initial value formulation of the underlying (quantum) framework to implement retarded boundary conditions and describe these real-time dissipative processes. We also demonstrate why the usual scattering formalism of quantum field theory inadequately accounts for these. The methods discussed here should be useful for deriving real-time quantities (including radiation reaction forces and gravitational wave emission) and hereditary terms in the post-Newtonian approximation (including memory, tail and other causal, history-dependent integrals) within the EFT approach. We also provide a consistent formulation of the radiation sector in the equivalent effective field theory approach of Kol and Smolkin., Comment: 23 pages, 8 figures

We provide a detailed analysis of several aspects of the turduckening technique for evolving black holes. At the analytical level we study the constraint propagation for a general family of BSSN-type formulation of Einstein's field equations and identify under what conditions the turducken procedure is rigorously justified and under what conditions constraint violations will propagate to the outside of the black holes. We present high-resolution spherically symmetric studies which verify our analytical predictions. Then we present three dimensional simulations of single distorted black holes using different variations of the turduckening method and also the puncture method. We study the effect that these different methods have on the coordinate conditions, constraint violations, and extracted gravitational waves. We find that the waves agree up to small but non-vanishing differences, caused by escaping superluminal gauge modes. These differences become smaller with increasing detector location., Comment: Minor changes to match the final version to appear in PRD

In order to generate initial data for nonlinear relativistic simulations, one needs to solve the Einstein constraints, which can be cast into a coupled set of nonlinear elliptic equations. Here we present an approach for solving these equations on three-dimensional multi-block domains using finite element methods. We illustrate our approach on a simple example of Brill wave initial data, with the constraints reducing to a single linear elliptic equation for the conformal factor $\psi$. We use quadratic Lagrange elements on semi-structured simplicial meshes, obtained by triangulation of multi-block grids. In the case of uniform refinement the scheme is superconvergent at most mesh vertices, due to local symmetry of the finite element basis with respect to local spatial inversions. We show that in the superconvergent case subsequent unstructured mesh refinements do not improve the quality of our initial data. As proof of concept that this approach is feasible for generating multi-block initial data in three dimensions, after constructing the initial data we evolve them in time using a high order finite-differencing multi-block approach and extract the gravitational waves from the numerical solution., Comment: Changes made to match the version to appear in Classical and Quantum Gravity

We investigate the late time behavior of a scalar field on a fixed Kerr background using a 2+1 dimensional pseudospectral evolution code. We compare evolutions of pure axisymmetric multipoles in both Kerr-Schild and Boyer-Lindquist coordinates. We find that the late-time power-law decay rate depends upon the slicing of the background, confirming previous theoretical predictions for those decay rates. The accuracy of the numerical evolutions is sufficient to decide unambiguously between competing claims in the literature.

Many systems of interest in general relativistic astrophysics, including neutron stars, accreting compact objects in X-ray binaries and active galactic nuclei, core collapse, and collapsars, are assumed to be approximately spherically symmetric or axisymmetric. In Newtonian or fixed-background relativistic approximations it is common practice to use spherical polar coordinates for computational grids; however, these coordinates have singularities and are difficult to use in fully relativistic models. We present, in this series of papers, a numerical technique which is able to use effectively spherical grids by employing multiple patches. We provide detailed instructions on how to implement such a scheme, and present a number of code tests for the fixed background case, including an accretion torus around a black hole., Comment: 26 pages, 20 figures. A high-resolution version is available at http://www.cct.lsu.edu/~bzink/papers/multipatch_1.pdf

to turducken (turduckens, turduckening, turduckened, turduckened) [math.]: To stuff a black hole. We analyze and apply an alternative to black hole excision based on smoothing the interior of black holes with arbitrary - possibly constraint violating - initial data, and solving the vacuum Einstein evolution equations everywhere. By deriving the constraint propagation system for our hyperbolic formulation of the BSSN evolution system we rigorously prove that the constraints propagate causally and so any constraint violations introduced inside the black holes cannot affect the exterior spacetime. (This does not follow from the causal structure of the spacetime as is often assumed.) We present numerical evolutions of Cook-Pfeiffer binary black hole initial configurations showing that these techniques appear to work robustly for generic data. We also present numerical evidence from spherically symmetric evolutions that for the gauge conditions used the same stationary end-state is approached irrespective of the choice of initial data and smoothing procedure., Comment: 5 pages, 4 figures. Uses improved smoothing procedure for initial data. A free turduckening implementation is available from the authors by anonymous CVS from :pserver:cvs_anon@cvs.aei.mpg.de:/numrelcvs as AEIThorns/NoExcision; use the password "anon"

We present a method for extracting gravitational waves from numerical spacetimes which generalizes and refines one of the standard methods based on the Regge--Wheeler--Zerilli perturbation formalism. [abridged] We then present fully nonlinear three-dimensional numerical evolutions of a distorted Schwarzschild black hole in Kerr--Schild coordinates with an odd parity perturbation and analyze the improvement we gain from our generalized wave extraction, comparing our new method to the standard one. [abridged] We find that, even with observers as far out as $R=80 M$--which is larger than what is commonly used in state-of-the-art simulations--the assumption in the standard method that the background is close to having Schwarzschild-like coordinates increases the error in the extracted waveforms considerably. Even for our coarsest resolutions, our new method decreases the error by between one and two orders of magnitudes. Furthermore, we explicitly see that the errors in the extracted waveforms obtained by the standard method do not converge to zero with increasing resolution. [abridged] In a general scenario, for example a collision of compact objects, there is no precise definition of gravitational radiation at a finite distance, and gravitational wave extraction methods at such distances are thus inherently approximate. The results of this paper bring up the possibility that different choices in the wave extraction procedure at a fixed and finite distance may result in relative differences in the waveforms which are actually larger than the numerical errors in the solution., Comment: 19 pages

We present numerical results from three-dimensional evolutions of scalar perturbations of Kerr black holes. Our simulations make use of a high-order accurate multi-block code which naturally allows for fixed adaptivity and smooth inner (excision) and outer boundaries. We focus on the quasinormal ringing phase, presenting a systematic method for extraction of the quasinormal mode frequencies and amplitudes and comparing our results against perturbation theory. The amplitude of each mode depends exponentially on the starting time of the quasinormal regime, which is not defined unambiguously. We show that this time-shift problem can be circumvented by looking at appropriately chosen relative mode amplitudes. From our simulations we extract the quasinormal frequencies and the relative and absolute amplitudes of corotating and counterrotating modes (including overtones in the corotating case). We study the dependence of these amplitudes on the shape of the initial perturbation, the angular dependence of the mode and the black hole spin, comparing against results from perturbation theory in the so-called asymptotic approximation. We also compare the quasinormal frequencies from our numerical simulations with predictions from perturbation theory, finding excellent agreement. Finally we study under what conditions the relative amplitude between given pairs of modes gets maximally excited and present a quantitative analysis of rotational mode-mode coupling. The main conclusions and techniques of our analysis are quite general and, as such, should be of interest in the study of ringdown gravitational waves produced by astrophysical gravitational wave sources., Comment: 21 pages, 9 figures

We describe a generic infrastructure for time evolution simulations in numerical relativity using multiple grid patches. After a motivation of this approach, we discuss the relative advantages of global and patch-local tensor bases. We describe both our multi-patch infrastructure and our time evolution scheme, and comment on adaptive time integrators and parallelisation. We also describe various patch system topologies that provide spherical outer and/or multiple inner boundaries. We employ penalty inter-patch boundary conditions, and we demonstrate the stability and accuracy of our three-dimensional implementation. We solve both a scalar wave equation on a stationary rotating black hole background and the full Einstein equations. For the scalar wave equation, we compare the effects of global and patch-local tensor bases, different finite differencing operators, and the effect of artificial dissipation onto stability and accuracy. We show that multi-patch systems can directly compete with the so-called fixed mesh refinement approach; however, one can also combine both. For the Einstein equations, we show that using multiple grid patches with penalty boundary conditions leads to a robustly stable system. We also show long-term stable and accurate evolutions of a one-dimensional non-linear gauge wave. Finally, we evolve weak gravitational waves in three dimensions and extract accurate waveforms, taking advantage of the spherical shape of our grid lines., Comment: 18 pages. Some clarifications added, figure layout improved

We construct new, efficient, and accurate high-order finite differencing operators which satisfy summation by parts. Since these operators are not uniquely defined, we consider several optimization criteria: minimizing the bandwidth, the truncation error on the boundary points, the spectral radius, or a combination of these. We examine in detail a set of operators that are up to tenth order accurate in the interior, and we surprisingly find that a combination of these optimizations can improve the operators' spectral radius and accuracy by orders of magnitude in certain cases. We also construct high-order dissipation operators that are compatible with these new finite difference operators and which are semi-definite with respect to the appropriate summation by parts scalar product. We test the stability and accuracy of these new difference and dissipation operators by evolving a three-dimensional scalar wave equation on a spherical domain consisting of seven blocks, each discretized with a structured grid, and connected through penalty boundary conditions., Comment: 16 pages, 9 figures. The files with the coefficients for the derivative and dissipation operators can be accessed by downloading the source code for the document. The files are located in the "coeffs" subdirectory

During the last few years progress has been made on several fronts making it possible to revisit Cauchy-perturbative matching (CPM) in numerical relativity in a more robust and accurate way. This paper is the first in a series where we plan to analyze CPM in the light of these new results. Here we start by testing high-order summation-by-parts operators, penalty boundaries and contraint-preserving boundary conditions applied to CPM in a setting that is simple enough to study all the ingredients in great detail: Einstein's equations in spherical symmetry, describing a black hole coupled to a massless scalar field. We show that with the techniques described above, the errors introduced by Cauchy-perturbative matching are very small, and that very long term and accurate CPM evolutions can be achieved. Our tests include the accretion and ring-down phase of a Schwarzschild black hole with CPM, where we find that the discrete evolution introduces, with a low spatial resolution of \Delta r = M/10, an error of 0.3% after an evolution time of 1,000,000 M. For a black hole of solar mass, this corresponds to approximately 5 s, and is therefore at the lower end of timescales discussed e.g. in the collapsar model of gamma-ray burst engines. (abridged), Comment: 14 pages, 20 figures

The need to smoothly cover a computational domain of interest generically requires the adoption of several grids. To solve the problem of interest under this grid-structure one must ensure the suitable transfer of information among the different grids involved. In this work we discuss a technique that allows one to construct finite difference schemes of arbitrary high order which are guaranteed to satisfy linear numerical and strict stability. The technique relies on the use of difference operators satisfying summation by parts and {\it penalty techniques} to transfer information between the grids. This allows the derivation of semidiscrete energy estimates for problems admitting such estimates at the continuum. We analyze several aspects of this technique when used in conjuction with high order schemes and illustrate its use in one, two and three dimensional numerical relativity model problems with non-trivial topologies, including truly spherical black hole excision., Comment: 35 pages

Outer boundary conditions for strongly and symmetric hyperbolic formulations of 3D Einstein's field equations with a live gauge condition are discussed. The boundary conditions have the property that they ensure constraint propagation and control in a sense made precise in this article the physical degrees of freedom at the boundary. We use Fourier-Laplace transformation techniques to find necessary conditions for the well posedness of the resulting initial-boundary value problem and integrate the resulting three-dimensional nonlinear equations using a finite-differencing code. We obtain a set of constraint-preserving boundary conditions which pass a robust numerical stability test. We explicitly compare these new boundary conditions to standard, maximally dissipative ones through Brill wave evolutions. Our numerical results explicitly show that in the latter case the constraint variables, describing the violation of the constraints, do not converge to zero when resolution is increased while for the new boundary conditions, the constraint variables do decrease as resolution is increased. As an application, we inject pulses of ``gravitational radiation'' through the boundaries of an initially flat spacetime domain, with enough amplitude to generate strong fields and induce large curvature scalars, showing that our boundary conditions are robust enough to handle nonlinear dynamics. We expect our boundary conditions to be useful for improving the accuracy and stability of current binary black hole and binary neutron star simulations, for a successful implementation of characteristic or perturbative matching techniques, and other applications. We also discuss limitations of our approach and possible future directions., Comment: 31 pages, updated version, to appear in Journal of Hyperbolic Differential Equations

Combining deeper insight of Einstein's equations with sophisticated numerical techniques promises the ability to construct accurate numerical implementations of these equations. We illustrate this in two examples, the numerical evolution of ``bubble'' and single black hole spacetimes. The former is chosen to demonstrate how accurate numerical solutions can answer open questions and even reveal unexpected phenomena. The latter illustrates some of the difficulties encountered in three-dimensional black hole simulations, and presents some possible remedies., Comment: Contribution to appear in the Proceedings of the March-2004 Heraeus Seminar in Bad Honnef (Eds D. Giulini, J. Frauendiener and V. Perlick)

The energy method can be used to identify well-posed initial boundary value problems for quasi-linear, symmetric hyperbolic partial differential equations with maximally dissipative boundary conditions. A similar analysis of the discrete system can be used to construct stable finite difference equations for these problems at the linear level. In this paper we apply these techniques to some test problems commonly used in numerical relativity and observe that while we obtain convergent schemes, fast growing modes, or ``artificial instabilities,'' contaminate the solution. We find that these growing modes can partially arise from the lack of a Leibnitz rule for discrete derivatives and discuss ways to limit this spurious growth., Comment: 18 pages, 22 figures

We present three-dimensional simulations of Einstein equations implementing a symmetric hyperbolic system of equations with dynamical lapse. The numerical implementation makes use of techniques that guarantee linear numerical stability for the associated initial-boundary value problem. The code is first tested with a gauge wave solution, where rather larger amplitudes and for significantly longer times are obtained with respect to other state of the art implementations. Additionally, by minimizing a suitably defined energy for the constraints in terms of free constraint-functions in the formulation one can dynamically single out preferred values of these functions for the problem at hand. We apply the technique to fully three-dimensional simulations of a stationary black hole spacetime with excision of the singularity, considerably extending the lifetime of the simulations., Comment: 21 pages. To appear in PRD

We discuss finite difference techniques for hyperbolic equations in non-trivial domains, as those that arise when simulating black hole spacetimes. In particular, we construct dissipative and difference operators that satisfy the {\it summation by parts} property in domains with excised multiple cubic regions. This property can be used to derive semi-discrete energy estimates for the associated initial-boundary value problem which in turn can be used to prove numerical stability., Comment: Some minor changes in presentation, to appear in CQG. 20 pages

I present a new, simple method to dynamically control the growth of the discretized constraints during a free evolution of Einstein's equations. During an evolution, any given family of formulations is adjusted off the constraints surface in a way such that, for any chosen numerical method and arbitrary but fixed resolution, the constraints growth can be minimized with respect to the freedom allowed by the formulation. In particular, provided there is enough freedom, the discretized constraints can be maintained close to its initial truncation value for all times, or decay from it. No a priori knowledge of the solution is needed, and the method can be applied to any formulation of Einstein's equations without affecting hyperbolicity. This method is independent of the numerical algorithm and accounts for constraint violating modes introduced both by continuum instabilities of the formulation and by the numerical method.

We use rigorous techniques from numerical analysis of hyperbolic equations in bounded domains to construct stable finite-difference schemes for Numerical Relativity, in particular for their use in black hole excision. As an application, we present 3D simulations of a scalar field propagating in a Schwarzschild black hole background., Comment: 4 pages, RevTex, 2 figures

In the Cauchy problem of general relativity one considers initial data that satisfies certain constraints. The evolution equations guarantee that the evolved variables will satisfy the constraints at later instants of time. This is only true within the domain of dependence of the initial data. If one wishes to consider situations where the evolutions are studied for longer intervals than the size of the domain of dependence, as is usually the case in three dimensional numerical relativity, one needs to give boundary data. The boundary data should be specified in such a way that the constraints are satisfied everywhere, at all times. In this paper we address this problem for the case of general relativity linearized around Minkowski space using the generalized Einstein-Christoffel symmetric hyperbolic system of evolution equations. We study the evolution equations for the constraints, specify boundary conditions for them that make them well posed and further choose these boundary conditions in such a way that the evolution equations for the metric variables are also well posed. We also consider the case of a manifold with a non-smooth boundary, as is the usual case of the cubic boxes commonly used in numerical relativity. The techniques discussed should be applicable to more general cases, as linearizations around more complicated backgrounds, and may be used to establish well posedness in the full non-linear case., Comment: 15 pages, RevTeX, no figures

It is often the case in numerical relativity that schemes that are known to be convergent for well posed systems are used in evolutions of weakly hyperbolic (WH) formulations of Einstein's equations. Here we explicitly show that with several of the discretizations that have been used through out the years, this procedure leads to non-convergent schemes. That is, arbitrarily small initial errors are amplified without bound when resolution is increased, independently of the amount of numerical dissipation introduced. The lack of convergence introduced by this instability can be particularly subtle, in the sense that it can be missed by several convergence tests, especially in 3+1 dimensional codes. We propose tests and methods to analyze convergence that may help detect these situations., Comment: 5 pages, RevTex, 4 figures

We present new many-parameter families of strongly and symmetric hyperbolic formulations of Einstein's equations that include quite general algebraic and live gauge conditions for the lapse. The first system that we present has 30 variables and incorporates an algebraic relationship between the lapse and the determinant of the three metric that generalizes the densitized lapse prescription. The second system has 34 variables and uses a family of live gauges that generalizes the Bona-Masso slicing conditions. These systems have free parameters even after imposing hyperbolicity and are expected to be useful in 3D numerical evolutions. We discuss under what conditions there are no superluminal characteristic speeds.

We study the stability properties of the Kidder-Scheel-Teukolsky (KST) many-parameter formulation of Einstein's equations for weak gravitational waves on flat space-time from a continuum and numerical point of view. At the continuum, performing a linearized analysis of the equations around flat spacetime, it turns out that they have, essentially, no non-principal terms. As a consequence, in the weak field limit the stability properties of this formulation depend only on the level of hyperbolicity of the system. At the discrete level we present some simple one-dimensional simulations using the KST family. The goal is to analyze the type of instabilities that appear as one changes parameter values in the formulation. Lessons learnt in this analysis can be applied in other formulations with similar properties., Comment: 17 Pages, RevTex, 16 figures

We discuss an equivalence between the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) formulation of the Einstein evolution equations, a subfamiliy of the Kidder--Scheel--Teukolsky formulation, and other strongly or symmetric hyperbolic first order systems with fixed shift and densitized lapse. This allows us to show under which conditions the BSSN system is, in a sense to be discussed, hyperbolic. This desirable property may account in part for the empirically observed better behavior of the BSSN formulation in numerical evolutions involving black holes., Comment: 7 Pages, RevTex, no figures

This is the first paper in a series aimed to implement boundary conditions consistent with the constraints' propagation in 3D numerical relativity. Here we consider spherically symmetric black hole spacetimes in vacuum or with a minimally coupled scalar field, within the Einstein-Christoffel symmetric hyperbolic formulation of Einstein's equations. By exploiting the characteristic propagation of the main variables and constraints, we are able to single out the only free modes at the outer boundary for these problems. In the vacuum case a single free modes exists which corresponds to a gauge freedom, while in the matter case an extra mode exists which is associated with the scalar field. We make use of the fact that the EC formulation has no superluminal characteristic speeds to excise the singularity. We present a second-order, finite difference discretization to treat these scenarios, where we implement these constraint-preserving boundary conditions, and are able to evolve the system for essentially unlimited times. As a test of the robustness of our approach, we allow large pulses of gauge and scalar field enter the domain through the outer boundary. We reproduce expected results, such as trivial (in the physical sense) evolution in the vacuum case (even in gauge- dynamical simulations), and the tail decay for the scalar field.