Evolutions of Nearly Maximally Spinning Black Hole Binaries Using the Moving Puncture Approach

TL;DR

This paper demonstrates that moving puncture numerical relativity codes can successfully evolve nearly maximally spinning black hole binaries, producing highly accurate waveforms with overlaps exceeding 0.999 compared to established simulations.

Contribution

It introduces a new evolution of nearly maximally spinning black hole binaries using the moving puncture approach with improved data and gauge conditions, achieving highly accurate waveforms.

Findings

Waveform overlap exceeds 0.999 with previous results

Successful evolution of black hole binaries with spin chi=0.99

Use of HiSpID data and modified lapse enhances stability

Abstract

We demonstrate that numerical relativity codes based on the moving punctures formalism are capable of evolving nearly maximally spinning black hole binaries. We compare a new evolution of an equal-mass, aligned-spin binary with dimensionless spin chi=0.99 using puncture-based data with recent simulations of the SXS Collaboration. We find that the overlap of our new waveform with the published results of the SXS Collaboration is larger than 0.999. To generate our new waveform, we use the recently introduced HiSpID puncture data, the CCZ4 evolution system, and a modified lapse condition that helps keep the horizon radii reasonably large.

Tables

Table 1. Table 1: Initial data parameters for a χ = 0.99 𝜒 0.99 \chi=0.99 highly spinning binary. The two spins are given by S → 1 , 2 = ( 0 , 0 , S ) subscript → 𝑆 1 2 0 0 𝑆 \vec{S}_{1,2}=(0,0,S) and the two momenta are P 1 , 2 = ± ( 0 , P , 0 ) subscript 𝑃 1 2 plus-or-minus 0 𝑃 0 P_{1,2}=\pm(0,P,0) . The parameter ℳ ℳ {\cal M} is the mass of the two black holes. Unlike for Bowen-York data, the momenta and spins cannot be specified exactly. However, the mass ℳ ℳ {\cal M} is very close to the measured horizon mass m H subscript 𝑚 𝐻 m_{H} . Quantities denoted by “init” were measured at t = 0 𝑡 0 t=0 , while quantities denoted by “equi” are averaged over the several orbits. A FE subscript 𝐴 FE A_{\rm FE} , s FE subscript 𝑠 FE s_{\rm FE} , r min subscript 𝑟 min r_{\rm min} , r max subscript 𝑟 max r_{\rm max} , and s far subscript 𝑠 far s_{\rm far} are attenuation parameters. m H subscript 𝑚 𝐻 m_{H} , S 𝑆 S , χ 𝜒 \chi are masses, spin angular momenta, and dimensionless spins, respectively, of the two black holes. M rem subscript 𝑀 rem M_{\rm rem} and χ rem subscript 𝜒 rem \chi_{\rm rem} are the remnant mass and dimensionless spin. The quantity r H subscript 𝑟 𝐻 r_{H} is the polar coordinate radius of the horizons. Finally, M ADM subscript 𝑀 ADM M_{\rm ADM} and J ADM subscript 𝐽 ADM J_{\rm ADM} are the ADM masses and spins.

$ℳ/M=0.505570$	$P/M=0.09675$
$S/M^2=0.253045$	$A_{\mathrm{FE}}=0.99$
$s_{\mathrm{FE}}=1.7$	$r_{\min }=0.01$
$r_{\max }=0.40$	$s_{\mathrm{far}}=10$
$J_{\mathrm{ADM}}/M^2=1.42621$	$M_{\mathrm{ADM}}/M=0.99998$
$m_{H\mathrm{init}}/M=0.50555$	$m_{H\mathrm{equi}}/M=0.5072\pm 0.0004$
$S_{\mathrm{init}}/M^2=0.2529$	$S_{\mathrm{equi}}/M^2=0.2547\pm 0.0004$
${\chi }_{\mathrm{init}}=0.9897$	${\chi }_{\mathrm{equi}}=0.9903\pm 0.0002$
$r_{H\mathrm{init}}/M=0.44$	$r_{H\mathrm{equi}}/M=0.082\pm 0.001$
$M_{\mathrm{rem}}/M=0.898\pm 0.001$	${\chi }_{\mathrm{rem}}=0.949\pm 0.001$

