Numerical relativity in spherical coordinates with the Einstein Toolkit
Vassilios Mewes Yosef Zlochower Manuela Campanelli Ian Ruchlin Zachariah B. Etienne Thomas W. Baumgarte




 Center for Computational Relativity and Gravitation, 
Rochester Institute of Technology, Rochester, NY 14623
 Center for Computational Relativity and Gravitation, 
Rochester Institute of Technology, Rochester, NY 14623
 Center for Computational Relativity and Gravitation, 
Rochester Institute of Technology, Rochester, NY 14623
 Department of Mathematics, West Virginia University, 
Morgantown, West Virginia 26506, USA
 Department of Mathematics, West Virginia University, 
Morgantown, West Virginia 26506, USA Center for Gravitational Waves and 
Cosmology, West Virginia University, Chestnut Ridge Research Building, Morgantown, 
West Virginia 26505, USA
 Department of Physics and Astronomy, 
Bowdoin College, Brunswick, Maine 04011, USA




Numerical relativity codes that do not make assumptions on spatial symmetries most 
commonly adopt Cartesian coordinates.  While these coordinates have many attractive 
features, spherical coordinates are much better suited to take advantage of approximate 
symmetries in a number of astrophysical objects, including single stars, black holes 
and accretion disks.  While the appearance of coordinate singularities often spoils 
numerical relativity simulations in spherical coordinates, especially in the absence 
of any symmetry assumptions, it has recently been demonstrated that these problems 
can be avoided if the coordinate singularities are handled analytically.  This is 
possible with the help of a reference-metric version of the Baumgarte-Shapiro-Shibata-Nakamura 
formulation together with a proper rescaling of tensorial quantities.  In this paper we 
report on an implementation of this formalism in the .  We adapt the infrastructure, originally designed for Cartesian coordinates, to handle spherical 
coordinates, by providing appropriate boundary conditions at both inner and outer 
boundaries.  We perform numerical simulations for a disturbed Kerr black hole, 
extract the gravitational wave signal, and demonstrate that the noise in these 
signals is orders of magnitude smaller when computed on spherical grids rather 
than Cartesian grids.  With the public release of our new thorns, our methods 
for numerical relativity in spherical coordinates will become available to the 
entire numerical relativity community.


 04.25.D-, 
  04.30.-w, 
  04.70.Bw, 
  95.30.Sf, 
  97.60.Lf 


Numerical relativity in spherical coordinates with the Einstein Toolkit
    Thomas W. Baumgarte
    December 30, 2023
=======================================================================







§ INTRODUCTION



LIGO <cit.> and Virgo's <cit.>
first direct detections of gravitational waves (GWs) 
from both binary black hole (BBH) and binary neutron star (BNS) 
mergers <cit.> open 
a new window for observations of the Universe. Moreover, the simultaneous detection 
of GWs and electromagnetic (EM) radiation from the BNS merger GW170817 
has launched the new field of EM-GW multi-messenger astronomy <cit.>.
Ever since the breakthrough simulations of BBHs in numerical 
relativity about a decade ago <cit.>, increasingly more
accurate models of BBH merger waveforms across the source
parameter space have been generated <cit.>.  
Together with approximate gravitational wave-form models (see, e.g., 
<cit.>), these numerical relativity simulations 
played a crucial role in the parameter estimation of GWs <cit.>
by LIGO-VIRGO <cit.>.

Among the missions of current and future GW detectors are tests of General 
Relativity (GR) <cit.>. 
While the remnant black hole (BH) mass and spin can be estimated from the 
inspiral phase <cit.>, measuring the quasinormal 
ringdown <cit.> 
of the remnant BH in the GW signal 
will provide an independent measurement of its mass and spin <cit.>, 
as well as tests of the no-hair theorem and GR <cit.>. 
Accurate modeling of the ringdown of a highly distorted remnant Kerr BH after merger is 
only possible using numerical relativity simulations of BBH coalescence through merger.
 
Arguably, one of the most widely used evolution schemes for these type of simulations
is the Baumgarte-Shapiro-Shibata-Nakamura (BSSN) 
formulation <cit.>.[Sometimes referred to as BSSNOK, 
because it is based on the strategy of <cit.> to simplify the spatial Ricci tensor.]
It is based on the Arnowitt-Deser-Misner (ADM) formulation of Einstein's 
equations <cit.>, and, like the ADM formulation, 
adopts a 3+1 foliation of spacetime <cit.>. Unlike the ADM formulation 
it also introduces a conformal-traceless decomposition as well as conformal connection 
functions (see also <cit.> for a textbook introduction).

To date, most numerical codes that adopt the BSSN formulation use finite 
differencing as well as Cartesian coordinates.
While Cartesian coordinates offer distinct advantages (most importantly, they are 
regular everywhere and do not feature any coordinate singularities),  
there are also several shortcomings: BHs, neutron stars, accretion
disks, etc., are often approximately 
spherical or axisymmetric, and Cartesian coordinates are not well suited 
to take advantage of these approximate symmetries. Furthermore, Cartesian
coordinates over-resolve angular directions at large distances,
which leads to the necessity of employing box-in-box mesh refinement. 

