To fully unlock the scientific potential of upcoming gravitational wave GW interferometers numerical relativity NR simulation accuracy will need to be greatly enhanced. We present three infrastructureagnostic improvements to the movingpuncture approach for binary black hole BBH simulations aimed at greatly reducing constraint violation and improving GW predictions. Although these improvements were developed within the highly efficient NR code BlackholesHome we demonstrate their effectiveness in the widely adopted Einstein TOOLKITCarpet adaptive mesh refinement framework. Our improvements include a modified KreissOliger dissipation prescription a Hamiltonian constraint damping adjustment to the BaumgarteShapiroShibataNakamura BSSN equations and an extra term to the inlineformulammlmath displayinlinemmlmrowmmlmn1mmlmnmmlmommlmommlmilogmmlmimmlmrowmmlmathinlineformula lapse evolution equation that slows the development of the sharp lapse feature which dominates numerical errors in BBH simulations. With minimal increase in computational cost these improvements greatly reduce GW noise enabling the extraction of highorder GW modes previously obscured by numerical noise. They also improve convergence properties near and inside the convergent regime reduce Hamiltonian momentum constraint violations in the strongfield region by roughly 2 3 orders of magnitude and in the GWextraction zone by 5 2 orders of magnitude. To promote community adoption we have opensourced the improved Einstein TOOLKIT thorn BaikalVacuum used in this work. Although our focus is on BBH evolutions and the BSSN formulation these improvements may also benefit compact binary simulations involving matter and other formulations a focus for future investigations.
['General relativity', 'alternative theories of gravity', 'General Relativity and Quantum Cosmology', 'Astrophysics - High Energy Astrophysical Phenomena']
Improved movingpuncture techniques for compact binary simulations
{'Improved Moving-Puncture Techniques for Compact Binary Simulations': 'Zachariah B. Etienne 1, 2, 3, ∗  \n1  \nDepartment of Physics, University of Idaho, Moscow, ID 83843, USA 2 Department of Physics and Astronomy, West Virginia University, Morgantown, WV 26506, USA 3 Center for Gravitational Waves and Cosmology, West Virginia University, Chestnut Ridge Research Building, Morgantown, WV 26505, USA (Dated: September 17, 2024)  \nTo fully unlock the scientific potential of upcoming gravitational wave (GW) interferometers, numerical relativity (NR) simulation accuracy will need to be greatly enhanced. We present three infrastructure-agnostic improvements to the moving-puncture approach for binary black hole (BBH) simulations, aimed at greatly reducing constraint violation and improving GW predictions. Although these improvements were developed within the highly efficient NR code BlackHoles@Home , we demonstrate their effectiveness in the widely-adopted Einstein Toolkit / Carpet AMR framework. Our improvements include a modified Kreiss-Oliger dissipation prescription, a Hamiltonian-constraintdamping adjustment to the BSSN equations, and an extra term to the 1+log lapse evolution equation that slows the development of the sharp lapse feature, which dominates numerical errors in BBH simulations. With minimal increase in computational cost, these improvements greatly reduce GW noise, enabling the extraction of high-order GW modes previously obscured by numerical noise. They also improve convergence properties near and inside the convergent regime, reduce Hamiltonian (momentum) constraint violations in the strong-field region by roughly two (three) orders of magnitude, and in the GW-extraction zone by five (two) orders of magnitude. To promote community adoption, we have open-sourced the improved Einstein Toolkit thorn BaikalVacuum used in this work. Although our focus is on BBH evolutions and the BSSN formulation, these improvements may also benefit compact binary simulations involving matter and other formulations, a focus for future investigations.', 'I. INTRODUCTION': "Numerical relativity (NR) simulations of compact binary mergers play a crucial role in gravitational wave (GW) data analysis, enabling reliable estimation of physical parameters from the late-inspiral, merger, and postmerger phases of GW signals. To date, detector noise in GW observations has generally surpassed the numerical errors intrinsic to NR simulations and GW approximants derived from them, so GW parameter estimation errors have generally been dominated by noise in the GW signal. However, improved sensitivities offered by current and forthcoming GW detectors threaten to change this dynamic. Third-generation (3G) detectors will necessitate a roughly tenfold improvement in NR simulation accuracy [1, 2]. Therefore, enhancing the accuracy of NR codes is vital for unlocking the full scientific potential of GW astronomy moving forward.  \nImproved accuracy can be achieved through various means. One widely adopted approach involves building next-generation NR codes upon more scalable infrastructures. Increased scalability enables simulations to be performed at higher resolutions and with greater physical realism. While running simulations at higher resolutions generally reduces numerical errors, it does not directly address one of the primary sources of error in existing NR codes: noise associated with the selection of numerical grids.  \nCartesian adaptive mesh refinement (AMR) grids, the most popular gridding approach in NR, span the many decades in physical scales intrinsic to compact binaries with nested Cartesian grid patches. The highest resolution patches are placed where the spacetime fields are sharpest: at the compact objects themselves. These patches are embedded in increasingly coarser patches, extending to the outer boundary. Cartesian AMR grids are both numerically convenient and highly scalable with modern AMR algorithms.  \nWhen it comes to solving Einstein's equations in NR, the moving-puncture approach [3, 4] is the most popular method. This approach generally combines the BSSN [5, 6], CCZ4 [7] or Z4c [8] 3+1 decompositions of Einstein's equations with the moving-puncture gauge conditions [3, 4, 9], enabling the evolution of black holes (BHs) that traverse the numerical grid without the need to excise their interiors. However, mixing the moving-puncture technique with Cartesian AMR grids results in significant numerical noise that negatively impacts the accuracy of compact binary simulations.  \nAs shown in Fig. 1, during the initial phase of a movingpuncture binary black hole (BBH) evolution, the lapse rapidly settles to new values in the strong-field region, releasing a sharp, superluminal gauge wave packet that radiates outward from the BHs (top panel) at a speed of √ 2 c -the characteristic speed of the 1+log lapse [10] and the only superluminal characteristic speed [11]. Given our standard choice of initial ('precollapsed') lapse, in the weak field the lapse tends to a value of 1 -α ∝ r -1 , both before and after the wave packet passes. Thus away  \n<!-- image -->  \nFIG. 1. Development of sharp moving-puncture lapse feature and its aftermath. Top left : r (1 -α ) colormap on the orbital plane, for r > 25 M ; data for r ≤ 25 M were excised (represented by a black dot). Bottom left : r (1 -α ) data along the x -axis, for x > 25 M . Top right : H constraint violation, plotted as log 10 |H| on the orbital plane. The dashed orange curve corresponds to the data in the top-left plot, and the blue solid curve corresponds to the data at t = 15 . 1 M later. The black arrow denotes the propagation distance and direction of a wave propagating outward at √ 2 c , from t/M = 52 . 1 to t/M = 67 . 2. Bottom right : log 10 |H| along the x -axis at t = 67 . 2 M , for x > 25 M .  \n<!-- image -->  \nfrom the wave packet, the r (1 -α ) plotted in Fig. 1 tends to constant values in the weak field.  \nAs the sharp lapse wave packet travels outward from the strong-field region, it crosses multiple AMR refinement boundaries. At each boundary, the numerical resolution drops discontinuously by a factor of two [12]. This crossing partially reflects the wave and generates substantial numerical noise in its wake.  \nAlthough this is a gauge wave, errors in its evolution directly contaminate the evolution of physical quantities, particularly the extrinsic curvature, which depends on second-order derivatives of the lapse. As Einstein's equations are highly nonlinear, errors in the lapse's evolution contaminate all physical quantities, including the Einstein constraints. The right panel of Fig. 1 illustrates the roughly seven-decade jump in the Hamiltonian-constraint violation H from its initial values following the passage of this wave. Elevated levels of H violation persist near refinement boundaries, evident as square shapes in the colormap.  \nMeanwhile, in the strong-field region, H violation acts as an effective energy density, directly influencing the dynamics of the compact binary and, thus, one of the core outputs of the NR simulation: the GW prediction. Regarding GW predictions, noise appears as spurious  \nhigh-frequency oscillations in the outgoing-radiation Weyl scalar ψ 4 . ψ 4 is decomposed into spin-weight s = -2 spherical harmonics and integrated twice in time to compute the GW strain. While the noise is evident in the ( ℓ = m = 2) dominant mode of ψ 4 , it completely obscures modes with ℓ > 4. Such high-order modes may be detectable by next-generation GW observatories, even in near-symmetric binaries [13]. Thus minimizing this noise is crucial for precise BBH waveform predictions for upcoming GW detectors.  \nWe emphasize that this discussion has been limited to Cartesian AMR grids not set to refine at the sharp gauge features, as this approach is widely used in the field. Improved refinement algorithms or grid structures may result in significantly less constraint violation and GW noise than what is reported here.  \nThis paper introduces improvements to the movingpuncture approach, aimed at curtailing the development of sharp features in the lapse and mitigating their impact. Our first two improvements focus on decreasing the strength of low-pass filters near punctures-where sharp fields are desirable-and increasing the strength far away.  \nThe first improvement multiplies the strength parameter of Kreiss-Oliger (KO) dissipation by e -4 ϕ , where ϕ is the standard BSSN conformal factor, and increases its  \nstrength from typical values by 2 to 5 times.  \nThe second improvement modifies a powerful but not widely used low-pass parabolic filter, which adds a term proportional to H to the right-hand side of the ϕ evolution equation. In the same vein as [14], we set the strength of the filter based on the maximum value allowed by the CFL condition, according to the local numerical grid spacing. This enables us to increase its strength by more than an order of magnitude on coarser numerical grids in the region far from the punctures.  \nOur third improvement adds a term to the 1+log lapse evolution equation that slows the evolution of α in the area immediately outside the BH horizon. This arrests the formation of the sharp lapse wave packet, both reducing its magnitude and greatly smoothing it.  \nModifying the moving-puncture approach in CartesianAMR BBH evolutions to mitigate the impact of the sharp lapse feature is not a new idea; [10] introduced several techniques to this end, achieving about a 40% reduction in GW noise (in ψ 4 ) and about 1.3 orders of magnitude reduction in |H| violations outside the horizons. While these improvements were effective, their recipes depend on careful tuning to the specific grid structure.  \nWhile we adopt the BSSN formalism (like [10]) and employ a constraint-damping term, the moving-puncture approach is also widely used within the CCZ4 and Z4c formalisms, which contain constraint-damping terms. [15] performed evolutions of puncture BBHs in the Z4c formalism, reporting 1-2 order of magnitude reductions in Hamiltonian constraint violation outside the punctures, and GW phase and amplitude errors up to a factor of two smaller as compared to the standard BSSN movingpuncture approach.  \nBy contrast, the three improvements presented in this work, when combined, reduce Hamiltonian (momentum) constraint violations in the strong-field region by about two (three) orders of magnitude and in the GW-extraction zone by five (two) orders of magnitude. Further, constraint violation convergence is generally improved, implying that these gaps widen with increased resolution. Coordinate eccentricity is slightly decreased, and numerical noise in the dominant ℓ = m = 2 mode of ψ 4 is reduced by an average factor of 4.3 across resolutions. Fidelity of higher-order modes is greatly improved, enabling the ℓ = m = 6 mode, previously obscured by noise, to become visible.  \nAlthough our focus is primarily on the impact of these improvements on BBH evolutions with the BSSN formalism, we expect that they will also be beneficial in other compact binary evolutions and formalisms that use the moving-puncture technique. Determining whether our improved techniques significantly improve the fidelity of binary neutron star (BNS) simulations will be a focus of an upcoming paper. In addition, formulations like Z4c have proven quite effective in reducing constraint violations in BNS simulations [15]. Comparing BNS evolutions performed with our improved BSSN moving-puncture approach to Z4 techniques widely used across the community  \nwill also be a focus of future work.  \nThe rest of the paper is structured as follows. Section II introduces the basic equations, including the BSSN evolution, moving-puncture gauge, and constraint equations. Building on this, Sec. III presents prescriptions for all three of our improved techniques. Next, Sec. IV describes our algorithmic approach for simulating the BBH mergers presented in this work. Results from these simulations, discussed in Sec. V, demonstrate the impacts of our improved techniques. Finally, we conclude and describe avenues for future work in Sec. VI.", 'II. BASIC EQUATIONS': 'We adopt G = c = 1 units, and set our mass scale M = 1 to be the sum of the individual puncture masses reported by our initial data solver TwoPunctures . In these units the total ADM mass is 0.989946. Further we adopt standard Einstein summation rules on repeated indices, and the convention that Greek and Latin indices represent 4- and 3-dimensional quantities, respectively.', 'A. BSSN Evolution Equations': "In this paper, we investigate the numerical evolutions of BBHs using the covariant 'Lagrangian' BSSN formulation in vacuum, as presented in [16]. This formulation has been adapted for spherical coordinate reference metrics by [17] and for more general reference metrics by [18].  \nStarting with the ADM line element,  \nds 2 = -α 2 dt 2 + γ ij ( dx i + β i dt ) ( dx j + β j dt ) , (1)  \nγ ij , β i , and α denote the spatial 3-metric, shift vector, and lapse, respectively.  \nWe then define the conformal metric ¯ γ ij as  \n¯ γ ij = ( ¯ γ γ ) 1 3 γ ij = e -4 ϕ γ ij , (2)  \nwhere ϕ represents the conformal factor. Instead of evolving ϕ , we follow previous works (e.g., [19, 20]) and evolve W = e -2 ϕ , as W exhibits greater smoothness near puncture BHs, resulting in reduced truncation errors (also see discussion in [3]). To this end, derivatives of the conformal factor ϕ appearing in the BSSN evolution equations are calculated with respect to the evolved variable W via ∂ µ W = -2 e -2 ϕ ∂ µ ϕ ⇒ ∂ µ ϕ = -∂ µ W/ (2 W ), and so forth for second derivatives.  \nWe also introduce a reference metric ˆ γ ij , defined by  \n¯ γ ij = ˆ γ ij + ε ij . (3)  \nIn this work, the reference metric corresponds to the Cartesian metric ˆ γ ij = δ ij , which remains constant in both space and time. Thus, all derivatives of the conformal metric ¯ γ ij are evaluated directly with respect to ε ij ,  \nsince ∂ µ ¯ γ ij = ∂ µ ε ij . Although analytically equivalent, computing finite-difference derivatives of diagonal metric components with respect to ε ij instead of ¯ γ ij = 1 + ε ij results in reduced roundoff error in the weak-field region.  \nBuilding on these definitions, we adopt the following time-evolution equations, replacing all occurrences of ϕ and its derivatives with equivalent expressions in terms of W :  \n∂ t ε ij = L β ε ij + 2 3 ¯ γ ij ( α ¯ A k k -¯ D k β k ) -2 α ¯ A ij (4) ∂ t ¯ A ij = L β ¯ A ij -2 3 ¯ A ij ¯ D k β k -2 α ¯ A ik ¯ A k j + α ¯ A ij K + e -4 ϕ { -2 α ¯ D i ¯ D j ϕ +4 α ¯ D i ϕ ¯ D j ϕ +4 ¯ D ( i α ¯ D j ) ϕ -¯ D i ¯ D j α + α ¯ R ij } TF (5)  \n∂ t K = L β K + 1 3 αK 2 + α ¯ A ij ¯ A ij -e -4 ϕ ( ¯ D i ¯ D i α +2 ¯ D i α ¯ D i ϕ ) (7)  \n∂ t ϕ = L β ϕ + 1 6 ( ¯ D k β k -αK ) (6)  \n∂ t ¯ Λ i = L β ¯ Λ i + ¯ γ jk ˆ D j ˆ D k β i + 2 3 ∆ i ¯ D j β j + 1 3 ¯ D i ¯ D j β j -2 ¯ A ij ( ∂ j α -6 ∂ j ϕ ) +2 ¯ A jk ∆ i jk 4 α ¯ γ ij ∂ j K (8)  \n-3  \nwhere the TF superscript indicates the trace-free part, and L β denotes the Lie derivative along the shift vector β i . The covariant derivatives with respect to the reference metric ˆ γ ij and the barred spatial 3-metric ¯ γ ij are represented by ˆ D j and ¯ D j , respectively. The tensor ∆ i jk is formed from the difference between the barred and hatted Christoffel symbols:  \n∆ i jk = ¯ Γ i jk -ˆ Γ i jk ,  \nwhere ∆ i is defined ∆ i ≡ ¯ γ jk ∆ i jk . The conformal, tracefree extrinsic curvature, ¯ A ij , is given by:  \n¯ A ij = e -4 ϕ ( K ij -1 3 γ ij K ) ,  \nwhere K represents the trace of the extrinsic curvature K ij .", 'B. Moving-Puncture Gauge Evolution Equations': 'To close the system of time-evolution equations, we choose the moving-puncture gauge conditions [3, 4, 9], which include the 1+log slicing condition for the lapse and the Γ-driving (here, Λ takes the place of Γ) condition for the shift:  \n∂ t α = β j ∂ j α -2 αK (9)  \n∂ t β i = β j ∂ j β i + B i (10)  \n∂ t B i = β j ∂ j B i + 3 4 ( ∂ t ¯ Λ i + β j ∂ j ¯ Λ i ) -ηβ i , (11)  \nwhere the Γ-driving shift damping parameter is set in this paper to η = 1 /M .  \nVariations in moving-puncture gauge conditions, such as a first-order implementation of the shift condition, do exist. Although we did not explore these alternatives in our study, we direct the reader to van Meter et al . [11] for a comprehensive analysis of different implementations. Our choice of gauge conditions was primarily motivated by their finding that incorporating advection terms into the lapse and shift conditions leads to well-behaved characteristic speeds.', 'C. BSSN Constraint Equations': 'Analytically (as opposed to numerically), the BSSN Hamiltonian constraint H = 0, so any expression for H can be multiplied by an arbitrary prefactor and still remain perfectly valid. As such, two expressions for the BSSN Hamiltonian constraint H are commonly found in the literature, differing only in their prefactors.  \nIn this work, we adopt H with a negative prefactor in front of ¯ D 2 ϕ , its principal part. To differentiate our definition of H from other choices, we will refer to it as H -when necessary:  \nH -= 2 3 K 2 -¯ A ij ¯ A ij + e -4 ϕ ( ¯ R -8 ¯ D i ϕ ¯ D i ϕ -8 ¯ D 2 ϕ ) = H = 0 . (12)  \nThe second convention (e.g., [21]) multiplies both sides by -e 5 ϕ / 8. This results in a sign flip on H , and thus the prefactor on ¯ D 2 ϕ becomes positive. We therefore refer to this form as H + :  \nH + = -e 5 ϕ 8 H -. (13)  \nAlthough the two expressions H + = H -= 0 are equivalent, it is important to distinguish between them. As we will soon see, the specific implementation of Hamiltonianconstraint damping (Sec. III B) will depend on the choice of this prefactor.  \nRegarding the physical implications of constraint violation, if we write H -in its general form with the matter source term included:  \nH m -= H --16 πρ = 0 , (14)  \nthen we may conclude that, in vacuum evolutions, positive Hamiltonian-constraint violations H -behave as a positive effective energy density ρ eff , and negative violations as a negative ρ eff :  \nH -= 16 πρ eff . (15)  \nAs for the momentum constraint, we adopt the following form:  \nM i = e -4 ϕ ( ˆ D j ¯ A ij +2 ¯ A k ( i ∆ j ) jk +6 ¯ A ij ∂ j ϕ -2 3 ¯ γ ij ∂ j K ) = 0 , (16)  \nand define the scalar M 2 = ¯ γ ij M i M j , which we use for diagnostics in this work.  \nAs Hamiltonian- and momentum-constraint violation are defined at each of the > ∼ 10 7 points in our computational domain, for convenience we will often compute them as a volume-averaged L 2 norm:  \n||H|| = √ ∫ V H 2 d 3 x ∫ V d 3 x (17)  \nwhere the volume V will generally exclude small spheres around the BHs and be tailored to either a strong-field or weak-field region.  \n√ ||M 2 || = 4 √ ∫ V (¯ γ ij M i M j ) 2 d 3 x ∫ V d 3 x , (18)', 'III. IMPROVED TECHNIQUES': "We present three techniques that lower numerical errors in moving-puncture simulations, with minimal additional computational cost. These include curvature-adjusted Kreiss-Oliger dissipation (CAKO), coarse-grid-adjusted Hamiltonian-constraint damping (CA H D), and the slowstart lapse (SSL) technique.  \nThese improvements are meant to be both easy to implement and agnostic to numerical grid structures. To this end, while we demonstrate their effectiveness in the widely used Einstein Toolkit [22, 23] Carpet [24, 25] AMR infrastructure, the new techniques were originally devised and tuned during the development of a new, superefficient NR code called BlackHoles@Home .  \nIllustrated in Fig. 2, BlackHoles@Home adopts a unique multi-patch grid structure that combines grids in sphericallike and Cartesian-like coordinates. Einstein's equations are solved in the corresponding spherical-like and Cartesian-like bases, leveraging recent developments in covariant BSSN [16], BSSN with a reference metric [17, 18], and NR code generation with NRPy [18, 26].  \nWith these new grid structures, BlackHoles@Home performs full, 3+1 NR BBH evolutions using about 3GB of RAM or less, enabling its use on consumer-grade desktop computers. As such, it aims to generate the largest-ever full-NR BBH GW catalog through an upcoming volunteercomputing project of the same name.  \nOwing to BlackHoles@Home 's high efficiency, the improvements presented here were tuned over approximately 10 simulations of a single BBH scenario, and the tuned values were proven similarly effective across more than 100 different BBH scenarios simulated with BlackHoles@Home . These scenarios spanned mass ratios from q = 1 to q = 6 and included arbitrarily directed spin vectors with dimensionless magnitudes up to 0.5.  \nBuilding on the success of these BlackHoles@Home -tuned improvements, this work focuses on evaluating their performance on Cartesian AMR grids. To this end, we implement the improvements within our BaikalVacuum  \nFIG. 2. Schematic of BBH inspiral numerical meshes in BlackHoles@Home , the code used to develop and tune the improvements presented here.  \n<!-- image -->  \nEinstein Toolkit thorn, and release the improved version of BaikalVacuum used in this paper open source [27]. This updated version of BaikalVacuum and its sister (nonvacuum BSSN) code Baikal will be incorporated into a future Einstein Toolkit release.  \nWe next describe each improvement in detail.", 'A. Curvature-Adjusted Kreiss-Oliger dissipation (CAKO)': "Kreiss-Oliger (KO) dissipation [28], a widely adopted low-pass filter in numerical relativity (NR), reduces numerical noise. In Cartesian coordinates, the addition of the n th-order KO dissipation operator (where n is odd) to the right-hand side of each evolution equation is formulated as:  \nD n KO ⃗ f = ϵ KO ( -1) ( n +3) / 2 2 n +1 (∆ x i ) n ∂ ( n +1) i ⃗ f, (19)  \nwhere ⃗ f denotes the vector of evolved functions on a numerical grid, ϵ KO is the dimensionless KO strength factor, and ∂ ( n +1) i represents the ( n +1)st partial derivative in the x i direction. This derivative is generally implemented by a centered finite-difference operator with the smallest possible stencil. For example, first-order KO dissipation includes the second derivative, which is represented by the standard 3-point stencil: ∂ 2 x f = ( f i -1 -2 f i + f i +1 ) / (∆ x ) 2 . The (∆ x i ) n prefactor ensures that the operator converges to zero in the continuum limit; KO dissipation is purely a numerical filter acting to damp high-frequency oscillations.  \nWhen solving a PDE with a dominant finite-difference truncation error of j th order ( j even), the n = j +1-order  \nKO dissipation operator is most often chosen, with the KO operator stencil comprising j +3 grid points. In this study, we use j = 8, and apply 9th-order KO dissipation to the central point on an 11-point stencil.  \nLarger KO strength factors ϵ KO lead to stronger filtering; however, ϵ KO is CFL-limited to values of approximately one or less when typical CFL factors and numerical integration schemes are chosen for evolving Einstein's equations in the moving-puncture formalism [29]. This said, so long as the CFL-limited value is not exceeded, ϵ KO need not be the same for all evolved grid functions.  \nIn [30], we found that setting ϵ KO , gauge = 0 . 9 for the gauge quantities and ϵ KO , other = 0 . 3 for the other spacetime quantities enhances stability when evolving boosted and unboosted near-extremal-spin ( a/M = 0 . 99) BHs. Several years later, in [10], we refined this approach in the context of a mildly spinning BBH evolution. Building on the idea that the sharp fields inside and around puncture BHs are natural and should not be smoothed, we made ϵ KO , gauge a function of time and the distance from the center-of-mass r , such that it increased gradually over time from 0.3 to 0.98 in the weak-field region. Additionally, we adjusted ϵ KO , gauge to drop to zero at punctures by multiplying by the conformal factor W raised to some fractional power.  \nWhile this refinement was quite effective in the context of the BBH evolution presented in that paper, ϵ KO , gauge ( t, r, W ) would need to be modified and potentially retuned for BBH scenarios at different initial separations and potentially different AMR grid structures.  \nHere we introduce a much simpler adjustment to ϵ KO , which we call 'CAKO': Curvature-Adjusted KO dissipation. It multiplies the KO dissipation strength by the conformal factor W = e -2 ϕ , which ensures the KO strength is nearly zero near punctures, tending to ϵ KO , CA in the weak field:  \nϵ KO ( W ) = Wϵ KO , CA . (20)  \nSimilar to our past work, we set ϵ KO , CA to 0 . 99 (near the CFL-limited maximum) for gauge quantities and to 0 . 3 for all other BSSN evolved variables. But unlike our past work, this prescription has been validated using both BlackHoles@Home grid structures across more than 100 unique scenarios across BBH parameter space, and (as shown in this work) AMR grid structures within the Einstein Toolkit in the context of a GW150914-like BBH system at multiple numerical resolutions.", 'B. Coarse-grid-Adjusted H Damping (CA H D)': 'In 2002, [31] found that analytically, the following modification to the right-hand side (RHS) of the conformal factor evolution equation causes one constraint-violating mode to exponentially damp:  \n∂ t ϕ = [Eq. (6) RHS] + κ ϕ H α H -, (21)  \nwhere κ ϕ H ≤ 0 for stability.  \nAs H = 0 in the continuum limit, this modification drops to zero with increasing resolution. Further, for unit consistency κ ϕ H must have units of mass M (i.e., the same as units of time and space in G = c = 1 units).  \nWhen analyzing the principal part of the adjusted ∂ t ϕ equation in the weak-field limit, we find  \n∂ t ϕ = H -, pp = -8 κ ϕ H ∂ i ∂ i ϕ, (22)  \nwhich is a parabolic PDE. If this parabolic equation is solved explicitly with a 3D Laplacian [10, 32] in flat space, then at a numerical grid point with ∆ s = min(∆ x i ) across all three spatial dimensions i , the CFL-constrained timestep ∆ t is given by  \n-8 κ ϕ H ∆ t (∆ s ) 2 ≤ 1 6 = ⇒ ∆ t ≤ (∆ s ) 2 -48 κ ϕ H . (23)  \nThe (∆ s ) 2 on the right-hand side is undesirable. For instance, suppose a simulation is performed at the maximum stable κ ϕ H at a given resolution ∆ s and corresponding CFL-limited ∆ t ∝ ∆ s . If a higher resolution (smaller ∆ s ) is found necessary for a convergence study, κ ϕ H would need to be adjusted downward in order for the CFL condition to be satisfied. This would result in an inconsistency in the strength parameter across different resolutions, potentially invalidating the convergence study.  \nTo address this, [21] included the simulation timestep in the prefactor:  \n∂ t ϕ = [Eq. (6) RHS] + K ∆ t H + , (24)  \nwhere K is dimensionless and K ≥ 0 for stability. Numerically, the term ∆ t ∝ ∆ s decreases the strength of the parabolic operator in direct proportion to the CFL limit, making this modification stable at all choices of grid resolution provided it is stable at a given choice of K . From an analytical perspective, both the ∆ t prefactor and H drop to zero in the continuum limit. Choosing K = 0 . 02-0 . 06, they showed that this damping term reduced constraint violations by almost an order of magnitude in the context of BNS simulations.  \nMore recently, [32] explored the effect of three different coefficient prescriptions in the context of BNS evolutions. These include: Eq. (24) with K set to a constant across all AMR grids; Eq. (24) with K decreased on coarser grids according to ∆ s min / ∆ s n (where ∆ s n is the ∆ s on the n th refinement level), and the third method,  \n∂ t ϕ = [Eq . (6) RHS] + c H H + , (25)  \nwhere c H is set to a fixed dimensionful constant (independent of ∆ t ). They found the third approach preferable because it improved the convergence properties of the post-merger phase. Additionally, they attributed the decreasing effectiveness of the [21] approach at higher resolutions to its relatively poor convergence behavior. As discussed, setting the coefficient on the parabolic operator to a fixed value results in instability if the evolution is repeated with sufficiently high resolution, so this approach  \nis best applied when the highest resolution can be known a priori .  \nOur improvement to the approach inverts the second prescription of [32]. Instead of decreasing the strength K on coarser grids, we increase the strength according to the maximum allowed timestep on the refinement level, similar to the approach of [14] for accelerating hyperbolic relaxation waves across a curvilinear numerical grid. Specifically, we adjust Eq. (24) as follows:  \n∂ t ϕ = [Eq. (6) RHS] -C ∆ t n [ CFL 0 ∆ s n ∆ t n ] H -, (26)  \nwhere the n subscript refers to the local value of a given quantity on grid n , and CFL 0 is the CFL factor on numerical grid n = 0, which contains the smallest timestep. We find C = 0 . 15 results in stable evolutions across BBH parameter space explored by BlackHoles@Home , including nonzero spins and mass ratios from q = 1 to q = 6, and we generally set CFL 0 to 90% of its highest possible value as a safety margin.  \nOn the finest AMR level ( n = 0), the bracketed term reduces to 1. This approach maximally benefits standard (2:1 refined) AMR grids with a global timestep ∆ t n = ∆ t 0 = ∆ t , as it exponentially increases the damping strength with n ; in this limit the bracketed term reduces to 2 n . For the ∆ t n used in this work (Sec. IV), the bracketed term increases the damping strength by a factor of 2 5 = 32 on the coarsest grid.', 'C. Slow-Start Lapse': "Our third and final improvement focuses on rapid changes in the lapse near puncture BHs at the start of a spacetime evolution. As illustrated in Fig. 1, these changes send a sharp pulse outward at √ 2 c . To minimize the early dynamics, the initial lapse α 0 is typically set to a 'pre-collapsed' value of α 0 = W 0 , where W 0 = e -2 ϕ 0 is the conformal factor W at t = 0.  \nOur slow-start lapse (SSL) technique modifies the 1+log lapse evolution condition Eq. (9) as follows:  \n∂ t α = [Eq. (9) RHS] -W [ he -t 2 / (2 σ 2 ) ] ( α -W ) . (27)  \nSSL consists of three terms multiplied together. First, -( α -W ) exponentially damps α to the (evolved) conformal factor W . By convention, at t = 0, α is set to W , rendering this term exactly zero at t = 0. Second, the Gaussian prefactor in square brackets sets the inverse damping timescale, such that the damping timescale is short and effective early in the evolution but grows to a very large and inert value later on. Third, to ensure that the wormhole-to-trumpet [33] transition inside the horizon is not affected, the damping strength is multiplied by an overall prefactor of W , which drops to zero at the puncture and smoothly increases away from the puncture.  \nIn short, the SSL modification acts to quench rapid deviations of α from the quantity to which it is initially set, W , just outside the horizon, arresting the development of the sharp lapse feature in this region.  \nAfter many experiments across ∼ 100 unique scenarios in BBH parameter space with BlackHoles@Home , we find that a Gaussian height of h = 3 5 M and a standard deviation of σ = 20 . 0 M are particularly effective in mitigating the initial lapse pulse. With these choices, SSL damps on a timescale of e t 2 / (2 σ 2 ) /h = 5 3 M at t = 0, ∼ 10 2 M at t = 57 M , ∼ 10 5 M at t = 94 M , and ∼ 10 16 M at t = 170 M . In this way, SSL impacts only the formation of the sharp lapse feature and is rendered completely ineffective shortly thereafter.", 'IV. ALGORITHMIC APPROACH': "We present reproductions of the GW150914 gallery example [34] from the Einstein Toolkit [22, 23], incorporating the aforementioned improvements in a systematic manner. This gallery example is notable, as it is perhaps the most comprehensively documented BBH NR simulation available; all NR codes, results, and plotting scripts required to reproduce the original simulation are openly accessible on Zenodo [34].  \nThe physical parameters of the simulation are inspired by the first-ever directly detected GW event, GW150914 [35]. In particular, it uses the best-fit parameters for GW150914 [36]: a mass ratio of q = 36 29 , with dimensionless spin parameters of ⃗ χ M = { χ x M , χ y M , χ z M } = { 0 , 0 , 0 . 31 } for the more-massive BH and ⃗ χ m = { χ x m , χ y m , χ z m } = { 0 , 0 , -0 . 46 } for the lessmassive one. The initial separation is set to 10 M , which places the BHs approximately 5 orbits before commonhorizon formation.  \nAs in the gallery example, we use the Einstein Toolkit [22, 23] and its Carpet [24, 25] AMR infrastructure to perform these simulations. The numerical parameters for the evolution are largely the same: 8thorder finite differencing, 5th-order prolongation in space, 2nd-order prolongation in time, RK4 timestepping. The TwoPunctures thorn is used to construct initial data [37], AHFinderDirect [38] for apparent-horizon diagnostics, QuasiLocalMeasures [39] for isolated horizon diagnostics, and the Kranc [40]-generated WeylScal4 for ψ 4 extraction. However, there are some notable differences.  \nFirst, instead of the ML BSSN [40-42] thorn, the NRPy -generated [18] BaikalVacuum thorn, modified as needed with the improved techniques, is used for the spacetime evolution. BaikalVacuum and its sister code Baikal are BSSN evolution codes in the Einstein Toolkit : BaikalVacuum solves the BSSN equations without stressenergy source terms, and Baikal includes these source terms. The latter was used in [43] and [44] for the evolving the spacetime fields in BNS simulations.  \nSecond, the TwoPunctures grid resolution was increased from its default values of N A × N B × N ϕ = 30 × 30 × 16  \nto 48 × 48 × 20, as the default values result in initial data constraint violations dominating the early evolution. Also, the initial lapse was chosen to be W = e -2 ϕ , instead of the ' TwoPunctures-averaged ' value, which sets the lapse to α ( r ) = 1 2 ∑ 2 i =1 (1 -m i / (2 r i )) for each puncture i of mass m i at distance r from r i . The former approach has the advantage of exhibiting the same asymptotic behavior at large radius and late times: α ∞ = 1 -M/r for total mass M . Meanwhile, the latter approach tends to 1 -M/ (2 r ) in the weak field. While TwoPunctures has an option to set the initial lapse to any power of ψ = e ϕ , it chooses a value of ψ that does not include the correction term u (Eq. 30 in [37]). To set the initial value of the lapse to the evolved quantity W , as needed for SSL, we modified TwoPunctures to include this option. Our modified version may be found in this GitHub repository [27].  \nThird, the Llama [45] cubed-spheres grid in the gravitational wave zone was replaced with box-in-box Cartesian adaptive-mesh refinement (AMR) grids of comparable or higher resolution. In addition, we choose a default resolution of all AMR grids 50% higher, so that, for example, the resolution at the punctures is M/ 42 instead of the gallery example value of M/ 28. As we will find in Sec. V, even M/ 42 lies slightly outside the convergent regime for modeling this physical scenario. For completeness, our AMR grid hierarchy comprises 11 levels of overlapping grid cubes, centered on and comoving with each puncture, with a standard factor-of-2 refinement across adjacent AMR levels (i.e., there are 10 levels of refinement built upon the coarsest level). The AMR grid cubes have half-sidelengths of { 2 3 M, 4 3 M, 8 3 M, 16 3 M, 64 3 M, 128 3 M,..., 4096 3 M } . As our default choice we choose grid resolutions of { M 42 , M 21 , M 10 . 5 , ..., 24 . 38 M } with CFL factors of 1 640 { 288 , 288 , 144 , 72 , 36 , 18 , 9 , 9 , 9 , 9 , 9 } , respectively.  \nFinally, the imposed symmetry across the z -axis was removed as BaikalVacuum does not currently support symmetries.", 'V. RESULTS': 'To investigate the individual and combined effectiveness of the improvements introduced in this work, we conducted five simulations based on the GW150914-inspired scenario detailed in Sec. IV at the default resolution. These simulations are denoted as Control, CAKO, CA H D, SSL, and All3. Control serves as the no-improvement baseline, while the CAKO, CA H D, and SSL each integrate a single enhancement: curvature-adjusted Kreiss-Oliger, H -damping, and slow-start lapse, respectively. All3 combines all three enhancements. For convergence analyses, we executed the Control and All3 simulations at resolutions 1.25x, 1.5x, and 1.71x higher than the default resolution (1.0x) [46].', 'A. Effect of Improvements on Trajectories and Lapse Evolution': "FIG. 3. Trajectory comparison: Control vs. All3 . Top panels, upper curves : distance from the origin of the more massive BH d M , comparing Control (blue solid) and All3 (orange dashed) over the first three orbits. Top panels, lower curves : least-squares cubic polynomial fit P 3 ( t ) to All3 (red dotted) and Control (green dotted), minus 0.1 to improve visibility. Bottom panels : difference between data and cubic polynomial fit.  \n<!-- image -->  \nAs the three improvements modify the evolution equations in the finite-resolution limit, as well as the 1+log lapse condition, it is reasonable to be concerned that All3 may influence the trajectories of the BHs, especially at lower resolutions, possibly inducing unwanted eccentricity. To investigate this, in Fig. 3 a cubic polynomial P 3 ( t ) is fit to d M ( t )-the distance of the more massive BH from the center of mass (the origin)-over the first three orbits (0 ≤ t ≤ 600 M ).  \nd M ( t ) and P 3 ( t ) -0 . 1 are plotted in the top panels of each plot. Comparing the top and bottom plots, the P 3 ( t ) fits-and thus the trajectories-appear to converge with increased resolution. In addition, the sharp non-  \nconvergent oscillation in Control at t ≈ 20 M is eliminated in All3. Comparing the effects of the individual improvements, we find that the SSL enhancement is responsible for the removal of this sharp feature. As SSL most significantly impacts the early gauge evolution, this is not surprising.  \nRegarding possible eccentricity enhancement in All3, the residual d M ( t ) -P 3 ( t ) (bottom panels of Fig. 3) shows the opposite is the case: at both default and highest resolutions, the coordinate eccentricity over the first three orbits is slightly lowered in All3 as compared to Control. That is to say, the amplitude of the residual oscillation is reduced in All3 as compared to Control. Further analysis reveals that All3 reduces the average absolute value of the residual by 33% at low resolution, from 0.0088 in Control to 0.0059 in All3; and by 25% at the highest resolution, from 0.0059 to 0.0044. We stress that this is a measure of the smoothness of d M ( t ), not eccentricity.  \nWe next turn our attention to further analyzing the effects of SSL. As discussed in Sec. III C, SSL slows the deviation of the lapse away from α = W outside the horizon at early times. Given that we adopt the standard 'pre-collapsed' lapse choice α = W at t = 0, SSL has no impact on the evolution at t = 0.  \nFigure 4 illustrates how the evolution of the lapse deviates from the conformal factor shortly after the simulation begins, comparing α -W in Control to SSL at times t/M = 1 . 37, 5.49, and 10.97. For comparison, the first orbital period completes at t ≈ 222 M .  \nThe sharp lapse feature propagates outward and induces a great deal of constraint violation in its wake, as discussed in Sec. I. Figure 5 shows how each improvement impacts the sharp lapse feature evolution in the wavezone. As in Fig. 1, r (1 -α ) is plotted, which should tend to a constant value in the limit r →∞ .  \nThe more massive puncture is initially situated at x = 4 . 46, moving off the x -axis shortly after the evolution begins. Thus the top panel shows the evolution very close to the puncture, and lower panels depict slices at larger distances from the puncture and at later times. The conformal factor along the axis smooths significantly away from the puncture, so the oscillations in α -W at t ≥ 5 . 49 M are primarily due to the lapse. While the effects of SSL are barely noticeable at t = 1 . 37 M , they become quite apparent starting at t = 5 . 49 M ; SSL causes the oscillations in the lapse to be greatly reduced in amplitude, both at t ≈ 5 . 49 M and t ≈ 10 . 97.  \nAt t = 175 . 5 M , the wavefront in Control, CA H D, and CAKO has propagated to roughly x = 250 M , consistent with a wave propagating from near the origin at t = 0 at a speed of √ 2 c . Notice this feature partially reflects off of the x = 170 . ¯ 6 M AMR grid boundary, due to undersampling across this factor-of-two drop in resolution. Meanwhile, the delayed evolution of the lapse away from W early in the SSL and All3 evolution acts to delay the wavefront by approximately t = 26 M ; these wavefronts at t = 201 . 6 M overlap those of runs without SSL at t = 175 . 5 M .  \nFIG. 4. Effects of SSL on Early Lapse Evolution : α -W at various times during the early evolution along the x -axis, in the vicinity of the more massive BH, comparing Control (blue solid) with SSL (orange dashed).  \n<!-- image -->  \nFIG. 5. Impact of improvements on the sharp lapse feature : r (1 -α ) is plotted at t/M = 175 . 5 along x -axis for Control (thin solid), CA H D (thin dotted; overlaps Control), and CAKO (thin dashed); and at t/M = 201 . 6 for SSL (thick solid) and All3 (thick dotted; overlaps SSL). Dashed vertical lines demarcate AMR boundaries along the x -axis at x = 85 . ¯ 3, 170 . ¯ 6, and 341 . ¯ 3. The transmitted and reflected wave pulse propagation directions in Control, CA H D, and CAKO across the 170 . ¯ 6 AMR grid boundary are denoted by the horizontal arrows.  \n<!-- image -->  \nAs the CA H D and Control results are indistinguishable, we conclude that CA H D has no significant impact on the evolution of the lapse itself. On the other hand, CAKO acts as a strong low-pass filter in this weak-field region, significantly damping oscillations in the sharp lapse feature and its reflections over the standard KO  \nprescription adopted in Control and CA H D.  \nMost notably, the SSL modification to the 1+log lapse condition significantly smooths the outgoing lapse pulse, resulting in virtually no reflection across the x = 170 . ¯ 6 M grid boundary. While this grid boundary generates minimal noise in SSL and All3 (which perfectly overlap in this plot), partial reflection is observed in SSL and All3 when the wavefront crosses the x = 341 . ¯ 3 M grid boundary. After partially reflecting, All3 exhibits less noise than SSL due to the enhanced KO dissipation implemented in CAKO.", 'B. Strong-Field Region: Effect of Improvements on Constraint Violations': 'As discussed in [10] and illustrated in Fig. 1, the outward-propagating sharp lapse pulse fills the computational domain with significant constraint violations in its wake. The top plot in Fig. 6 compares how our improvements, both individually and collectively, reduce the volume-averaged L 2 norm of the Hamiltonian-constraint violation ||H|| in the strong-field region.  \nCA H Dhas the most significant influence on ||H|| , reducing it by between 2-3 orders of magnitude after t = 10 M . This is consistent with its intended effect of causing one Hamiltonian-constraint-violating mode to exponentially damp. Further, All3 reduces ||H|| in the strong field by an additional 0.2-1.2 orders of magnitude. The spike in constraint violation observed at t ≈ 930 M is coincident with common horizon formation. While the same spike is only observed in CA H D and All3, the Hamiltonianconstraint violation in these cases never exceeds Control, CAKO, or SSL.  \nMeanwhile, as shown in the bottom plot of Fig. 6, at most times SSL greatly reduces momentum-constraint violation √ ||M 2 || in the strong field at early times, decreasing it by between 0-3.5 orders of magnitude. After the merger, CAKO proves to be most effective at reducing momentum-constraint violation. Once more, the cumulative impact of all improvements (All3) surpasses that of any individual improvement.  \nNext we analyze the convergence behavior of the constraint violations. Assuming we are in the convergent regime (i.e., the region where truncation, and not roundoff or undersampling error dominates), we would expect the constraints to converge to zero between second and eighth order in our mixed-order scheme; second order corresponds to the order at which data are interpolated (prolongated) in time at AMR grid boundaries, and eighth order corresponds to the finite-differencing order for spatial derivatives. Witnessing inconsistent convergence order in such a complex system typically implies that some physical lengthscale is not being properly resolved (undersampling error).  \nGiven that constraint violations converge to zero in the continuum limit, in the truncation-error-dominated regime, we should find that the constraint violation C i at', '|| || in the range r 20 M': "Figure 7 plots the base-10 logarithm L 2 norms of Hamiltonian ( ||H|| ) and momentum ( √ ||M 2 || ) constraint violations in the strong-field region, through merger and ringdown, up to t = 1000 M . After this, the key payload from the simulation-the GW signal-has largely propagated away and no longer depends on data in this  \n<!-- image -->  \n2  \nFIG. 6. Effect of improvements: Constraint violation in strong-field region : Volume-averaged L 2 norm of ||H|| ( top plot ) and √ ||M 2 || ( bottom plot ) inside spherical region of radius 20 M , centered at the origin (center of mass), excluding two spheres of radius 1 M centered at each BH. Results from five runs are compared: Control (thick solid blue), CAKO (dash-dotted orange), CA H D (dashed green), SSL (dotted red), and All3 (thin solid purple). Volume excludes spheres of radius r = 1 M centered at each BH.  \n<!-- image -->  \nresolution ∆ x i and convergence order n should approximately obey the equation:  \nC i = A (∆ x i ) n , (28)  \nwhere A is related to the ( n +1)st derivative of the term in C i that dominates the truncation error. So long as C i is known at at least two resolutions, both n and A can be computed.  \nregion. Data from Control and All3 at the default (1.0x) resolution, as well as 1.5x and 1.71x higher resolutions are shifted on this logarithmic scale such that the difference with the corresponding 1.25x case is minimized; from this difference an implied convergence order n can be immediately computed using Eq. 28:  \nn = log 10 C i -log 10 C 1 . 25 x log 10 ∆ x i -log 10 ∆ x 1 . 25 x . (29)  \nThe top plot of Fig. 7 indicates that through most of the evolution, Hamiltonian-constraint violation in All3 remains about 100 times lower than Control at 1.25x resolution. More strikingly, the implied convergence orders at higher resolutions are consistently higher in All3 (4.6-6.4) compared to Control (3.5-4.2), indicating that this gap only widens at higher resolutions.  \nAll3 is even more effective at reducing momentumconstraint violations in the strong-field (bottom plot of Fig. 7), dropping them by about 3-4 orders of magnitude as compared to Control at 1.25x resolution. Unlike the implied convergence order from Hamiltonian-constraint violations, the momentum constraint convergence order is nearly doubled at higher resolutions in All3 (5.9-7.1) compared to Control (3.0-3.2), reflecting that at higher resolutions, the 3-4 decade gap between Control and All3 also widens. Finally, 1-2 order-of-magnitude bumps in √ ||M 2 || appear periodically as reflections from distant AMR boundaries enter this region. These reflections strongly impact All3's convergence around t/M = 200 and 400, but not at other times.  \nAnalyzing implied convergence orders in Fig. 7, two results stand out: All3's Hamiltonian-constraint violation exhibits n = 9 . 6 (1.0x vs 1.25x) and Control's momentumconstraint violation exhibits n = 1 . 6 (1.0x vs 1.25x). These implied convergence orders lie outside the expected range for our mixed-order scheme of n = 2 to n = 8, indicating that our chosen default resolution may be outside the convergent regime.  \nUnusual implied convergence orders warrant further investigation. In Fig. 8, we examine Hamiltonian-constraint violation in the region where the binary orbits, along the x -axis at various times when the BHs are close to the y -axis. Control is characterized by sharp, wavelike features about H = 0, while All3 exhibits a purely negative (at the lowest resolution) or mostly positive (at higher resolutions) systematic bias, albeit with constraint violation amplitudes roughly 1 / 100th that of Control (consistent with Fig. 7). For scale, All3 data would appear nearly indistinguishable from the horizontal dashed lines marking H = 0 in the top plot.  \nWhile the amplitude of constraint violation indeed drops with increased resolution in both Control and All3, there is a clear zero crossing in All3 when comparing lowest and 1.25x higher resolution cases, consistent with the lowest resolution not yet being in the convergent regime. Further, the 1.25x resolution does not exhibit a clear bias in All3, possibly implying that it lies very close to the  \nFIG. 7. Convergence of constraint violation in strongfield region : Convergence to zero of volume-averaged L 2 norms of Hamiltonian ( top plot ) and momentum ( bottom plot ) with increased resolution in a spherical region of radius 20 M , centered at the origin (center of mass), excluding two spheres of radius 1 M centered at each BH. Control and All3 results are compared. The lowest resolution, as well as 1.5x and 1.71x higher resolution results have been shifted up or down to maximally overlap the 1.25x higher resolution results. Shifts and implied convergence orders ( n ) are reported in the legend for each case.  \n<!-- image -->  \nboundary between the convergent and non-convergent regimes  \nAs indicated by Eq. (15) and surrounding discussion, H violations act as an effective energy density ρ eff ; thus, the consistent ρ eff < 0 in 1.0x All3 implies that Hamiltonianconstraint violation acts to continually push the BHs away from one another. If this were the dominant error, we would expect a delayed merger. Following the same reasoning, the consistent ρ eff > 0 in 1.5x and 1.71x resolutions of All3 imply an early merger, and the 1.25x resolution somewhere in-between. As will be discussed in our analysis of the extracted waveforms (Sec. V C), this is exactly what we observe. Such a clean analysis  \nis impossible in Control due to the extremely high level of noise, and implies that clean waveform convergence in Control might depend on near-exact cancellation of H waves.  \nFigure 9 shows that, after the early settling of M irr , the irreducible mass remains constant to within a couple of parts in 10 5 across all cases. Comparing 1.0x with  \n<!-- image -->  \n-  \nFIG. 8. Hamiltonian-constraint violation in the binary region during inspiral : Convergence of H (thin solid curves) in Control with increased resolution ( top plot ) and All3 ( bottom plot ) along the positive x -axis at times when the BHs are very close to the y -axis, starting at approximately 0 . 75 orbits: t/M ≈ 168, 272, 371, 464, 549, 628, 701, 765, and 828. A thick, gray, dashed line denotes y = 0. Note that both plots zoom out the top panel to better fit the data. The BHs are initially placed at x = 4 . 46 and x = 5 . 54, and the origin is the center of mass. Merger occurs at roughly t/M = 930.  \nWe conclude that examining H in the binary region serves some use in identifying whether BBH evolutions are in the convergent regime, at least in the All3 cases. Regarding Control we found a low implied convergence order of ||M 2 || when comparing 1.0x and 1.25x resolutions. Further analysis of the irreducible mass (Fig. 9) of the less massive, faster-spinning BH lends additional insight.  \n1.25x and 1.5x resolutions, Control exhibits a negative implied convergence order beyond t ≈ 150 M , providing further evidence that the 1.0x resolution data in the strong-field are not in the convergent regime. Measuring n directly for Control is quite challenging due to numerical noise, though given the stochastic behavior observed in H near the punctures (Fig. 8), this is unsurprising. That said, comparing 1.25x, 1.5x, and 1.71x resolution data yields an implied convergence order of roughly n = 1 at t = 400 M , but the 1.25x and 1.5x data completely overlap at t ≈ 600 M (implying n = 0). Similar convergence behavior was observed in this quantity in [10].  \nMeanwhile, All3 exhibits far less noise and cleaner convergence in this diagnostic, with an implied convergence order of n ≈ 4 . 7 in the range 200 < t/M < 600 when analyzing the 1.0x, 1.25x, and 1.5x data, and n ≈ 2 . 0 when analyzing the 1.25x, 1.5x, and 1.71x data. Although the All3 data have implied convergence orders consistent with the range expected from our mixed-order scheme (2 ≤ n ≤ 8), Fig. 8 has established that constraint violations in this region in the 1.0x All3 simulation are not in the convergent regime.  \nFIG. 9. Convergence behavior of M irr , for smaller, faster-spinning BH : Comparison of irreducible mass for Control ( top panel ) and All3 ( bottom panel ) cases at all four resolutions. Thicker lines imply coarser resolution. Data beyond t = 800 M are removed, as this is close to the time of common horizon formation; merger occurs at t ≈ 910 M .  \n<!-- image -->", 'C. Weak-Field Region': "While numerical noise in Control significantly influences the constraints and M irr in the strong field, it is most prominent in the weak field region where GW contenttypically the core payload of BBH simulations-is extracted. As our improvements act to mitigate this noise, they are particularly effective in the weak-field region. Figure 10 demonstrates that throughout much of the in-  \nFigure 12 illustrates our procedure for quantifying noise in the dominant, ℓ = m = 2, spin-weight s = -2 mode of ψ 4 . In the top panel, we plot the amplitude of this mode extracted from the simulation at a radius R ext = 100 M , for both Control and All3 at 1.25x higher resolution than the lowest resolution. The vertical dashed lines denote the region where a quartic polynomial was fit through the amplitude data. The middle panel plots the absolute difference between the data and the quartic polynomial fit (i.e., the 'noise'), and the bottom panel the noise integrated over time. Averaged over the time interval of the fit, our combined improvements reduce numerical noise in the ℓ = m = 2 mode of ψ 4 by a factor of 4.0 at 1.25x resolution. Applying the same technique to other resolutions, we find noise reduction factors of 6.0 (1.0x), 4.3 (1.25x), and 3.1 (1.71x). This underscores the fact that this noise converges away extremely slowly at higher resolutions.  \n<!-- image -->  \nFIG. 10. Effect of improvements: Constraint violation in gravitational wavezone : Volume-averaged L 2 norm of ||H|| ( left ) and √ ||M 2 || ( right ) within a thick spherical shell of 20 ≤ r/M ≤ 500, centered at the origin (center of mass). Results from five runs are compared: Control (thick solid blue), CAKO (dash-dotted orange), CA H D (dashed green), SSL (dotted red), and All3 (thin solid purple).  \n<!-- image -->  \nspiral ( t < ∼ 1000 M ), Hamiltonian-constraint violations in the GW extraction region (the 'wavezone') are reduced by roughly four to six orders of magnitude in All3, largely a consequence of SSL. This contrasts with the strong-field region (Fig. 6), where 2-3 orders of magnitude reductions are observed. Meanwhile, during the inspiral, All3 reduces momentum-constraint violations by between 2-3 orders of magnitude, similar to the strong-field region.  \nAt the merger, a spike in Hamiltonian- and momentumconstraint violations propagates from the strong-field region into the wavezone. In Control, this spike is observed only in momentum-constraint violations, likely because noise levels in this region are far higher. Post-merger, Hamiltonian-constraint violation reductions in All3 equilibrate to values roughly 3-4 orders of magnitude below Control, with near-equal contributions from SSL and  \nCA H D.  \nCareful analysis of momentum-constraint violations in this region after merger shows that CA H D and SSL overlap Control, making CAKO the dominant contributor to reducing momentum-constraint violations. While All3 always falls below Control, they are quite comparable as constraint violations associated with common horizon formation propagate outward. Later, the gap between Control and All3 slowly widens to about one order of magnitude by the end of the simulation.  \nTurning now to constraint convergence behavior in the wavezone at higher resolutions, Fig. 11 plots the constraint violation at all four resolutions, using the 1.25x resolution as the basis of comparison. In Control, neither Hamiltonian nor momentum-constraint violations converge to zero in the expected range of 2 ≤ n ≤ 8, while, remarkably, All3 exhibits implied convergence orders of 3.5-3.8 in the Hamiltonian-constraint violation. However, All3 falls outside the expected convergence-order range in momentum-constraint violation, and in fact is slightly anti-convergent. While this is rather troubling for both Control and All3, upon closer analysis, we find that the level of violation in All3 remains lower than Control at all resolutions, just as in Fig. 7. We also find that if the integral is extended to the outer boundary as in [10], we find a consistent result: both Control and All3 exhibit convergence orders in the expected range 2 ≤ n ≤ 8, implying convergent behavior outside of the wavezone.  \nConsidering the GWs themselves, while GW detectors measure variation in GW strain, NR simulations typically compute the Weyl scalar ψ 4 . As ψ 4 is the second time derivative of the strain in the weak-field limit-dependent on second spatial derivatives of the metric-it amplifies intrinsic numerical noise and thus provides an excellent noise diagnostic.  \nSuch significant noise reductions in the dominant ψ 4 mode are quite important, as GWs are usually the core payload of NR BBH simulations. Looking ahead, the improved sensitivities of current and future GW observatories will require significantly improved accuracy of  \nFIG. 11. Convergence of constraint violation in gravitational wavezone : Convergence to zero of volume-averaged L 2 norms of Hamiltonian ( top plot ) and momentum ( bottom plot ) with increased resolution in a thick spherical shell of inner radius 20 M and outer radius 500 M , centered at the origin (center of mass). Lowest resolution as well as 1.5x and 1.71x higher resolution results have been shifted up or down to maximally overlap the 1.25x higher resolution results. Shifts and implied convergence orders ( n ) are reported for each case.  \n<!-- image -->  \nNR GW predictions [1, 2], including higher-order modes. In anticipation of this, Fig. 13 shows All3's impact on reducing noise in subdominant, even ℓ = m modes at the highest (1.71x) resolution. The noise at ℓ = m = 4 is greater in both All3 and Control than at ℓ = m = 2, yet noise is still significantly reduced in All3. Turning to the ℓ = m = 6 subdominant mode, we find a truly exciting result. This mode is completely obscured by noise through the entire inspiral in Control, yet is mostly visible in All3, enabling the strain prediction to be extracted!  \nStill, there are limits to All3's effectiveness: the ℓ = m = 8 mode cannot be extracted from either All3 or Control. This underscores the need for further research into improved techniques and numerical grid structures  \nFIG. 13. Subdominant modes : Amplitude of spin-weight s = -2 subdominant even ℓ = m modes of ψ 4 up to 8, versus retarded time, at highest resolution. Top panel : ℓ = m = 4 mode, Middle panel : ℓ = m = 6 mode, Bottom panel : ℓ = m = 8 mode.  \n<!-- image -->  \n-  \nFIG. 12. GW noise analysis: dominant mode . Top panel : Amplitude of dominant ( ℓ = m = 2) spin-weight s = -2 mode of ψ 4 versus retarded time. Vertical dashed lines: region where a quartic polynomial is fit to amplitude data. Middle panel : Absolute difference between data and fit. Bottom panel : Cumulative absolute difference between data and fit.  \nsuch as those that Llama or BlackHoles@Home provide.  \n<!-- image -->  \n-  \nWe next turn our attention to the convergence behavior of ψ 4 , splitting our analysis into the phase and amplitude of the dominant ℓ = m = 2 mode. The phase in Control (top two panels of Fig. 14) is too noisy to measure n until a retarded time of about 780 M , when monotonic convergence behavior starts to emerge. However, although the implied convergence order of 2.3 between 1.0x, 1.25x,  \nand 1.5x resolutions falls within the expected range of 2 ≤ n ≤ 8, when comparing differences in 1.25x, 1.5x, and 1.71x resolutions, n is roughly 12 before merger and 23 after merger, reflecting serious problems with convergence. We attribute this unusual behavior to the high-amplitude, stochastic oscillations in H between the BHs illustrated in Fig. 8. As discussed in Sec. II C, these oscillations equate to spurious positive and negative effective energy densities, which, when they do not cancel, will impose an effective attraction or repulsion between the BHs, respectively. We fear that non-cancellations may directly impact the GW phase in potentially unpredictable ways.  \nReturning to Fig. 8, H is roughly 100x lower between the BHs in All3 as compared to Control. In spite of this, the persistent negative effective energy density at 1.0x resolution likely causes a significantly delayed merger in All3 (middle panel, Fig. 14). While this may seem troubling for All3, it is even more troubling for Control, as it implies that the stochastic behavior in H between the BHs would need to cancel to better than one part in 100 over the inspiral in order to exhibit similar or lower phase errors. Given the small number of waves observed in the binary region, such cancellation cannot be guaranteed, and is likely the cause of the inconsistent convergence in Control.  \nConsidering All3 at 1.25x and higher resolutions, Fig. 8 indicates that All3 exhibits slight negative oscillation about H = 0 in 1.25x, but a consistently positive bias at 1.5x and 1.71x. Given inconsistent convergence observed in Control due to oscillations, as well as the clear negative bias of H in 1.0x All3 results, caution is clearly warranted when considering All3 1.25x ψ 4 phase results. The fourth panel from the top of Fig. 14 indicates that while the difference between 1.25x and 1.5x phase results is smaller in All3 than Control (second panel from the top), the difference flips sign close to merger, as compared to all other resolution pairs (All3 or Control).  \nThis likely implies that 1.25x resolution in All3 and Control are slightly outside the convergent regime. While this analysis of H in the binary region (Fig. 8) to explain non-convergent behavior has been useful, it remains to be seen whether it will be predictive of convergence in other BBH scenarios.  \nNext, we assess whether Control and All3 waveforms converge to the same result. To this end, in the bottom panel of Fig. 14, we fix the resolution and simply plot the phase differences between Control and All3. Notice that phase differences between Control and All3 drop to tiny fractions of a radian with increasing resolution, indeed consistent with them converging to the same result.  \nFor completeness, we present the same analysis for the rescaled amplitude of the extracted GW signal: R ext | ψ 4 | , for the dominant mode in Fig. 15. Perhaps unsurprisingly, the conclusions drawn from our phase analysis hold for the amplitude as well, though the significantly delayed merger in the 1.0x resolution All3 case is much more evident when considering the amplitude.  \nWe conclude that improvements do not appear to neg-  \nFIG. 14. GW phase convergence analysis: ℓ = m = 2. Top panel : accumulated phase | ϕ | of dominant ( ℓ = m = 2) mode of ψ 4 at all four resolutions in Control. Second panel from top : Accumulated phase difference across neighboring resolutions in Control. Third panel from top : Same as top panel, but for All3. Fourth panel from top : Same as second panel, but for All3. Bottom panel : Accumulated phase difference between Control and All3 at the three highest resolutions. All data were extracted from the simulation at extraction radius R ext = 100 M .  \n<!-- image -->  \n-  \natively impact the dominant mode of the GW prediction; both amplitude and phase appear to converge to the same values as Control. That said, 1.25x resolution All3 amplitude and phase results are much closer to the converged solution than in 1.25x Control, implying smaller errors at resolutions bordering the non-convergent regime. Perhaps most importantly, we find that All3's enormous reduction in Hamiltonian constraint noise in the binary region and the associated purely positive or negative biases, may enable us to predict for other BBH scenarios whether the actual merger time should be earlier or later than the given resolution. If so, Hamiltonian constraint information in the binary region, combined with All3, may be used to improve GW predictions from moving-puncture BBH merger simulations in the future.", 'VI. CONCLUSIONS AND FUTURE WORK': 'We have introduced three improvements to the movingpuncture technique, aimed at improving the accuracy of BBH evolutions. Our results show that the proposed improvements effectively reduce numerical errors and noise in different parts of the computational domain, from the strong-field area near the black holes to the weak-field  \nFIG. 15. GW amplitude convergence analysis . Same as Fig. 14, except instead of phase, we consider the rescaled amplitude R ext | ψ 4 | of the dominant ( ℓ = m = 2) mode of ψ 4 .  \n<!-- image -->  \n-  \nregion where GW predictions are extracted.  \nWhen combining all three improvements (All3), performance is significantly better than any single enhancement, cutting down Hamiltonian (momentum) constraint violations in the strong-field area by roughly two (three) orders of magnitude and in the GW-extraction zone by five (two) orders of magnitude. Further, the convergence of constraint violations improves near and inside the convergent regime, suggesting that these improvements become even more pronounced at higher resolutions.  \nThe improvements also slightly reduce coordinate eccentricity and lower numerical noise in the dominant ℓ = m = 2 mode of ψ 4 by an average factor of 4.3 across different resolutions. Additionally, the fidelity of higherorder modes increases significantly, making the ℓ = m = 6 mode, which was previously hidden by noise, detectable. This is particularly encouraging given the demands of  \n- [1] M. Purrer and C.-J. Haster, Phys. Rev. Res. 2 , 023151 (2020), arXiv:1912.10055 [gr-qc].\n- [2] D. Ferguson, K. Jani, P. Laguna, and D. Shoemaker, Phys. Rev. D 104 , 044037 (2021), arXiv:2006.04272 [gr-qc].\n- [3] M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlochower, Phys. Rev. Lett. 96 , 111101 (2006), arXiv:grqc/0511048.\n- [4] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and J. van Meter, Phys. Rev. Lett. 96 , 111102 (2006), arXiv:grqc/0511103.\n- [5] T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D 59 , 024007 (1998), arXiv:gr-qc/9810065.  \nnext-generation GW observatories.  \nLooking forward, we plan to explore the impacts of these improvements on BNS simulations and compare their efficacy to Z4-based formalisms. However, there are limitations to the wider applicability of these improvements: As the development of CA H D was built upon characteristic analysis of BSSN, it is unlikely to be useful in Z4-based formalisms, which exhibit a different characteristic structure and include damping terms already. Further, the initial gauge pulse is likely significantly less in BNS evolutions, meaning that SSL may not be particularly helpful in this case.  \nMoving beyond Cartesian AMR-based grids, it will be interesting to learn the impacts of these improvements on mixed AMR/cubed-spheres grid structures like Llama [45], as they may act to eliminate reflections and have already proven quite effective in BlackHoles@Home . Meanwhile, work will continue to improve upon the moving-puncture technique further, with particular emphasis on addressing spikes in constraint violation associated with merger.', 'VII. ACKNOWLEDGMENTS': "We would like to thank S. Cupp, P. Diener, R. Haas, B. J. Kelly, and L. R. Werneck for their comments on early drafts of this manuscript. We would also like to thank J. G. Baker, W. K. Black, B. Brugmann, W. East, E. Hirschmann, C. O. Lousto, D. Neilsen, V. Paschalidis, and Y. Zlochower for helpful discussions in the preparation of this manuscript. Z.B.E. gratefully acknowledges support from NSF awards PHY-2110352, OAC-2227105, AST2227080, as well as NASA awards ISFM-80NSSC18K0538 and TCAN-80NSSC18K1488. This research made use of the resources of the High Performance Computing Center at Idaho National Laboratory, which is supported by the Office of Nuclear Energy of the U.S. Department of Energy and the Nuclear Science User Facilities under Contract No. DE-AC07-05ID14517. In addition, it made use of the Falcon [47] supercomputer, operated by the Idaho C3+3 Collaboration.  \n- [6] M. Shibata and T. Nakamura, Phys. Rev. D 52 , 5428 (1995).\n- [7] D. Alic, C. Bona-Casas, C. Bona, L. Rezzolla, and C. Palenzuela, Phys. Rev. D 85 , 064040 (2012), arXiv:1106.2254 [gr-qc].\n- [8] S. Bernuzzi and D. Hilditch, Phys. Rev. D 81 , 084003 (2010), arXiv:0912.2920 [gr-qc].\n- [9] M. Alcubierre, B. Bruegmann, P. Diener, M. Koppitz, D. Pollney, E. Seidel, and R. Takahashi, Phys. Rev. D 67 , 084023 (2003), arXiv:gr-qc/0206072.\n- [10] Z. B. Etienne, J. G. Baker, V. Paschalidis, B. J. Kelly, and S. L. Shapiro, Phys. Rev. D 90 , 064032 (2014), arXiv:1404.6523 [astro-ph.HE].\n- [11] J. R. van Meter, J. G. Baker, M. Koppitz, and D.-I. Choi, Phys. Rev. D 73 , 124011 (2006), arXiv:gr-qc/0605030.\n- [12] In virtually all AMR implementations, and even for those that offer refinement factors other than two, two is generally used without exception.\n- [13] P. Amaro-Seoane et al. (LISA), arXiv preprints (2017), 1702.00786 [astro-ph.IM].\n- [14] T. Assumpcao, L. R. Werneck, T. P. Jacques, and Z. B. Etienne, Phys. Rev. D 105 , 104037 (2022), arXiv:2111.02424 [gr-qc].\n- [15] D. Hilditch, S. Bernuzzi, M. Thierfelder, Z. Cao, W. Tichy, and B. Bruegmann, Phys. Rev. D 88 , 084057 (2013), arXiv:1212.2901 [gr-qc].\n- [16] J. D. Brown, Phys. Rev. D 79 , 104029 (2009), arXiv:0902.3652 [gr-qc].\n- [17] T. W. Baumgarte, P. J. Montero, I. Cordero-Carrion, and E. Muller, prd 87 , 044026:1 (2013), arXiv:1211.6632 [gr-qc].\n- [18] I. Ruchlin, Z. B. Etienne, and T. W. Baumgarte, Phys. Rev. D 97 , 064036 (2018), arXiv:1712.07658 [gr-qc].\n- [19] W. Tichy and P. Marronetti, Phys. Rev. D 76 , 061502 (2007), arXiv:gr-qc/0703075.\n- [20] P. Marronetti, W. Tichy, B. Bruegmann, J. Gonzalez, and U. Sperhake, Phys. Rev. D 77 , 064010 (2008), arXiv:0709.2160 [gr-qc].\n- [21] M. D. Duez, P. Marronetti, S. L. Shapiro, and T. W. Baumgarte, Phys. Rev. D 67 , 024004 (2003), arXiv:grqc/0209102.\n- [22] F. Loffler, J. Faber, E. Bentivegna, T. Bode, P. Diener, R. Haas, I. Hinder, B. C. Mundim, C. D. Ott, E. Schnetter, G. Allen, M. Campanelli, and P. Laguna, Class. Quantum Grav. 29 , 115001:1 (2011), arXiv:1111.3344 [gr-qc].\n- [23] 'Einstein Toolkit Consortium Homepage,' (2024), http: //einsteintoolkit.org/ .\n- [24] E. Schnetter, S. H. Hawley, and I. Hawke, Class. Quantum Grav. 21 , 1465 (2004), gr-qc/0310042.\n- [25] 'Carpet: Adaptive mesh refinement for the Cactus Framework,' (2024), http://www.carpetcode.org/ .\n- [26] 'NRPy webpage,' (2024), https://github.com/nrpy/ nrpy .\n- [27] Z. Etienne, 'Updated BaikalVacuum with all improvements implemented,' (2024), available at https:// github.com/zachetienne/baikalvacimproved .\n- [28] H.-O. Kreiss and J. Oliger, Methods for the approximate solution of time dependent problems , GARP Publications  \n- Series No. 10 (World Meteorological Organization and International Council of Scientific Unions, Geneva, 1973).\n- [29] G. Bozzola and V. Paschalidis, Phys. Rev. D 104 , 044004 (2021), arXiv:2104.06978 [gr-qc].\n- [30] Y. T. Liu, Z. B. Etienne, and S. L. Shapiro, Phys. Rev. D 80 , 121503 (2009), arXiv:1001.4077 [gr-qc].\n- [31] G. Yoneda and H.-a. Shinkai, Phys. Rev. D 66 , 124003 (2002), arXiv:gr-qc/0204002.\n- [32] C. A. Raithel and V. Paschalidis, Phys. Rev. D 106 , 023015 (2022), arXiv:2204.00698 [gr-qc].\n- [33] J. D. Brown, Phys. Rev. D 80 , 084042 (2009), arXiv:0908.3814 [gr-qc].\n- [34] B. Wardell, I. Hinder, and E. Bentivegna, 'Simulation of GW150914 binary black hole merger using the Einstein Toolkit,' (2016).\n- [35] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 116 , 061102 (2016), arXiv:1602.03837 [gr-qc].\n- [36] B. P. Abbott et al. (LIGO Scientific, Virgo), Phys. Rev. Lett. 116 , 241102 (2016), arXiv:1602.03840 [gr-qc].\n- [37] M. Ansorg, B. Brugmann, and W. Tichy, Phys. Rev. D 70 , 064011 (2004), gr-qc/0404056.\n- [38] J. Thornburg, Classical and Quantum Gravity 21 , 743 (2004), gr-qc/0306056.\n- [39] O. Dreyer, B. Krishnan, D. Shoemaker, and E. Schnetter, Phys. Rev. D 67 , 024018 (2003), arXiv:gr-qc/0206008.\n- [40] S. Husa, I. Hinder, and C. Lechner, Comput. Phys. Commun. 174 , 983 (2006), arXiv:gr-qc/0404023.\n- [41] J. D. Brown, P. Diener, O. Sarbach, E. Schnetter, and M. Tiglio, Phys. Rev. D 79 , 044023 (2009), arXiv:0809.3533 [gr-qc].\n- [42] McLachlan, 'McLachlan, a public BSSN code,' (2008).\n- [43] F. G. L. Armengol et al. , Phys. Rev. D 106 , 083015 (2022), arXiv:2112.09817 [astro-ph.HE].\n- [44] L. R. Werneck et al. , Phys. Rev. D 107 , 044037 (2023), arXiv:2208.14487 [gr-qc].\n- [45] D. Pollney, C. Reisswig, E. Schnetter, N. Dorband, and P. Diener, Phys. Rev. D 83 , 044045 (2011), arXiv:0910.3803 [gr-qc].\n- [46] AMR grid structures and resolutions were chosen such that each resolution maintains consistent physical boundaries for AMR grids in our moving-box AMR algorithm.\n- [47] Idaho C3+3 Collaboration, 'Falcon: High performance supercomputer,' (2022)."}