Equations

$${\gamma }_{ij}={\psi }^4{\tilde{\gamma }}_{ij},$$

$$K_{ij}={\psi }^{-2}{\tilde{A}}_{ij}+\frac{1}{3}K{\gamma }_{ij},$$

$${\tilde{A}}_{ij}={\tilde{M}}_{ij}+(\tilde{\mathbb{L}}b{)}_{ij},$$

$${\tilde{\gamma }}_{ij}={\tilde{\gamma }}_{ij}^{(+)}+{\tilde{\gamma }}_{ij}^{(-)}-{\delta }_{ij},$$

$$K=K^{(+)}+K^{(-)},$$

$$M_{ij}={\left[{\tilde{A}}_{ij}^{(+)}+{\tilde{A}}_{ij}^{(-)}\right]}^{\mathbf{TF}},$$

$${\psi }_0={\psi }_{(+)}+{\psi }_{(-)}-1,$$

$$r_{\mathrm{QI}}=r_{\mathrm{FE}}[1-A_{\mathrm{FE}}\exp (-r_{\mathrm{FE}}^2/{s_{\mathrm{FE}}}^2)],$$

$${\tilde{D}}^{2}u-g\frac{\psi \tilde{R}}{8}-g\frac{{\psi }^5{K}^2}{12}+g\frac{{\tilde{A}}_{ij}{\tilde{A}}^{ij}}{8{\psi }^7}+g{\tilde{D}}^2\left({\psi }_{(+)}+{\psi }_{(-)}\right)$$

$${\tilde{\Delta }}_{\mathbb{L}}{b}^i+g{\tilde{D}}_j{\tilde{M}}^{ij}-g\frac{2}{3}{\psi }^6{\tilde{\gamma }}^{ij}{\tilde{D}}_{j}K$$

$$g$$

$$g_{\pm }$$

$$\mathcal{G}(r_{\pm })$$

$${\tilde{\gamma }}_{ij}\to {\delta }_{ij}+g({\tilde{\gamma }}_{ij}-{\delta }_{ij}),$$

$${\tilde{\Gamma }}_{ij}^k\to g{\tilde{\Gamma }}_{ij}^k.$$

$${\tilde{\gamma }}_{ij}^{(\pm )}\to f(r_{\pm })({\tilde{\gamma }}_{ij}^{(\pm )}-{\delta }_{ij})+{\delta }_{ij},$$

$$K^{(\pm )}\to f(r_{\pm })K^{(\pm )},$$

$$\left({\psi }_{(\pm )}-1\right)\to f(r_{\pm })\left({\psi }_{(\pm )}-1\right),$$

$$({\partial }_t-{\beta }^i{\partial }_i)\alpha$$

$${\partial }_t{\beta }^a$$

$$\eta (\vec{r})=({\eta }_c-{\eta }_o)\exp (-(r/{\eta }_s{)}^4)+{\eta }_o,$$

$$\underset{t_0}{\mathrm{MAX}}\frac{\left|⟨R(t),S(t+t_0)⟩\right|}{\sqrt{\left|⟨R(t),R(t)⟩\right|\left|⟨S(t+t_0),S(t+t_0)⟩\right|}},$$

$$⟨a(t),b(t)⟩={\int }_0^{t_f}\bar{a}(t)b(t)dt,$$

Full text

I Introduction

Since the breakthroughs in numerical relativity of 2005 Pretorius (2005); Campanelli et al. (2006a); Baker et al. (2006) it is possible to accurately simulate moderate-mass-ratio and moderate-spin black-hole binaries. State of the art numerical relativity codes now routinely evolve binaries with mass ratios as small as $q\lesssim 1/20$ Gonzalez et al. (2009); Lousto et al. (2010); Lousto and Zlochower (2011); Sperhake et al. (2011); Chu et al. (2016); Jani et al. (2016), and are pushing towards much smaller mass ratios. Indeed, there have been several explorations of $q=1/100$ binaries Lousto and Zlochower (2011); Sperhake et al. (2011).