In large part, the main target of vacuum numerical relativity simulations are BBH
mergers, whose remnants are
Kerr <cit.> BHs. Being axisymmetric or nearly axisymmetric,
these merger remnants are prime targets for evolutions in spherical
coordinates. 
In <cit.>, the authors implemented the BSSN equations
in spherical coordinates without regularization in spherical symmetry.
An important ingredient in obtaining stability was the use of the partially 
implicit Runge-Kutta methods developed in <cit.>, even 
though it became clear later that stability can also be achieved with higher-order 
fully explicit Runge-Kutta methods <cit.>.
In <cit.>, the authors extended the evolution system described 
in <cit.> to full 3D and performed the first numerical relativity 
simulations in spherical coordinates without the assumption of any symmetries. 
The key idea in this approach is to handle the 
coordinate singularities at the origin and on the axis analytically, rather than 
numerically, which can be achieved with the help of a reference-metric formulation 
of the BSSN equations <cit.> 
together with a proper rescaling of all tensorial variables. The same methods can also be 
applied to relativistic hydrodynamics <cit.>. Several examples of vacuum and 
hydrodynamics simulations in spherical coordinates, including an off-center BH head-on 
collision, can be found in <cit.>, and simulations of critical collapse in 
the absence of spherical symmetry in <cit.>.
This approach has been generalized in the code <cit.> for 
various other curvilinear coordinate systems.

The use of spherical coordinates has clear advantages. Most importantly, the grid can 
take advantage of the approximate symmetries of the astrophysical objects to be simulated.  
Also, the number of angular grid points is independent of radius, while in Cartesian 
coordinates the number of points per great circle grows with distance from the origin. 
The unigrid (i.e. single computational domain without mesh refinement)
character of a spherical mesh does not produce short-wavelength noise as is the case 
for simulations with mesh refinement boundaries <cit.>. 
From a computational standpoint, our unigrid implementation in the offers 
another advantage: It is well documented that mesh refinement codes do not scale as well
as unigrid codes (see <cit.> for comparing scaling properties of 
unigrid and mesh refinement in the ).

However, these advantages compared to Cartesian coordinates come at a price: Spherical
coordinates have a well know limitation in the form of severely shorter time steps
due to the Courant-Friedrichs-Lewy (CFL) condition, as the cell volumes are not constant,
but decrease with increasing latitude towards the pole and decreasing radius towards the 
origin. A related issue is that the coordinate system becomes singular both at the 
origin and the polar axis, where coordinate values become multivalued. 

An approach to combine the best of both worlds is the use of multipatch
computational domains, in which the domain is broken up into several overlapping
patches locally adapted to the underlying symmetries of the physical system and
free of coordinate singularities and the time step limitations of spherical unigrid 
meshes <cit.>.
The SpeC code <cit.> uses such a 
multipatch grid structure <cit.>, but in the context
of a pseudospectral evolution scheme. Other  techniques in numerical
relativity codes for dealing with
the polar singularities  are the use of stereographic angular
coordinates, coupled to the eth formalism <cit.>, as is done in the PITT Null
code <cit.>, and the use of cubed spheres, as is
done in the Llama infrastructure <cit.>.
The use of multipatch grids, however, is not free of caveats either: Interpatch 
boundaries require interpolation of fields in ghost zones
which might introduce similar numerical noise as Cartesian mesh-refinement boundaries.

In this work, we report on an implementation of the BSSN equations in spherical coordinates 
described in <cit.> as a thorn called SphericalBSSN in the publicly available 
 <cit.>, using code for the BSSN equations provided by <cit.>.  
The was designed with Cartesian coordinates in mind, so that we had to adapt 
our implementation of spherical coordinates to its infrastructure in some regards. 
We first identify the x, y and z coordinates 
defined in the with r, θ and φ. The uses a vertex-centered grid
for finite differencing, meaning that grid points are placed on the edges of the 
physical domain.  This is not 
desirable in spherical coordinates, because grid points at the origin or on the axis 
would be singular.  We therefore move both the r and θ axes by half a grid point, 
so that, effectively, we implement a cell-centered grid in these directions (compare 
Fig. 1 in <cit.>). While, in Cartesian coordinates, the domain boundaries 
in x, y and z all correspond to outer boundaries, only the upper r domain 
boundary corresponds to an outer boundary in spherical coordinates. All other boundary 
conditions are inner boundary conditions. For the φ coordinate, these 
boundary conditions result from periodicity, while for the θ direction as well 
as at r = 0, the boundary conditions result from parity across the pole or the origin.  
For all inner boundaries, the ghost zones are filled in using properly identified 
interior grid points (see again Fig. 1 in <cit.> for an illustration), 
taking into account the parity of tensorial quantities.

In MPI-parallelized domain decompositions, the inner boundary conditions can require 
information from different processes and are therefore more difficult to implement than
in the context of OpenMP-parallelized, single-domain implementations. Using the 
Slab thorn <cit.>, we have
implemented the inner boundary conditions in an MPI-parallelized way, allowing
for arbitrary MPI domain decompositions. 
We also made several changes to existing diagnostics in the so that they can be 
used for evolutions in spherical coordinates, specifically the apparent horizon (AH) 
finder <cit.> and a thorn that 
computes quasilocal quantities <cit.> on AHs.
We test the new thorn, together with the changes in the existing diagnostics, 
for a single Bowen-York spinning BH <cit.>, which is equivalent to a Kerr BH 
with an axisymmetric Brill wave.

Throughout this paper we use units in which c=G=1, Latin indices run 
from 1 to 3, and the Einstein summation convention is used. 




§ IMPLEMENTATION IN THE EINSTEIN TOOLKIT



