Sparse Spectral-Element Methods for the Helically Reduced Einstein Equations

We describe ongoing work towards construction—via multidomain, modal, spectral methods—of helically symmetric spacetimes representing binary neutron stars. Adopting “particle” models, we focus here on solution of the helically reduced Einstein equations. These models allow us to remove the complication of a stellar fluid. However, the choice of inner boundary conditions for these models is also a stubborn issue. We examine this issue and its effect on the “harmonic gauge”.


Introduction
To model the inspiral and merger of binary objects (blackholes or neutron stars), many researchers have been solving the Einstein equations numerically. Such simulation involves both the construction of gravitational initial data at time t 0 and its subsequent evolution to a final time t F t 0 . Interpretation of experimental detections of gravitational waves relies on numerical simulation. Moreover, detection of weak signals is facilitated by statistical techniques alongside "template banks" of numerically generated signals. We consider a nonstandard problem, solution of the Einstein equations reduced by helical symmetry, as described by Beetle, Bromley, Hernández, and Price (BBHP) [1,2]. Heuristically, helical reduction is a data+evolution synthesis. Although solutions to the BBHP equations are ultimately unphysical, they may approximate the early phase of inspiral and serve as reduced order models. Moreover, they may provide excellent "trial data" (the starting point for the construction of t = t 0 initial data). Finally, they would address bewitching mathematical issues concerning exact helical symmetry in general relativity.
We consider a spectral element approach [3][4][5][6][7][8] for solving the BBHP equations. Although the equations involve a mixed-typed operator L, we solve them via relaxation using a Broyden-Krylov approach. The computational domain which surrounds the compact objects is split into 11+ subdomains (blocks, spherical shells, and cylindrical shells with classical spectral expansions thereon). To rapidly solve the linear systems arising in our scheme, we have developed sparse modal methods based on the application of spectral integration matrices, extending ideas originally described in the 1990s by Coutsias, Hagstrom, Hesthaven, and Torres.
We use preconditioned GMRES to solve these systems, with standard domain decomposition methods. In addition, we have developed fast methods for inversion of subdomain approximations of L, either via modal-based preconditioning or direct schemes.

Background
This section first describes the helically reduced wave equation (HRWE), a model for the helically reduced Einstein equations. Our HRWE description fixes ideas.