However, when it comes to highly-spinning binaries, prior to the work of Lovelace et al. (2008) of the SXS Collaboration 111https://www.black-holes.org, it was not even possible to construct initial data for binaries with spins larger than $\sim 0.93$ Cook and York (1990). This limitation was due to the use of conformally flat initial data. Conformal flatness is a convenient assumption because the Einstein constraint system takes on particularly simple forms. Indeed, using the puncture approach, the momentum constraints can be solved exactly using the Bowen-York ansatz Bowen and York (1980). There were several attempts to increase the spin of black hole, while still preserving conformal flatness Dain et al. (2002); Lousto et al. (2012), but these introduced negligible improvements. Lovelace et al. Lovelace et al. (2008) were able to overcome these limitations by choosing the initial data to be a superposition of conformally Kerr black holes in the Kerr-Schild gauge. Using these new data, they were soon able to evolve binaries with spins as large as 0.97 Lovelace et al. (2012) and, later, spins as high as 0.994 Scheel et al. (2015).

While spins of $0.92$ may seem reasonably close to 1, the scale is misleading. The amount of rotational energy in a black hole with spin 0.9 is only 52% of the maximum. Furthermore, particle limit and perturbative calculations show even more extreme differences between spins of 1 and spins only slightly smaller. For example, Yang et al. Yang et al. (2015) studied an analog to turbulence in black-hole perturbation theory. For spins close to 1, there is an inverse energy cascade from higher azimuthal ($m$) modes to lower ones for $\ell$ modes that obey $\epsilon=|1-\chi|\lesssim\ell^{-2}$. This gives hints that a more useful measure of the spin is actually $1/\epsilon$. Similarly, both the analysis of Kerr geodesics Schnittman (2014); Berti et al. (2015) and particle-limit calculations of recoils Hirata (2011); van de Meent (2014), indicate that the dynamics of nearly extremal-spin black holes cannot be elucidated with any degree of certainty using lower spin simulations.

Another area of interest is the use of numerical relativity waveforms in the detection and parameter estimation of gravitational wave signals as observed by LIGO and other detectors Abbott et al. (2016). This important region of parameter space with highly spinning binaries is currently poorly covered and will benefit from new and accurate simulations.

Recently, we introduced a version of highly-spinning initial data, also based on the superposition of two Kerr black holes Ruchlin et al. (2017); Healy et al. (2016), but this time in a puncture gauge. The main differences between the two approaches is how easily the latter can be incorporated into moving-punctures codes. In Ref. Ruchlin et al. (2017), we were able to evolve an equal-mass binary with aligned spins, and spin magnitudes of $\chi=0.95$, using this new data and compare with the results of the Lovelace et al..

Prior to our work, Hannam et al. Hannam et al. (2007) considered the case of non-boosted, highly spinning black holes. Similar to what we see here, they found that removing the conformally flat ansatz greatly reduces the amount of unphysical radiation.

In this paper, we show the results of a simulation of an equal-mass binary with aligned spins of $\chi=0.99$. We compare the $(\ell=2,m=2)$ and $(\ell=3,m=m2)$ modes of the waveform with those previously published by the SXS Collaboration in Scheel et al. (2015). This comparison allows us to assess the errors in these waveforms and to gain confidence about reaching the required accuracy for use in gravitational wave astronomy.

We use the following standard conventions throughout this paper. In all cases, we use geometric units where $G=1$ and $c=1$. Latin letters ($i$, $j$, $\ldots$) represent spatial indices. Spatial 3-metrics are denoted by $\gamma_{ij}$ and extrinsic curvatures by $K_{ij}$. The trace-free part of the extrinsic curvature is denoted by $A_{ij}$. A tilde indicates a conformally related quantity. Thus $\gamma_{ij}=\psi^{4}\tilde{\gamma}_{ij}$ and $A_{ij}=\psi^{-2}\tilde{A}_{ij}$, where $\psi$ is some conformal factor. We denote the covariant derivative associated with $\gamma_{ij}$ by $D_{i}$ and the covariant derivative associated with $\tilde{\gamma}_{ij}$ by $\tilde{D}_{i}$. A lapse function is denoted by $\alpha$, while a shift vector by $\beta^{i}$.