The is an open source code suite for relativistic astrophysics 
simulations. It uses the modular Cactus framework <cit.> 
(consisting of general modules called thorns) 
and provides adaptive mesh refinement (AMR) via the Carpet
driver <cit.>. Here we describe the implementation of the BSSN
evolution code in spherical coordinates within the Toolkit. All codes
mentioned here are either publicly available already, or are in the process of being 
released. 




 §.§ Evolution system



As outlined in <cit.>, the key 
idea in allowing stable evolutions of
the BSSN <cit.> equations in spherical coordinates 
is to treat the coordinate singularities analytically rather than numerically.  
Specifically, the equations contain terms that diverge with O(r^-2) close to 
the origin and O(sin^-2θ) close to the axis. Adopting a reference-metric 
formulation of the BSSN equations <cit.> 
together with a proper rescaling of all tensorial quantities, these terms can be 
differentiated analytically, and, for regular spacetimes, all code variables remain finite.
We also assume absence of conical singularities, which is sometimes 
referred to as  elementary flatness <cit.>.
This approach has been generalized in <cit.> for a larger number of 
curvilinear coordinate systems, such as spherical coordinates with a sinh(r) radial
coordinate and cylindrical coordinates, among others. 
We give a summary of the evolution system below, and refer the reader to the full details
in <cit.>.

Central to the method is  the conformally related spatial metric

    γ̅_ij = e^-4ϕγ_ij,

where γ_ij is the physical spatial metric, and ϕ the conformal factor 

    e^4ϕ = (γ / γ̅)^1/3,

where γ and γ̅ are the determinants of the physical and conformally 
related metric, respectively. In order to make the conformal rescaling unique, we adopt 
Brown's Lagrangian choice <cit.>

    ∂_t γ̅ = 0,

fixing γ̅ to its initial value throughout the evolution.
Similarly, the conformally related extrinsic curvature is defined as

    A̅_ij = e^-4ϕ(K_ij-1/3γ_ijK),

where K_ij is the physical extrinsic curvature and K=γ^abK_ab its trace.

The main idea is to write the conformally related metric as the sum of the flat background
metric plus perturbations (which need not be small)

    γ̅_ij = γ̂_ij + ϵ_ij,

where γ̂_ij is  the reference metric in spherical coordinates,

    γ̂_ij = [            1            0            0;            0          r^2            0;            0            0 r^2  sin^2 θ ],

and the corrections ϵ_ij are given by

    ϵ_ij = [               h_rr            r  h_rθ      r  sinθ  h_rφ;            r  h_rθ          r^2  h_θθ    r^2  sinθ  h_θφ;      r  sinθ  h_rφ    r^2  sinθ  h_θφ r^2  sin^2 θ  h_φφ ],

where h_ij is the rescaled evolved metric. This idea is similar to 
bimetric formalisms <cit.> in GR, in which reference metrics are
employed to give physical meaning to pseudotensors in curvilinear coordinates,
or in the integration of the Ricci scalar on a hypersurface <cit.>.
The evolved conformally rescaled extrinsic 
curvature a_ij is similarly related to the conformally related extrinsic curvature 
A̅_ij

    A̅_ij = [               a_rr            r  a_rθ      r  sinθ  a_rφ;            r  a_rθ          r^2  a_θθ    r^2  sinθ  a_θφ;      r  sinθ  a_rφ    r^2  sinθ  a_θφ r^2  sin^2 θ  a_φφ ].
 
The conformal connection coefficients Λ̅^i are treated as independent variables that
satisfy the initial constraint

    Λ̅^i - Δ^i = 0.
 
Here 

    Δ^i ≡γ̅^abΔ^i_ab

and Δ^i_jk is the difference between the Christoffel symbols of the
conformally rescaled and flat reference metric,

    Δ^i_jk≡Γ̅^i_jk - Γ̂^i_jk.

The conformal connection coefficients Λ̅^i therefore transform like vectors 
in the reference-metric formalism. Similar to our treatment of the metric and the extrinsic 
curvature we write 

    Λ̅^i = [            λ^r;        λ^θ / r; λ^φ / (r sinθ) ]

and evolve the variables λ^i in our code. We refer the reader 
to <cit.> for the full details of the evolution 
system.

The physical metric γ_ij and the physical extrinsic curvature 
K_ij can be reconstructed from the evolved variables h_ij and a_ij as follows:

    γ_ij = e^4ϕ[                 1+h_rr                r  h_rθ          r  sinθ  h_rφ;                r  h_rθ          r^2  (1+h_θθ)        r^2  sinθ  h_θφ;          r  sinθ  h_rφ        r^2  sinθ  h_θφ r^2  sin^2 θ  (1+h_φφ) ],

and

    K_ij = e^4ϕA̅_ij + 1/3γ_ij K.

Together with the lapse α and the shift β^i, this set of the 3+1 
variables {α,β^i,γ_ij,K_ij}, expressed in spherical coordinates,
is stored in the thorn ADMBase to interface with existing diagnostics in 
the  <cit.>. 

Numerical code for the evolution system is provided by the code, and the time
integration is performed with the method of lines as implemented in the 
MoL <cit.> thorn.




 §.§ Spherical parity boundary conditions



There is no global and regular one-to-one map from spherical to Cartesian coordinates.  
Instead, at least two charts are needed to cover an entire sphere of a given radius. 
Ultimately, this is due to the fact that θ and φ are multivalued at the 
coordinate origin, and φ is multivalued at the polar axis. As a result, 
the Jacobian from spherical to Cartesian coordinates diverges both at the origin and 
polar axis. 
In our implementation in the we use the existing Cartesian grid infrastructure as 
spherical coordinates by implementing the internal boundary conditions in a way that uses
the underlying topologically Cartesian grid.