Helically Reduced Wave Equation
The wave equation is ψ = 0, where = −∂ 2 t + Δ x . Assume that ψ rotates rigidly with rate Ω. With the z-axis as the rotation axis, ψ then depends on time t only through ϕ = φ − Ωt, where φ is the azimuthal angle. Via the We adopt a "2-center domain" D, a 3d ball with two smaller 3d balls excised from it; see Fig. 1. Its boundary ∂D = ∂S − I ∪ ∂S − II ∪ ∂S + out is the union of two inner spheres (the −'s) and one outer sphere (the +). We consider the mixed-type problem For simplicity here, we have put on ∂S + out a simple Sommerfeld condition; in practice, f + is a nonlocal function of ψ enforcing an exact Dirichlet-to-Neumann map [4].

Helically Reduced Einstein Equations
We consider the vacuum Einstein equations in Landau-Lifshitz form. Write the (densitized contravariant) metric tensor as g μν = η μν − h μν , where η μν = diag(−1, 1, 1, 1) is the flat metric and h μν is the metric perturbation. Assume the harmonic gauge condition ∂ ν h μν = 0. Then the vacuum Einstein equations can be expressed as where S μνκβ τ φγ α (g) depends on g μν and its inverse g μν (but not on derivatives of either). Einsteins equations are then a constrained system of 10 nonlinear wave equations.
The BBHP reduction of the Einstein equations is similar to the one outlined for the wave equation. Technically, it assumes the existence of a Killing vector field, but we give a brief and heuristic description of their approach. The challenge is that the perturbations h μν themselves are not "helical scalars"; helical reduction is therefore not tantamount to the replacement h μν → Lh μν . However, BBHP have introduced helical scalars ψ A through which the reduction can be carried out. These are These 4 real and 3 complex quantities contain the 10 degrees of freedom in the h μν . We express this transformation as ψ A = e iμ(A)Ωt M A μν h μν , where A runs over the (tensor-spherical-harmonic) labels, and μ(A) is 0,1, or 2. Contraction of M A μν on (3), with subsequent helical reduction based on the action of on ψ A , yields Here L is the operator appearing in the HRWE. Similar to the boundary conditions appearing in (1), the boundary conditions we adopt for (5) are Again for simplicity, here we have a Sommerfeld condition on ∂S + out , but in practice use a nonlocal outgoing condition based on an exact Dirichlet-to-Neumann map.
Price has written down the analog of (2) for linearized gravity, i.e. (5) when the right-hand side of the equation is set to zero. This solution may be viewed as an exact solution to h μν = −16πT μν , where the stress energy tensor T μν corresponds to a massive point particle in an eternal circular orbit (as discussed, such a point source is excised from our domain D). Price's solution is analogous to the electromagnetic solution given by G. A. Schott, and it is given by the leading 1/R terms in the appendix expressions (12). Unfortunately, ∂ ν h μν = 0 for this solution.

Sparse Spectral Element Methods
This section summarizes our numerical methods. It is necessarily impressionistic, as even an incomplete presentation of details would take too much space.

Overview
We split D into subdomains, here with the minimal configuration of 11 subdomains: blocks B 1,3,5 ; cylinders C 1,2,3,4,5 ; an inner shell S I around "particle" I ; an inner shell S II around "particle" II ; and an outer shell S out . This corresponds to a "binary blackhole" (BBH) domain with two excised inner balls. For a "binary neutron star" (BNS) domain, we further split both S I and S II into two overlapping concentric spheres (a stellar surface then resides in each overlap [8]), and fill in the excised regions with two extra blocks B 2,4 . Figure 1 depicts a BNS domain. Pioneered by Pfeiffer et al. [9], such decompositions are used in the EllipticSolver of SpEC [10].
The unknowns in our approach are the modal expansion coefficients associated with subdomain expansions in terms of classical (Chebyshev, Fourier, and spherical harmonic) basis functions. As described below, we make extensive use of integration matrices to achieve sparse representations of the relevant operators. Before presenting details, we first address how we handle the nonlinearities in (5). Let ψ A (the vector of unknowns) be a concatenation of the modal coefficients from all 11 subdomains. Then, as sketched below, upon approximation (5) becomes where BL approximates the operator on the left-hand side of (5), with B representing the action of "integration preconditioning" (see below) on each subdomain. For linearized gravity, with the right-hand side of (5) set to zero, the vector g A is zero, save for select entries related to the inner Dirichlet values For the full Einstein equations g A depends on ψ A , and its evaluation relies on spectral analysis/synthesis (forward/backward transform) and numerical differentiation on each subdomain. We then view as a fixed-point equation, accelerating its convergence with the Broyden algorithm. This approach relies on approximation and inversion of the operator appearing on the left-hand side of (5). Reference [4] is a detailed account of the case μ(A) = 0, i.e. the HRWE. For μ(A) = 1 or 2, the operator mixes the U A and V A in ψ A = U A + iV A . We have not yet described our treatment of this scenario, but note that it relies on Schur-complement techniques. Here we describe only the μ(A) = 0 case.

Integration Preconditioning and Other Key Aspects
We use "integration preconditioning" [11] in order to achieve sparse linear systems. Especially for the cylinders and inner shells, the details are formidable. We convey the basic ideas with the Laplacian Δ, rather than L, on a cube (suppressing tildes on the co-rotating coordinates). Let u(x, y, z) ≈

. An approximation of the Poisson equation is
with u the vector of u ij k with appropriate ordering, and D 2 representing double differentiation in the modal Chebyshev basis. Let B 2 [2] represent double integration in the modal basis, where the [2] indicates that the first two rows have all zero entries. To (9) we apply the "preconditioner" [2] , thereby reaching [2] ) u = B g. (10) The system (10) is sparse, and the number of empty rows (all 0's) equals the number of tau-conditions to be enforced, auxiliary equations enforcing, say, u| ∂C = f . "Integration preconditioning" of (9) results in the sparse system (10), but the issue of condition number is subtle. Indeed, passage from (9)

Current Complexity Estimates
We aim to solve the linear systems (say approximating the HRWE or the Helmholtz equation posed on the 2-center domain D) arising in our problems at demonstrably sub-quadratic complexity; indeed, we believe that an order- 5 3 complexity is achievable. This is the complexity associated with matrix-vector multiplication; despite our sparse representations, the "gluing" of overlapping subdomains in our decomposition of D prevents realization of a linear-complexity matrix-vector product.
To document progress towards our goal, we summarize our current complexity estimates associated with solution of the HRWE on each subdomain type. These solves serve as part of our preconditioner for the global GMRES solution of the HRWE on D. For this discussion, let N represent the total number of modes on a given subdomain; e.g. N = (N x + 1)(N y + 1)(N z + 1) for the block considered above.
For S out let M be the matrix which represents r 2 L (the r 2 factor here is explained in [4]) and includes inserted tau-vectors to enforce boundary conditions. Ignoring tau-vectors, M is block diagonal in the spectral space of spherical harmonics indexed by ( , m). View its elements as M mk, m k , where mk is a "clumped index' and k, k are Chebyshev indices. Then, apart from tau-vectors, M mk, m k = 0, unless = and m = m . Moreover, each (N r + 1) × (N r + 1) block M mk, mk is itself sparse and banded. While these desirable structures are somewhat spoiled by the tau-conditions, through the use of the Woodbury formula and band solvers, for S out we are able to directly invert M (i.e. solve the HRWE on S out ) at O(N) cost.
For the inner shells, S I and S II , our representation M of r 2 L is only block banded. Indeed, the centers of these shells lie off the rotation axis, and r 2 L mixes spherical harmonic modes. Its spectral representation [4] is remarkably complicated, and relies on identities for spherical harmonics found in the treatise [12]. We solve the HRWE on inner shells via preconditioned GMRES, with a modal block-Jacobi preconditioner defined by inversion of the diagonal blocks M mk, mk . Apart from tau-conditions, these blocks are again sparse and banded, and therefore amenable to the fast methods alluded to in the last bullet. Construction and reuse of this block-Jacobi preconditioner therefore has O(N) cost. Moreover, we have empirically observed (see [4]) that such preconditioning yields low and essentially resolution independent iteration counts. While more analysis is needed, from a practical standpoint solution of the HRWE on an inner shell has an O(N) cost.
The situation on blocks is worse. Part of our global preconditioner for solving the HRWE on D involves inversion of the Poisson problem on blocks as an approximation to the HRWE (and inversion of the Helmholtz equation when the spin index μ(A) = 0). This works extremely well, likely due to the fact that Ω 1 and the blocks are close to the rotation center. In any case, we solve the Poisson/Helmholtz problem on a block via a direct approach [13] based on a rank-augmenting generalization of the Woodbury formula. This direct method is empirically well-conditioned and low-memory. Moreover, it has an O(N 2 ) set-up cost with a small constant, followed by an O(N 4/3 ) cost for subsequent solves. If possible, we hope to reduce the set-up cost to an O(N 5/3 ) complexity.
The situation on cylinders is worst of all, although to date we have not focused much attention on these subdomains. We solve the HRWE on cylinders (or the collection of cylinders) via GMRES, with modal-based preconditioners that empirically yield resolution-independent iteration counts. Application of the preconditioner currently involves an O(N 7/3 ) set-up cost, followed by an O(N 5/3 ) cost for subsequent solves. Here, we believe improvement is possible.

Numerical Tests
Our decomposition of D is from Table IV

Comparison with Exact Solutions
Our first test uses PriceLG, the aforementioned solution for linearized gravity due to Price. Here the solution is a superposition of two point-sources with the above masses. Each source point's contribution to the helical scalars is defined by the leading 1/R term in (12), and these expressions seed inner Dirichlet conditions on ∂S − I ∪ ∂S − II . The outer boundary conditions are nonlocal conditions which are exact for this solution. Table 1 lists errors for ψ (nn) ; errors for the other scalars are similar. These are relative L 2 -errors (against the exact solution) computed on both B 3 and S out via interpolation onto uniform reference grids. Since the problem is linear, each line of .0076e-16), with the rms calculation taken over the (relatively coarse) dual-nodal subdomain grids. The first and last components of ∂ ν h μν vanish for the exact PriceLG solution.
Our next test is SchwarzH, the Schwarzschild metric in harmonic coordinates: Here the radius r and the direction cosines ν = (sin θ cos φ, sin θ sin φ, cos θ) are chosen relative to a point (x 0 , y 0 , z 0 ) = (-0.9+1.37e-3,-1.6854e-4,2.9985 e-3) which is off-center but close to the center of ∂S − I . The mass is M = 0.05, and for this choice the horizon of the blackhole lies inside of (but is not concentric with) ∂S − I . For this test Ω = 0, and the exact solution (11) seeds inner Dirichlet boundary conditions on both ∂S − I and ∂S − II . On ∂S + out rather than radiation conditions, we adopt an inhomogeneous Neumann condition based on the exact solution. Table 2 lists Here relative L 2 errors are listed only for ψ (nn) . The lowest resolution run has zero initial iterate; afterwards the initial iterate stems from the previous solution. Respectively, tGMRES and iGMRES are the tolerance and iteration number for the GMRES solve 1.1266e-11 3.6899e-12 5.0e- 11 4 As for the PriceLG test, only errors for ψ (nn) are listed

553,149
----------4.0790e- 13 4 errors with the same meanings as before. Now ∂ ν h μν rms for each μ is comparable to the corresponding table errors; i.e. the gauge constraint converges to zero. Table 2 also lists the number iBROY of iterations performed by the Broyden solver to achieve the tolerance tBROY. Each Broyden iteration itself involves a linear solve via GMRES. Each GMRES tolerance is tBROY/10, the same as the corresponding line in Table 1.

Gauge Constraint Tests
Our next two tests explore to what extent the gauge constraint ∂ ν h μν = 0 is satisfied for the Einstein problem (5) and (6) in a binary scenario. For the first test we redo the PriceLG test, except now with the Einstein equations. The inner and outer boundary conditions are exactly as before. Table 3 again lists errors for ψ (nn) with the same meanings as before. Errors are computed against the finest-resolution numerical solution. The table also lists the L 2 -norm of the nonlinear residual. The last table line has ∂ ν h μν rms (3.2031e-03 5.1309e-03 4.5881e-03 4.6711e-03). That is, the gauge constraint does not converge to zero. Since the harmonic gauge is not satisfied, we cannot view these as solutions to the Einstein equations! Presumably, the violation of the harmonic gauge in the preceding example stems from the fact that ∂ ν h μν = 0 for the Price solution used to fix the inner boundary conditions. Beetle, Bromley, and, Hernández, and Price have given a refined set of inner boundary conditions, based on the near-field asymptotics of a moving Schwarzschild blackhole and meant to improve on the point-particle boundary conditions. The appendix lists (our understanding of) these conditions. With the hope that these refined boundary will result in lower gauge errors, we have performed the previous test with them. Convergence of the numerical solution is similar, but ∂ ν h μν rms has comparable size and still does not converge to zero.

Conclusions and Acknowledgments
Our tentative conclusion for helically symmetric BBH models is that violation of the gauge constraint stems from imperfect inner boundary conditions. We have also found a persistent gauge error in our BNS models, despite being several orders of magnitude smaller in that context. The BNS model, with stars in place of excised regions, involves no inner boundary conditions. We believe that in this context it is the outer boundary conditions which give rise to the constraint violation. Likely, at both the inner and outer boundaries some of the helical scalars need to be fixed using the gauge constraint itself (similar to "constraint preserving boundary conditions" used in evolution codes). The author is grateful for correspondence with Richard H. Price, and for assistance from UNM's Center for Advanced Research Computing.

Appendix: Beetle, Bromley, Hernández, and Price Inner Boundary Conditions
This appendix lists expansions for the helical scalars which somewhat generalize the ones in [1]. We need two expansions, one for a "particle" at (−a I , 0, 0), and one at (a II , 0, 0). Each has its own mass M I,II . The top choice of ± or ∓ refers to II and the bottom to I . For both a I and a II , we define v = aΩ, γ = (1 − v 2 ) −1/2 , R = γ λ ∓ vγρ sin(ϕ + Ωλ), K R = −ρ cos ϕ ± a cos(Ωλ) ± vγ R sin(Ωλ), and K I = ρ sin ϕ±a sin(Ωλ)∓vγ R cos(Ωλ). For λ see after (2). Assuming M/R 1, we have Worried about a possible sign discrepancy with the results in [1], we have also considered (12a-g) with all correction terms (those with R 4 in the denominator) flipped by a sign.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence and indicate if changes were made.
The images or other third party material in this chapter are included in the chapter's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.