This paper is organized as follows. In Sec. II.1, we provide a brief overview of how the initial data are constructed. In Sec. II.2 we describe the numerical techniques used to evolve these data. In Sec. III, we compare the new HiSpID waveform with a similar SXS waveform. In Sec. III.1, we analyze the various diagnostics to determine the accuracy of the simulation. Finally, in Sec. IV, we discuss our results.

II Numerical Techniques

II.1 Initial Data

We construct initial data for a black-hole binary with individual spins $\chi_{1,2}=0.99$ using the HiSpID code Ruchlin et al. (2017); Healy et al. (2016). The HiSpID code solves the four Einstein constraint equations using the conformal transverse traceless decomposition York (1999); Cook (2000); Pfeiffer and York (2003); Alcubierre (2008). In this approach, the spatial metric $\gamma_{ij}$ and extrinsic curvature $K_{ij}$ are given by

[TABLE]

where the conformal metric $\tilde{\gamma}_{ij}$, the trace of the extrinsic curvature $K$, and the trace-free tensor $\tilde{M}_{ij}$ are free data. The Einstein constraints then become a set of four coupled elliptical equations for the scalar field $u=\psi-\psi_{0}$ and components of the spatial vector $b^{i}$ ($\psi_{0}$ is a singular function specified analytically). The resulting elliptical equations are solved using an extension to the TwoPunctures Ansorg et al. (2004) thorn.

The free data are chosen by superimposing two boosted Kerr black holes, as described in more detail in Ruchlin et al. (2017). The superposition has the form

[TABLE]

where $(+)$ and $(-)$ refer to the two black holes, $\tilde{\gamma}^{(\pm)}_{ij}$ and $\tilde{A}_{ij}$ are the conformal metric and trace-free extrinsic curvatures for a boosted and rotated Kerr black hole, $K^{(\pm)}$ is the mean curvature, and the conformal factor $\psi_{(\pm)}$ is chosen such that $\psi_{(\pm)}=\sqrt[12]{{\rm det}(\gamma_{ij}^{(\pm)})}$ (where $\gamma_{ij}^{(\pm)}$ is the physical metric from a boosted and rotated Kerr black hole).

To get $\tilde{\gamma}_{ij}^{(\pm)}$, etc., we start with Kerr black holes in quasi-isotropic (QI) coordinates and perform a fisheye (FE) radial coordinate transformation (where $r_{\rm QI}=0$ is the location of the puncture),

[TABLE]

where $r_{\rm FE}$ is the fisheye radial coordinate, $r_{\rm QI}$ is the original QI radial coordinate, and $A_{\rm FE}$ and $s_{\rm FE}$ are parameters. These coordinates have the property that at large $r_{\rm FE}$, $r_{\rm QI}\approx r_{\rm FE}$, and at small $r_{\rm FE}$, $dr_{\rm QI}=(1-A)dr_{\rm FE}$ (i.e., $dr_{\rm QI}<dr_{\rm FE}$). The FE transformation is needed because it expands the horizon size from $r_{h}\approx 0.035$ to $r_{h}\approx 0.5$. We then transform the metric to Cartesian-like coordinates of the form $x=r\sin\theta\cos\phi$, $y=r\sin\theta\sin\phi$, $z=r\cos\theta$, where $r=r_{\rm FE}$. We then perform a Lorentz-like boost on this metric and, in the case of nonaligned spins, a rotation. The resulting 4-metric is then decomposed into a spatial metric $\gamma_{ij}$ and extrinsic curvature $K_{\ij}$.

We use the $g$ attenuation described in Ruchlin et al. (2017) to modify both the metric and elliptical equations inside the horizons. We briefly summarize the procedure here. The modified Hamiltonian and momentum constraint equations for the correction functions $u$ and $b^{i}$ are

[TABLE]

where $\tilde{\Delta}_{\mathbb{L}}b^{i}\equiv\tilde{D}_{j}(\tilde{\mathbb{L}}b)^{ij}$ is the vector Laplacian and $\tilde{R}$ is the scalar curvature associated with $\tilde{\gamma}_{ij}$, and where the attenuation function $g$ takes the form