As explained already in the Introduction, we start by identifying the internal 
(x,y,z) coordinate representation used in Carpet with the spherical coordinates 
(r,θ,φ).  Carpet uses a vertex-centered grid structure, meaning that 
grid points exist on the edges of the physical domains.  This is not desirable in 
spherical coordinates, because of the coordinate singularities at the origin, r = 0, 
and the poles at θ = 0 and θ = π.  We therefore shift both the r and 
θ axes by half a grid point. Therefore, the physical 3D domain has the following 
extents: 

    r    ∈   [dr/2,r_ max], 
    θ   ∈   [d θ/2,π-d θ/2], 
    φ   ∈    [0,2π - d φ].
  
Effectively, we therefore adopt a cell-centered grid in the r and θ directions, 
but maintain a vertex-centered grid in the φ direction.

Cartesian coordinates are topologically 𝐑^3, and all domain boundaries 
for large or small values of the coordinates x, y or z correspond to outer boundaries.  
In spherical coordinates, on the other hand, only r_ max corresponds to an 
outer boundary, while all other domain boundaries represent inner boundaries.   
At r = 0, for example, a radial grid line can be extended to negative values of r.  
We allow for n_ gzr ghost zone grid points at negative r; these ghost zone 
grid points correspond to interior grid points with positive r for some other values 
of the angles θ and φ (see Fig. 1 in <cit.> for an 
illustration). Specifically, we identify ghost zones at the origin with interior 
grid points at the coordinate locations

    r    →    - r 
    θ   →   π - θ, 
    φ   →   φ + π.

We can then fill these ghost zones by applying internal parity boundary conditions, 
which we explain in more detail below. Similarly, meridians, i.e., great circles 
of constant φ, can be extended across the pole, and the ghost zones there 
can again be identified with internal grid points. For θ_ min = 0 we have  

    θ   →    - θ, 
    φ   →   φ + π,

and for θ_ max = π 

    θ   →   π - θ, 
    φ   →   φ + π.

We also introduce ghost zones for 
φ, which can be filled by imposing periodicity. We note that application 
of this scheme requires an even number of grid points in the φ direction.

For scalar quantities, field values at interior grid points can be copied directly into 
the corresponding ghost zone grid points. For tensorial quantities, however, we 
have to take into account the fact that the direction of unit vectors changes when 
crossing the origin or pole (see <cit.> for more details).  
This observation leads to parity factors that arise in the application of the inner 
boundary conditions. We list these factors in Table <ref>.

By using the cell-centered grid in r and θ, and using the described 
internal boundary conditions, we are able to have a one-to-one mapping of the 
internal Cartesian coordinates in the , which are topologically 𝐑^3,
with the spherical coordinates used in the evolution.



Allowing for arbitrary MPI domain decompositions requires communication across
processes, as a given process might not be in possession of the point that is mapped to 
a ghost zone in its domain.  We have implemented these boundary conditions using the 
Slab thorn, which provides MPI-parallelized infrastructure to take 3D subarrays 
(slabs) of the 3D domain, manipulate them and then broadcast the manipulated 
slabs back to all processors that contain a part of it.  In what follows, we show how the 
internal boundary conditions are implemented as slab transfers using the
SLAB thorn (see Figs. <ref> and <ref>).






The source slab for the boundary condition at the origin contains the first n_ gzr physical 
points in r, where, again, n_ gzr is the number of ghost zones in the r direction, 
and all physical points 
in θ and φ. The operation θ→π - θ
is performed by inverting the slab in the θ direction, while 
the φ→φ + π operation corresponds to moving all points from 
[0,π] → [π,2 π] and, taking into account
the periodicity in φ, all points from [π,2 π] → [0,π]. The φ part of 
the operation is achieved using two separate calls to the slab transfer. Finally, the
source slab is inverted in r and the resulting slabs transferred into the ghost zones of
the domain. This is illustrated in Fig. <ref>. 
Note that the procedure only fills ghost zones that correspond to physical points in 
θ and φ, while double and triple ghost zones 
on the edges and corners of the computational grid need to be filled by subsequent
internal boundary conditions.

We then proceed by imposing the θ boundary conditions in a similar manner (see
schematic in Fig. <ref>), followed by 
applying periodic boundary conditions in φ using an existing thorn in the 
. This order ensures that all ghost zones 
that need to be specified by internal parity boundaries are filled correctly. 

In future applications that include magnetohydrodynamics and/or other matter fields, 
the same parity boundary conditions will apply to the matter fields as well.  
We have therefore implemented these boundary conditions in a separate thorn, 
SphericalBC, so that they are available for all evolved quantities, and not 
only for the spacetime evolution.




 §.§ Time Step considerations



It is well known that the time integration of hyperbolic PDEs in spherical coordinates 
suffers from severe CFL time step restrictions, as the cell sizes
become smaller and smaller with increasing latitude from the equator
towards the poles, and decreasing distance from the origin. In a spherical unigrid in flat spacetimes, 
the time step due to the CFL condition is given 
by <cit.>:

    dt = 𝒞  min[d r,d r/2 d θ,
    d r/2sin(d θ/2)d φ],
 
where the CFL factor 𝒞 is chosen between 0 and 1. The time step is therefore 
limited by the azimuthal distance between cells at the origin and polar axis.
Compared with Cartesian coordinates, where dt ≈ dx (when the same 
spatial resolution dx is used in all three coordinates), the time step in spherical
coordinates varies as dt ≈ d r   d θ  d φ.  
Thus high angular resolution will impose severe time step restrictions in spherical
coordinates. 

When using a high number of azimuthal cells, the CFL 
restriction might render the numerical integration 
computationally unfeasible. There are several approaches to mitigate this problem
(for an introduction, see e.g. <cit.>), from various multipatch 
approaches <cit.> to reducing the number of azimuthal cells at high 
latitudes as mesh coarsening in the azimuthal direction <cit.>, 
focusing resolution of the polar angle at the equator <cit.>,  
or the use of filters <cit.>, to name a few. 

To circumvent the severe time step limitation in cases when evolving BHs centered 
at the coordinate origin, we have devised a simple excision strategy in order to enlarge
time steps in these evolutions.
Specifically, we employ a radial extrapolation of all evolved variables deep inside 
the AH during the evolution, which essentially amounts to excising parts of the BH interior. 
Similar strategies have been employed and shown to work in the context of puncture
BH <cit.>.
Within a fixed number of radial points the BSSN variables are not evolved but rather
linearly extrapolated inwards radially from the first evolved points. When using this
technique, we get a dramatic increase in time step, which is now given by:

    dt = 𝒞  min[dr, A dr dθ,
    A dr sin(d θ/2) d φ],

where A ≡ n_ Ex-1+1/2 and n_ Ex is the number of
excised radial grid points. In the simulations presented in Sec. <ref>
below, we choose this parameter such that n_ Ex dr ≈ r_ AH/2 
initially, and a CFL factor 𝒞 = 0.4.  
We emphasize that we employ this excision for the purposes of speeding up the 
simulation only – it is not needed for stability. In Fig. <ref> below 
we compare simulations with and without this excision; we also note that the simulations
in <cit.>
did not use such an excision.




 §.§ Diagnostics



We use the AHFinderDirect thorn <cit.> to find
AHs <cit.> and the
QuasiLocalMeasures thorn <cit.> to calculate the angular momentum
of the apparent horizon during the evolution. The BH spin is measured using the 
flat space rotational Killing vector method <cit.> that was shown
to be equivalent to the Komar angular momentum <cit.> in foliations adapted
to the axisymmetry of the spacetime <cit.>.
Both thorns were written explicitly for the Cartesian coordinates employed in 
Carpet — interpolating entirely in the Cartesian basis both the 
ADMBase variables {γ_ij, K_ij, α, β^i}
and the partial derivatives of the spatial metric and extrinsic
curvature. As indicated in Sec. <ref> above, 
the ADMBase variables in the SphericalBSSN
thorn are the physical metric, extrinsic curvature, lapse and shift in spherical 
coordinates, which means we need to transform the ADMBase variables and their 
partial derivatives to Cartesian coordinates after interpolation. 
This is required because AHFinderDirect expects the computational domain to have
Cartesian topology (i.e., any surface with constant r, will not appear to be
closed to AHFinderDirect). In its original form, AHFinderDirect 
interacts with the rest of the Toolkit by requesting
the interpolation of the metric functions, and their derivatives, at
various points in Cartesian coordinates. To make this work with SphericalBSSN, we modify this behavior by transforming the Cartesian
coordinates to spherical (the necessary Jacobians are provided by aliased functions
defined in SphericalBSSN), and then after the interpolation step,
transforming the metric functions from spherical to Cartesian
coordinates using the Jacobian 
J^a_i≡∂ x^a/∂ x^i 
from spherical to Cartesian coordinates according to

    γ_ij   =    J^a_i  J^b_j γ_ab,
    γ_ij,k   =    J^a_i,c  J^b_j  J^c_k γ_ab,
       +    J^a_i  J^b_j,c  J^c_k γ_ab,
       +    J^a_i  J^b_j  J^c_k γ_ab,c,
    
     K_ij   =    J^a_i  J^b_j  K_ab,
    
     K_ij,k   =    J^a_i,c  J^b_j  J^c_k  K_ab,
       +    J^a_i  J^b_j,c  J^c_k  K_ab,
       +    J^a_i  J^b_j  J^c_k  K_ab,c,
    β^i    =    J^i_a β^a,

where we adopt the convention that indices a and b refer to spherical coordinates 
r, θ and φ, and indices i and j to the Cartesian 
coordinates x, y and z in coordinate transformations, and a comma 
indicates ordinary partial differentiation.
When the origin of the AHFinderDirect internal six-patch system coincides with the
coordinate origin, we add a small offset in θ at points located on the z-axis,
as the Jacobian diverges at those points.  

We extract GWs by computing the Weyl scalar Ψ_4 using the electric and magnetic 
parts of the Weyl tensor and constructing the numerical tetrad as described 
in <cit.> in spherical coordinates. The calculation of Ψ_4 and the
remaining Weyl scalars is contained in a new thorn called SphericalWeylScal.
The multipole expansion of the real and imaginary parts of Ψ_4 in spin-weighted 
spherical harmonics <cit.> is performed on the spherical grid used in the 
evolution. SphericalWeylScal performs the multipole expansion after the 
calculation of the Weyl scalars.
 