[TABLE]

$r_{\pm}$ is the coordinate distance to puncture $(+)$ or $(-)$, and the parameters $r_{\rm min}<r_{\rm max}$ are chosen to be within the horizon. In addition, we attenuate the background metric itself when calculating the $\tilde{D}^{2}u$ and $\tilde{\Delta}_{\mathbb{L}}b^{i}$. To do this, we take

[TABLE]

Note that the modified $\tilde{\Gamma}^{k}_{\,ij}$ is not consistent with the modified $\tilde{\gamma}_{ij}$. There is no advantage to making them consistent because the constraints will be violated in the attenuation zone regardless. By modifying the metric in this way, we can ensure that the elliptical system has exactly the form of the flat space Poisson system in the vicinity of the punctures.

Finally, far from the holes, we attenuate $\tilde{\gamma}_{ij}$, $K$, and $\psi_{0}$. This is achieved by consistently changing the metric fields and their derivatives so that

[TABLE]

where $f(r)=\exp(-r^{4}/s_{\rm far}^{4})$ and $r_{\pm}$ is the coordinate distance to puncture $(+)$ or $(-)$.

For compatibility with the original TwoPunctures code, we chose to set up HiSpID so that the parameters of the binary are specified in terms of momenta and spins of the two holes. However, unlike for Bowen-York data, the values specified are only approximate, as the solution vector $b^{i}$ can modify both of these. In practice, we find that the spins are modified by only a trivial amount while orbital angular momentum is reduced significantly. We compensate for this by choosing larger momentum parameters than those predicted by simple quasicircular conditions would imply Healy et al. (2017). All parameters for the $\chi=0.99$ run are given in Table 1. The quantity $r_{H}$ in the table is the polar coordinate radius (which is the smallest radius on each horizon). As this is gauge dependent, it can change arbitrarily during the evolution. However, large changes are generally undesirable. The size of $r_{H}$ is also directly related to the number of refinement levels required, and therefore to the computational cost. An ideal gauge would have $r_{H}$ settle to a moderate value and remain there. The initial size of the horizon is chosen to be large in order to speed up the convergence of the initial data solver (this is due to the scale set by the $g$ attenuation discussed above). However, the gauge conditions we use quickly drive $r_{H}$ towards smaller values. We note that in quasi-isotropic coordinates, the coordinate radius of a maximally spinning black hole is zero.

II.2 Evolution

We evolve black hole binary initial data sets using the LazEv Zlochower et al. (2005) implementation of the moving punctures approach for the conformal and covariant formulation of the Z4 (CCZ4) system (Ref. Alic et al. (2012)) which includes stronger damping of the constraint violations than the standard BSSNOK Nakamura et al. (1987); Shibata and Nakamura (1995); Baumgarte and Shapiro (1998) system. For the run presented here, we use centered, eighth-order accurate finite differencing in space Lousto and Zlochower (2008) and a fourth-order Runge-Kutta time integrator. Our code uses the Cactus/EinsteinToolkit cac ; ein infrastructure. We use the Carpet mesh refinement driver to provide a “moving boxes” style of mesh refinement Schnetter et al. (2004). Fifth-order Kreiss-Oliger dissipation is added to evolved variables with dissipation coefficient $\epsilon=0.1$. For the CCZ4 damping parameters, we chose $\kappa_{1}=0.2$, $\kappa_{2}=0$, and $\kappa_{3}=0$ (see Alic et al. (2012)), but found that these had to be modified during the evolution.

We locate the apparent horizons using the AHFinderDirect code Thornburg (2004) and measure the horizon spins using the isolated horizon (IH) algorithm Dreyer et al. (2003). We calculate the radiation scalar $\psi_{4}$ using the Antenna thorn Campanelli et al. (2006b); Baker et al. (2002). We then extrapolate the waveform to an infinite observer location using the perturbative formulas given in Ref. Nakano et al. (2015).

For the gauge equations, we use Alcubierre et al. (2003); Campanelli et al. (2006a); van Meter et al. (2006)