While the Jacobian for spherical coordinates is simple to implement
directly into the analysis thorns, we coded our
modification to AHFinderDirect  and QuasiLocalMeasures  so
that they call aliased functions. In this way, both codes can now work
with arbitrary coordinate systems, as the calculation of the Jacobians, etc.,
are handled by auxiliary routines.




 §.§ Initial Data



As a demonstration of our methods we show in Sec. <ref> below an evolution 
of a spinning Bowen-York BH <cit.>, which describes a perturbed Kerr BH <cit.>.

Bowen-York data are conformally flat, so that h_ij = 0 identically, as well as 
Λ̅^i = 0.  The data are also maximally sliced, so that K = 0.  
The momentum constraint can then be solved analytically for the conformally rescaled 
extrinsic curvature; for rotating Bowen-York BHs, the only non-vanishing component is 
the rφ component.  Given the analytical solution for this extrinsic curvature, 
the Hamiltonian constraint can then be solved numerically for the conformal factor e^ϕ.  
The only non-vanishing component of the extrinsic curvature variables a_ij 
defined in (<ref>) is then

    a_rφ = 3 e^-6 ϕ J sinθ/r^3,

where J is the magnitude of the BH's angular momentum 
(see also exercise 3.11 in <cit.>).

In order to allow for future applications with more general sets of initial data 
that may have been prepared in Cartesian coordinates, we do not implement the above 
results directly, but instead use the TwoPunctures thorn <cit.> to 
set up the data. 
This thorn uses spectral methods to solve the Einstein 
constraints, and interpolates the Cartesian 
ADMBase variables onto the computational mesh used 
in the simulation.  We have adapted the thorn to interpolate the Cartesian ADMBase 
variables onto the spherical grid points instead. Upon the completion of the interpolation, the metric
γ_ij and extrinsic curvature K_ij are transformed from Cartesian 
to spherical coordinates as described in Section <ref> above.  
The evolved metric variables h_ij are then computed from

    h_ij = e^-4ϕγ_ij ⊙ [                1              1/r      1/(r  sinθ);              1/r            1/r^2    1/(r^2  sinθ);      1/(r  sinθ)    1/(r^2  sinθ) 1/(r^2  sin^2 θ) ] - 1,

where ⊙ indicates the Hadamard product (element-wise matrix multiplication) and
1 the identity matrix, while the evolved extrinsic curvature variables a_ij are 

    a_ij =    e^-4ϕ(K_ij-1/3γ_ijK) 
       ⊙ [                1              1/r      1/(r  sinθ);              1/r            1/r^2    1/(r^2  sinθ);      1/(r  sinθ)    1/(r^2  sinθ) 1/(r^2  sin^2 θ) ].

We confirmed that these variables agree with the values listed above to within 
truncation error.

We complete the specification of the initial data with choices for the initial lapse 
and shift. We choose an initial shift β^i=0 and an initial lapse α = e^-2ϕ 
for the SphericalBSSN runs and β^i=0, α = 2 r/(M + 2r) for 
the comparison runs.




§ RESULTS








We perform simulations of a single, spinning, and initially conformally flat BH
(the Bowen-York solution <cit.>), with the initial data prepared as described 
in Section <ref> above. Since the Kerr <cit.> spacetime is not
conformally flat, these initial data
represent a spinning BH with gravitational wave content that will be radiated 
away <cit.>, allowing the BH to settle to the Kerr solution.
In all results presented here, dimensionful quantities are reported in terms of
M=1.  The BH has an initial spin J=0.8M^2 and an ADM mass <cit.> 
of M_ BH=1.18112M, giving a Kerr parameter a ≡ J/M_ BH = 0.677 M. 
We perform simulations of these initial data using our SphericalBSSN implementation.  
For comparison, we evolve the same initial data in Cartesian coordinates
using the McLachlan <cit.> thorn. McLachlan is a finite 
difference code generated using Kranc <cit.> 
that solves the BSSN equations as part of the . For our comparisons here we use 
4th-order spatial finite differences in both codes, but we note that and
McLachlan are capable of providing finite-difference stencils for the BSSN equations 
at arbitrary order. We summarize the details of the relevant simulation parameters in 
Table <ref>.




 §.§ BH mass and spin



In Figs. <ref> and <ref> we plot the evolution of the irreducible mass
of the BH and its angular momentum, respectively. During the first 35 M of the evolution 
the BH mass increases due to the absorption of some of the GW content in the 
spacetime (see Fig. <ref> below). We omit this initial time in 
Fig. <ref>, and instead show the long-term behavior after the BH has settled down.  
We show results for two different θ resolutions (N_θ = 32 and 64) with 
SphericalBSSN and two resolutions (dx=0.0125 and 0.02 on the finest mesh) using 
the in Cartesian coordinates with box-in-box mesh refinement. The radial resolution in 
both evolutions (dr=0.02 and dx=0.02 on the finest Cartesian mesh which covers the AH) 
gives approximately 25 radial points across the minimum diameter (0.25) of the AH initially. 
For the irreducible mass shown in Fig. <ref>, the results obtained with the 
higher resolution SphericalBSSN and the two Cartesian runs agree well, while 
the lower resolution SphericalBSSN run appears to be under-resolved, showing a 
linear growth in the irreducible mass that is unphysical. 
There is a notable absence of oscillations in the higher resolution SphericalBSSN 
run, which can be seen in both Cartesian runs (converging away with increasing 
resolution).
 