[TABLE]

Note that the lapse is not evolved with the standard 1+log form. Here we multiply the rhs of the lapse equation by an additional factor of $\alpha$. This has the effect of increasing the equilibrium (coordinate) size of the horizons. For the initial values of shift, we chose $\beta^{i}(t=0)=0$, while for the initial values of the lapse, we chose an ad-hoc function $\alpha(t=0)=\tilde{\psi}^{-2}$, where $\tilde{\psi}=1+{\cal M}/(2r_{1})+{\cal M}/(2r_{2})$ and $r_{i}$ is the coordinate distance to BH $i$. For the function $\eta$, we chose

[TABLE]

where $\eta_{c}=2.0/M$, $\eta_{s}=40.0M$, and $\eta_{o}=0.25/M$. With this choice, $\eta$ is small in the outer zones. As shown in Ref. Schnetter (2010), the magnitude of $\eta$ limits how large the timestep can be with $dt_{\rm max}\propto 1/\eta$. Since this limit is independent of spatial resolution, it is only significant in the very coarse outer zones where the standard Courant-Friedrichs-Lewy condition would otherwise lead to a large value for $dt_{\rm max}$.

We evolved the $\chi=0.99$ data using 11 levels of refinement, with the outermost grid extending to $400M$ with a gridspacing of $2.78M$. The gridspacing on the finest grid was $h=M/368.64$. The total cost of the simulation was 710 KSU.

One remarkable consequence of these superimposed Kerr data is how small the initial pulse of unphysical radiation is. As first seen in the nonboosted case by Hannam et al. Hannam et al. (2007), the initial pulse is roughly four times as large as the orbital signal at a separation of $D=10M$. While this may sound quite large, for a $\chi=0.9$ binary, the amount of unphysical radiation for a Bowen-York binary is six times more, and it rapidly increases with spin. The full waveform, including initial pulse, is shown in Fig. 1.

III Results

We performed a single simulation from a coordinate separation of $10M$ (proper separation of $12.2M$) through merger for an equal-mass binary where both spins are aligned with the orbital angular momentum and have dimensionless magnitudes of 0.99. We compare these with the BBH:0177 waveform Scheel et al. (2015); SXS . In order to compare the HiSpID and SXS waveform, we rescale the time coordinate by the ratio of the final masses and then introduce a constant phase and time translation to minimize the RMS difference between the two waveforms. We note that the SXS waveform is longer by about $t=5000M$.

In Figs. 2 and 3, we directly compare the new HiSpID waveform with the corresponding SXS waveform. We translate the HiSpID waveform to maximize the overlap. Hence, the “starting” time in the figures is $t\sim 5000M$. The overlap is defined by Campanelli et al. (2009)

[TABLE]

where

[TABLE]

an overbar denotes complex conjugation, and $t_{0}$ is chosen to maximize the result, while $t=0$ corresponds to the time just after the initial pulse has radiated away and $t_{f}$ to the last timestep in the HiSpID simulation. We find an overlap of $0.99975$ for the ($\ell=2,m=2)$ mode, which is quite good considering that the HiSpID waveform is eccentric ($e\sim 0.01$), while the SXS waveform is not. Note the phase agreement is within 0.25 rad across the entire waveform and the amplitude agreement is better than 4%. The agreement in frequency is even better, with a relative difference of lass than 2% across the entire waveform.

The next largest modes after the $(\ell=2,m=\pm 2)$ modes are the $(\ell=4,m=\pm 4)$ modes. However, in our simulations, these show significant effects of dissipation postmerger. We therefore compare the $(\ell=3,m=2)$ mode instead. As shown in Fig. 4, the agreement between HiSpID and SXS is quite good even for a higher-order mode. The overlap between the $(\ell=3,m=2)$ modes is 0.998 [the constant $t_{0}$ was fixed by maximizing the overlap of the $(\ell=2,m=2)$ modes, the maximum overlap of the $(\ell=3,m=2)$ is 0.9998].

III.1 Diagnostics

One of the most important diagnostics for a BHB simulation is the degree to which the constraints are satisfied and to what degree the horizon masses and spins are conserved. In Fig. 5, we show the individual horizon mass and (dimensionless spin). Note that prior to merger, the spins are within $\pm 0.001$ of $0.99$ and the masses change by less than 0.2%. In Fig. 6, we show the $L^{2}$ norm of the Hamiltonian and momentum constraints. Here the $L^{2}$ norm is over the region outside the two horizons (or common horizon) and inside a sphere of radius $30M$. Note how the constraints start small ($10^{-8}$) and quickly increase to $10^{-4}$. This increase is due to unresolved features in the initial data (i.e., the AMR grid cannot propagate high-frequency data accurately). The constraints then damp, as is expected for CCZ4. However, they start to exponentially blow up around $400M$. We found that the parameters $\kappa_{1}$ and $\kappa_{2}$ had to be fine-tuned to prevent this blow-up. We found that increasing the damping parameters can effectively drive the constraints smaller for a short time, but large values of $\kappa_{i}$ led to an exponential blowup of the constraints at later times. We used a trial-and-error approach to fine-tuning these parameters during the run. We show the values of $\kappa_{1}$ used during the evolution in the top of Fig. 6.

One challenge with the HiSpID data is obtaining low-eccentricity data without performing an iterative procedure where the initial data are evolved for a few orbits and then refined based on the measured orbital evolution Pfeiffer et al. (2007); Buonanno et al. (2011); Purrer et al. (2012); Buchman et al. (2012). In Healy et al. (2017), it was shown that relatively low-eccentricity initial data parameters can be obtained using higher-order post-Newtonian approximations. However, as shown in Table 1, unlike for Bowen-York data, here we cannot specify the initial momenta precisely. That is to say, the orbital angular momentum of the background (i.e., prior to the inclusion of corrections due to the fields $u$ and $b^{i}$) is significantly larger than the final orbital angular momentum of the initial data. We compensate for this by increasing the momentum parameters until the ADM angular momentum matches the expected value based on quasicircular orbits. However, we have no method of correcting for the radial momentum (other than using an iterative evolution procedure). Consequentially, the eccentricity of the initial data is relatively high at $e\approx 0.01$, as shown in Fig. 7. Of course, we can run the data for a few orbits and then refine the parameters, but such a procedure is computationally expensive. We are thus working on improving the evolution efficiency.

One method which we found was useful for increasing the run speed was to change the lapse condition. When using the standard 1+log lapse, the horizons are a factor of 0.625 as wide (see Fig. 8). Evolving the data with horizons this small requires roughly a factor of 2 more in terms of computational expense because an additional level of refinement is needed. Using harmonic slicing leads to still larger horizons, but this proved to be unstable. The rapid change in the gauge at early time, as is evident in the size and shape of the horizon (see Fig. 8) may be responsible for the initial jump in the constraint violations seen in Fig. 6.

IV Discussion

In this paper we demonstrated that it is possible to evolve black hole binaries with nearly maximal spin using the “moving puncture” formalism. This means that comparative studies of these challenging evolutions by the two main methods (the generalized harmonic approach used by SXS and various flavors of the “moving punctures” approach used by many other groups) to numerically solve the field equations of general relativity field equations can now be performed. Independent comparison, along the lines explored in Lovelace et al. (2016), have been very successful in demonstrating the accuracy and correctness of moderate-spin black hole simulations. These new techniques also open the possibility to explore a region of parameter space which is of high interest for both astrophysical and gravitational wave studies.

Acknowledgements.

The authors thank M.Scheel for careful reading of the manuscript and gratefully acknowledge the National Science Foundation (NSF) for financial support from Grants No. PHY-1607520, No. PHY-1707946, No. ACI-1550436, No. AST-1516150, No. ACI-1516125. This work used the Extreme Science and Engineering Discovery Environment (XSEDE) [allocation TG-PHY060027N], which is supported by NSF grant No. ACI-1548562. Computational resources were also provided by the NewHorizons and BlueSky Clusters at the Rochester Institute of Technology, which were supported by NSF grants No. PHY-0722703, No. DMS-0820923, No. AST-1028087, and No. PHY-1229173.