The evolution of the angular momentum of the AH, shown in Fig. <ref>, exhibits a 
similar behavior. The two high-resolution runs with McLachlan and SphericalBSSN
perform similarly, while the lower resolution runs show linear drifts in both Cartesian 
and spherical coordinates. The Cartesian simulations show larger initial oscillations
that do not seem to converge away with increasing resolution, likely due to reflections 
of short-wavelength modes across mesh boundaries <cit.>. 
Just as for the irreducible mass, the high-resolution SphericalBSSN simulations 
performs best. 

To test the effect of the excision described in Section <ref>, we plot
the initial evolution of the irreducible mass for a run with, and one
without, our excision procedure in Fig. <ref>.   
Evidently, the excision procedure does not have any visible effect on 
the accuracy or stability of our method.





 §.§ Gravitational waves



In Fig. <ref>, we plot the absolute value of the multipole
expansion in spin-weighted spherical harmonics _-2Y_lm of the Weyl scalar Ψ_4,
extracted at a radius of r=180, for the even l=2 through 8, m=0 multipoles
(by symmetry, only the m=0 modes are non-zero). The plot shows the 
Cartesian simulation with dx=0.02 and two different θ-resolutions for the 
spherical simulations. There are notable differences between the Cartesian and spherical
evolutions: In the spherical simulations, there is an absence of
initial noise pulse
before the radiation reaches it peak value, and the decay after the peak value proceeds
much cleaner and to orders of magnitude below the values attained in the Cartesian 
simulation. The reason for this difference in behavior is the
fact that there are partial reflections of the outgoing wave at each Cartesian mesh 
refinement boundary (see, e.g. <cit.>), causing the 
unphysical excitation of l≥4 multipole modes, as well as reflections in the initial 
noise pulse in seen in the l=2 mode.  These reflections affect strong-field quantities 
as well, as described in <cit.>. Contrary to this, the
spherical grid in SphericalBSSN is a single uniform grid, so there is a complete 
absence of these reflections (apart from reflections from the outer boundary), leading to much 
cleaner signals especially in the higher-order multipoles.  








 §.§ Kerr quasinormal modes




Given that the SphericalBSSN simulation provides us with very accurate
higher order modes, we turn our attention to an analysis of the BH's
quasinormal modes. As explained in the setup of the initial data, the spinning,
initially conformally flat BH should settle down to a Kerr BH via
ringdown of its quasinormal modes (QNM) (see <cit.> for a review).
In Fig. <ref>, we plot the l=2 through 8, m=0 modes of |Ψ_4|
for the high resolution N_θ=64 simulation alone.
We note that even (odd) l-modes contain only the real (imaginary) part
of Ψ_4. The l=2 mode follows a clear
exponential decay, but the higher-order modes exhibit a beating 
modulation on top of the exponential decay.  The reason for this behavior is
that the quasinormal modes for Kerr are defined in terms of spheroidal
harmonics, _sS_lm, while we decompose the waveform in terms of 
spin-weighted spherical harmonics _sY_lm.
Following <cit.>, we can decompose Ψ_4 in terms of spin-weighted 
spheroidal harmonics _sS_l'm with s = -2 according to

    Ψ_4 = ∑_l',m',n' A_l'm'n' _sS_l'm' e^-α_l'm'n' t e^i ω_l'm'n' t
 
where n is the overtone number of each mode, α_l'mn is its the decay rate, 
ω_l'mn its frequency, and the coefficients A_l'mn are the amplitudes of 
the individual modes.  
In particular, we see that each mode oscillates and decays at well-defined rates.  
In practice, however, Ψ_4 is projected into the spin-weighted spherical 
harmonics _sY_lm, i.e.

    Ψ_4^lm    =    ∫Ψ_4 _s Y^*_lm dΩ
        =    ∑_l',m',n' A_l'm'n' e^-α_l'm'n' t e^i ω_l'm'n' t∫_sS_l'm'_sY^*_lm dΩ
        =    ∑_l',m',n' A_l'm'n' e^-α_l'm'n' t e^i ω_l'm'n' tμ^*_mll'n'(a) δ_mm'
        =    ∑_l',n' A_l'mn' e^-α_l'mn' t e^i ω_l'mn' tμ^*_mll'n'(a).

Here the coefficients μ^*_mll'n'(a) describe the mixing between spin-weighted
spheroidal and spherical harmonics; they are defined in eq. (5) of <cit.> and depend on 
the black hole Kerr parameter a.

Ignoring higher-order overtone modes with n' > 0, which decay faster than 
the fundamental modes, we fit the l=3 through 8, m=0 spherical harmonic modes 
computed from our numerical data to the form

    Ψ^l0_4(t) ≈ ∑_l^'=2^l^'=10 
    A_l^'0 e^(-α_l^'00 t) sin (ω_l^'00 t + ϕ_l^'0),

where the unknowns A_l^'0 and
ϕ_l^'0 serve as 18 parameters corresponding to the amplitudes 
(including the mixing coefficients) and phases, respectively.
We fix the α_l^'00 and ω_l^'00 in the fit to be the
values of the decay rate and frequency that correspond to the Kerr BH in our
simulation (J/M_ BH^2=0.573 and M_ BH=1.18112), using the tabulated values 
and Mathematica notebooks to calculate QNMs <cit.> found at <cit.>. 
The results of the fit for the l=3 through 8 modes are shown in 
Fig. <ref> for a fitting window of t=220-340. 
The beating of the modes is very well captured by modeling a given l mode as the sum 
of the expected decay rates and frequencies calculated in spin-weighted spheroidal harmonics, 
showing that mode mixing is responsible for the observed beating. A similar type of 
equal m mode mixing has been observed in <cit.>.




§ DISCUSSION



We report on an implementation of the BSSN equations in spherical coordinates in 
the .  While Cartesian coordinates have advantages for many applications, 
spherical coordinates are much better suited to take advantage of the approximate 
symmetries in many astrophysical systems. The problems associated with the coordinate 
singularities that appear in curvilinear coordinates can be avoided if these singularities 
are treated analytically – which, in turn, is possible with the help of a reference-metric 
formulation of the BSSN equations <cit.> 
and a proper rescaling of all tensorial quantities <cit.>. 
We implement this formalism in the in an effort to make these techniques publicly 
available to the entire numerical relativity community and beyond.

Specifically, we adapt the infrastructure, which originally was designed 
for Cartesian coordinates, for spherical coordinates. In contrast to Cartesian 
coordinates, spherical coordinates feature inner boundary condition, where ghost 
zones are filled by copying interior data from other parts of the numerical grid, 
taking into account proper parity conditions. We implemented these boundary conditions, 
which may require communication across processors, within an MPI-parallelized 
infrastructure using the Slab thorn. Numerical code for the BSSN equations in 
spherical coordinates were provided by  <cit.>.

In order to test and calibrate our implementation we performed simulations of a single, 
spinning and initially conformally flat BH, and compared the evolution of BH mass, 
spin and GWs using our spherical BSSN (SphericalBSSN) and Cartesian AMR BSSN 
(McLachlan) code with comparable grid resolutions.
For sufficiently high resolutions, the evolution on a spherical mesh conserves
irreducible mass and angular momentum far better than with Cartesian
AMR. In particular, there are no reflections
of the initial junk radiation or outgoing initial gauge
pulse <cit.> at mesh refinement boundaries, causing
the evolution of irreducible mass and spin to be smoother in the unigrid spherical 
evolution. 
The advantage of using unigrid spherical coordinates over Cartesian coordinates with
box-in-box mesh refinement becomes particularly apparent when
analyzing the higher-order l-multipoles of the GW wave signal. These
signals are affected by partial reflections at mesh refinement 
boundaries, leading to a contamination of all higher order l,m modes that never fully
leaves the computational domain. This effect is completely absent in the simulations using
the SphericalBSSN thorn, where the quasinormal ringdown of the Kerr BH is 
observed to much smaller amplitudes than in the Cartesian simulations. We observe a 
significant beating of the exponential ringdown of multipoles with l > 2,
which can be explained by spheroidal-spherical multipole mode mixing. The accurate modeling
of the ringdown of higher order modes is necessary in order to provide
GW detectors with accurate templates <cit.>, as the
measurement of two quasinormal modes is needed to test 
the no hair theorem <cit.>.

The current SphericalBSSN thorn adopts uniform resolution in radius, which requires
a large number of points in order to place the outer boundary sufficiently far away and
to avoid contaminating the inner parts of the computational domain with noise 
from the outer boundary (i.e., to causally disconnect these inner parts from the outer 
boundary). Possible approaches to improve this is to adopt a non-uniform radial grid, 
e.g. a logarithmic grid as implemented in <cit.> 
or <cit.>, or to use more general 
radial coordinates. The code <cit.> allows for 
such generalized radial coordinates – a convenient choice is sinh(r) – and 
we plan to port these features into the SphericalBSSN thorn in the
future.

We also plan to supplement our current implementation of Einstein's vacuum equations 
in spherical coordinates with methods for relativistic hydrodynamics and magnetohydrodynamics
as another set of publicly available thorns for the .  
As shown in <cit.>, these equations can be expressed with the 
help of a reference-metric as well.  Further using a rescaling of all tensorial quantities 
similar to the rescaling of the gravitational field quantities in this paper, the evolution 
of hydrodynamical variables is unaffected by the coordinate singularities.  We hope that 
with these methods, and possibly implementations of microphysical processes like radiation 
transport and nuclear reaction chains, the in spherical coordinates will become a 
powerful and efficient community tool for fully relativistic simulations of a number of 
different objects, including rotating neutron stars, gravitational collapse, accretion 
disks, and supernova explosions.   
We believe that this will result in a new open source state-of-the-art code 
that will prove to be a valuable resource for a broad range of future simulations. 



  The authors would like to thank Emanuele Berti for useful discussions and
  Dennis B. Bowen for a careful reading of the manuscript.
  We gratefully acknowledge the National Science Foundation (NSF) for
  financial support from Grants No.  OAC-1550436, AST-1516150, PHY-1607520, PHY-1305730, 
  PHY-1707946,  PHY-1726215, to RIT, as well as Grants No.  PHYS-1402780 and PHYS-1707526 to Bowdoin College.
  V.M. also acknowledges partial support from AYA2015-66899-C2-1-P, and RIT for the 
  FGWA SIRA initiative.  
  This work used the Extreme Science and Engineering Discovery Environment (XSEDE) 
  [allocation TG-PHY060027N], which is supported by NSF grant No. ACI-1548562, and by the 
  BlueSky Cluster at RIT, which is supported by NSF grants AST-1028087, PHY-0722703, and PHY-1229173.
  Funding for  computer equipment to support the development of was provided in part by NSF EPSCoR Grant OIA-1458952 to West Virginia University.
  Computational resources were also provided by the Blue Waters sustained-petascale computing 
  NSF project OAC-1516125. 




apsrev4-1 


