$\mathtt{bimEX}$: A Mathematica package for exact computations in $3+1$ bimetric relativity

TL;DR

bimEX is a Mathematica package designed for exact symbolic and component-based computations in 3+1 bimetric relativity, facilitating the derivation of complex equations in this field.

Contribution

It introduces a new Mathematica package that streamlines handling multiple metrics and deriving bimetric equations in specific charts, enhancing computational efficiency and accuracy.

Findings

Successfully derived covariant BSSN bimetric equations

Reproduced 3+1 equations in spherical polar coordinates

Validated package performance in spherical symmetry scenarios

Abstract

We present $\mathtt{bimEX}$, a Mathematica package for exact computations in 3$+$1 bimetric relativity. It is based on the $\mathtt{xAct}$ bundle, which can handle computations involving both abstract tensors and their components. In this communication, we refer to the latter case as concrete computations. The package consists of two main parts. The first part involves the abstract tensors, and focuses on how to deal with multiple metrics in $\mathtt{xAct}$. The second part takes an ansatz for the primary variables in a chart as the input, and returns the covariant BSSN bimetric equations in components in that chart. Several functions are implemented to make this process as fast and user-friendly as possible. The package has been used and tested extensively in spherical symmetry and was the workhorse in obtaining the bimetric covariant BSSN equations and reproducing the bimetric 3$+$1 equations in the spherical polar chart.

Equations

$$\displaystyle{\color[rgb]{0,0,0}g}$$

$$\displaystyle{\color[rgb]{0,0,0}f}$$

$$\displaystyle{\color[rgb]{0,0,0}K}_{ij}\coloneqq-\dfrac{1}{2}\mathcal{L}_{{\color[rgb]{0,0,0}n}}\,{\color[rgb]{0,0,0}\gamma}_{ij},\quad{\color[rgb]{0,0,0}\mkern 1.0mu\widetilde{\mkern-1.0muK}}_{ij}\coloneqq-\dfrac{1}{2}\mathcal{L}_{{\color[rgb]{0,0,0}\mkern 1.0mu\widetilde{\mkern-1.0mun}}}\,{\color[rgb]{0,0,0}\varphi}_{ij},$$

$$\displaystyle{\color[rgb]{0,0,0}h}\coloneqq{\color[rgb]{0,0,0}g}\operatorname{\#}{\color[rgb]{0,0,0}f}={\color[rgb]{0,0,0}E}^{\mathtt{{\scriptscriptstyle\top}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{\eta}}}L{\color[rgb]{0,0,0}M_{\mathtt{o}}}={\color[rgb]{0,0,0}h}^{\mathtt{{\scriptscriptstyle\top}}},$$

$$\displaystyle L=\begin{pmatrix}{\color[rgb]{0,0,0}\lambda}&{\color[rgb]{0,0,0}\boldsymbol{\mathrm{p}}}^{\mathtt{{\scriptscriptstyle\top}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}}\\ {\color[rgb]{0,0,0}\boldsymbol{\mathrm{p}}}&{\color[rgb]{0,0,0}\boldsymbol{\Lambda}_{\mathrm{s}}}\end{pmatrix}\begin{pmatrix}1&0\\ 0&{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}}\end{pmatrix}=\begin{pmatrix}{\color[rgb]{0,0,0}\lambda}&{\color[rgb]{0,0,0}\boldsymbol{\mathrm{p}}}^{\mathtt{{\scriptscriptstyle\top}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}}\\ {\color[rgb]{0,0,0}\boldsymbol{\mathrm{p}}}&{\color[rgb]{0,0,0}\boldsymbol{\Lambda}_{\mathrm{s}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}}\end{pmatrix},$$

$$\displaystyle{\color[rgb]{0,0,0}\beta}={\color[rgb]{0,0,0}q}+\dfrac{{\color[rgb]{0,0,0}\alpha}}{{\color[rgb]{0,0,0}\lambda}}{\color[rgb]{0,0,0}e}^{-1}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{p}}},\qquad{\color[rgb]{0,0,0}\mkern 1.0mu\widetilde{\mkern-1.0mu\beta}}={\color[rgb]{0,0,0}q}-\dfrac{{\color[rgb]{0,0,0}\mkern 1.0mu\widetilde{\mkern-1.0mu\alpha}}}{{\color[rgb]{0,0,0}\lambda}}{\color[rgb]{0,0,0}m}^{-1}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{p}}},$$

$$\displaystyle{\color[rgb]{0,0,0}\chi}={\color[rgb]{0,0,0}e}^{\mathtt{{\scriptscriptstyle\top}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\boldsymbol{\Lambda}_{\mathrm{s}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}}{\color[rgb]{0,0,0}m_{\mathtt{o}}}=({\color[rgb]{0,0,0}e}^{\mathtt{{\scriptscriptstyle\top}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\boldsymbol{\Lambda}_{\mathrm{s}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}}{\color[rgb]{0,0,0}m_{\mathtt{o}}})^{\mathtt{{\scriptscriptstyle\top}}}={\color[rgb]{0,0,0}\chi}^{\mathtt{{\scriptscriptstyle\top}}},$$

$$\displaystyle{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}}$$

$$\displaystyle{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}}$$

$$\displaystyle{\color[rgb]{0,0,0}{\widebar{\gamma}}}_{ij}$$

$$\displaystyle{\color[rgb]{0,0,0}{\widehat{\varphi}}}_{ij}$$

$$\displaystyle{\color[rgb]{0,0,0}{\widebar{A}}}_{ij}$$

$$\displaystyle{\color[rgb]{0,0,0}{\widehat{A}}}_{ij}$$

$$\displaystyle{\color[rgb]{0,0,0}{\widebar{K}}}$$

$$\displaystyle{\color[rgb]{0,0,0}{\widebar{\Lambda}}}^{i}\coloneqq{\color[rgb]{0,0,0}\gamma}^{jk}{\color[rgb]{0,0,0}\bigtriangleup{\widebar{\Gamma}}}^{i}_{jk}={\color[rgb]{0,0,0}\gamma}^{jk}\left({\color[rgb]{0,0,0}{\widebar{\Gamma}}}^{i}_{jk}-{\color[rgb]{0,0,0}{\widebar{\Gamma}}_{\mathtt{B}}}^{i}_{jk}\right),$$

$$\displaystyle{\color[rgb]{0,0,0}{\widehat{\Lambda}}}^{i}\coloneqq{\color[rgb]{0,0,0}\varphi}^{jk}{\color[rgb]{0,0,0}\bigtriangleup{\widehat{\Gamma}}}^{i}_{jk}={\color[rgb]{0,0,0}\varphi}^{jk}\left({\color[rgb]{0,0,0}{\widehat{\Gamma}}}^{i}_{jk}-{\color[rgb]{0,0,0}{\widehat{\Gamma}}_{\mathtt{B}}}^{i}_{jk}\right),$$

$$\displaystyle{\color[rgb]{0,0,0}{\accentset{\circ}{{\chi}}}}_{ij}$$

$$\displaystyle{\color[rgb]{0,0,0}\phi},{\color[rgb]{0,0,0}\psi},{\color[rgb]{0,0,0}{\widebar{e}}}^{\textbf{a}}{}_{i},{\color[rgb]{0,0,0}{\widehat{m}}_{\mathrm{o}}}^{\textbf{a}}{}_{i},{\color[rgb]{0,0,0}{\widebar{A}}}^{i}{}_{j},{\color[rgb]{0,0,0}{\widehat{A}}}^{i}{}_{j},{\color[rgb]{0,0,0}{\widebar{\Lambda}}}^{i},{\color[rgb]{0,0,0}{\widehat{\Lambda}}}^{i},{\color[rgb]{0,0,0}\boldsymbol{\mathrm{p}}}^{\textbf{a}},{\color[rgb]{0,0,0}q}^{i},{\color[rgb]{0,0,0}{\widebar{\Gamma}}_{\mathtt{B}}}^{i}_{jk},{\color[rgb]{0,0,0}{\widehat{\Gamma}}_{\mathtt{B}}}^{i}_{jk},$$

$$\displaystyle{\color[rgb]{0,0,0}{\widebar{e}}}=\mathrm{e}^{-2{\color[rgb]{0,0,0}\phi}}{\color[rgb]{0,0,0}e},\quad{\color[rgb]{0,0,0}{\widehat{m}}}=\mathrm{e}^{-2{\color[rgb]{0,0,0}\psi}}{\color[rgb]{0,0,0}m}.$$

$$\underline{f\text{textunderscore}}:=f[⟨\mathtt{detected}\mathtt{or}\mathtt{specified}\mathtt{variables}⟩].$$

$$\displaystyle{\color[rgb]{0,0,0}{\widebar{\eulR}}}_{ij}$$

$$\displaystyle\quad\ -{\color[rgb]{0,0,0}{\widebar{\gamma}}}^{k\ell}{\color[rgb]{0,0,0}{\widebar{\gamma}}}_{m(i}{\color[rgb]{0,0,0}{\widebar{R}}_{\mathtt{B}}}_{j)k\ell}{}^{m}+{\color[rgb]{0,0,0}{\widebar{\gamma}}}^{\ell m}{\color[rgb]{0,0,0}\bigtriangleup{\widebar{\Gamma}}}^{k}_{\ell m}{\color[rgb]{0,0,0}\bigtriangleup{\widebar{\Gamma}}}_{(ij)k}$$

$$\displaystyle\quad\ +{\color[rgb]{0,0,0}{\widebar{\gamma}}}^{k\ell}\left(2{\color[rgb]{0,0,0}\bigtriangleup{\widebar{\Gamma}}}^{m}_{k(i}{\color[rgb]{0,0,0}\bigtriangleup{\widebar{\Gamma}}}_{j)m\ell}+{\color[rgb]{0,0,0}\bigtriangleup{\widebar{\Gamma}}}^{m}_{ik}{\color[rgb]{0,0,0}\bigtriangleup{\widebar{\Gamma}}}_{mj\ell}\right),$$

$$\displaystyle{\color[rgb]{0,0,0}{\widehat{\eulR}}}_{ij}$$

$$\displaystyle\quad\ -{\color[rgb]{0,0,0}{\widehat{\varphi}}}^{k\ell}{\color[rgb]{0,0,0}{\widehat{\varphi}}}_{m(i}{\color[rgb]{0,0,0}{\widehat{R}}_{\mathtt{B}}}_{j)k\ell}{}^{m}+{\color[rgb]{0,0,0}{\widehat{\varphi}}}^{\ell m}{\color[rgb]{0,0,0}\bigtriangleup{\widehat{\Gamma}}}^{k}_{\ell m}{\color[rgb]{0,0,0}\bigtriangleup{\widehat{\Gamma}}}_{(ij)k}$$

$$\displaystyle\quad\ +{\color[rgb]{0,0,0}{\widehat{\varphi}}}^{k\ell}\left(2{\color[rgb]{0,0,0}\bigtriangleup{\widehat{\Gamma}}}^{m}_{k(i}{\color[rgb]{0,0,0}\bigtriangleup{\widehat{\Gamma}}}_{j)m\ell}+{\color[rgb]{0,0,0}\bigtriangleup{\widehat{\Gamma}}}^{m}_{ik}{\color[rgb]{0,0,0}\bigtriangleup{\widehat{\Gamma}}}_{mj\ell}\right),$$

Full text

††journal: Computer Physics Communications

PROGRAM SUMMARY.

Program Title: $\mathtt{bimEX}$

Program Files doi: http://dx.doi.org/10.17632/2s5d7csc9w.1

Licensing provisions: GNU General Public License 3.0 (GPL)

Programming language: $\mathtt{Mathematica}$

Supplementary material:

README file, containing instructions about how to use the working example. 2. 2.

Working example, constituted by the notebooks:

1. (a)

bimEX_Working_Example.nb 2. (b)

bimEX_Decomposition_Lists_Loader.nb 3. (c)

bimEX_Decomposition_xAct_Loader.nb

*Nature of problem (approx. 50-250 words):

*Writing the bimetric covariant BSSN equations in any desired ansatz and chart.

*Solution method (approx. 50-250 words):

*Definition of functions within the $\mathtt{Mathematica}$ package $\mathtt{xAct}$, which computes all the components of the defined abstract tensors and reduce the abstract tensors to their representation in components.

*Additional comments including Restrictions and Unusual features (approx. 50-250 words):

*GitHub repository at https://github.com/nubirel/bimEX

1 Introduction

1.1 Motivation and general description

The Hassan–Rosen bimetric theory, or bimetric relativity (BR), is a nonlinear theory of interacting massless and massive spin-2 fields [1, 2, 3, 4]. As a theory of modified gravity, it has a rich phenomenology, and its spectrum of solutions contain both the general relativity (GR) solutions and novel solutions [5, 6]. We refer the reader to [7] for a review on the theory, to [8] for a more recent non-GR exact solution of the theory, and to [9] for a recent work affirming the compatibility of the theory with local tests of gravity. However, the number of non-GR solutions is modest, at present. This is due to the difficulties in solving the bimetric field equations (BFE), both analytically and numerically. Most of the non-GR solutions, see, e.g., [10, 11, 12, 13] and [14], have been obtained by integrating the BFE numerically after reducing them to ordinary differential equations.

In the search for solutions describing more realistic physical systems, e.g., spherically symmetric vacuum and non-vacuum solutions with non-trivial dynamics, one has to deal with a system of partial differential equations (PDEs). In GR, the Einstein field equations are also PDEs in the general case, and their numerical integration needs a recasting as a well-posed Cauchy problem. See [15] for a review on the history of the Cauchy problem in GR. The same thing holds in BR.

The recasting of the BFE as a Cauchy problem was established in [16]. However, this recasting does not result in a well-posed formulation. For this reason, following the road suggested by numerical relativity, it is desirable to recast the BFE in the covariant BSSN formalism (see [17, 18] and [19] for the BSSN formulation of the Einstein field equations), which is well-posed in GR if one chooses the standard gauge and satisfies other technical conditions [20, 21, 22]. In [23], the covariant BSSN formulation of the BFE is computed. Unfortunately, due to the particular interaction between the metrics, the well-posedness of the bimetric covariant BSSN has not been established yet.

Having found the bimetric covariant BSSN formulation, one would like to be able to compute it in any desired chart. This is the first step towards the numerical integration of the equations and the attainment of solutions to the BFE describing sensible bimetric physical systems, e.g., nonlinear bimetric gravitational collapse. $\mathtt{bimEX}$ was developed with this goal in mind, and was the workhorse in obtaining the results in [23, 24], and reproducing the results in [16]. The name of the package is an acronym for “bimetric exact computations”. The last part of the name, $\mathtt{EX}$, can be interpreted as “ex” in exact, or “e” in exact and “cs”($=$x) from computations.

In this communication we describe $\mathtt{bimEX}$, developing a working framework for $\mathtt{Mathematica}$ 11.0 [25] or later, and $\mathtt{xAct}$ 1.1.3 [26] or later, to perform computations in $3+1$ bimetric relativity. $\mathtt{xAct}$ already provides an excellent framework to handle both abstract and concrete computations; with “concrete” equations, we mean equations written in components in some chart. $\mathtt{bimEX}$ adds to it the definitions of several bimetric geometrical objects defined in [16] and [23], the definitions of the bimetric covariant BSSN constraint and evolutions equations introduced in [23] and the definitions of several ready-to-use functions that allow the user to obtain the concrete equations starting from the abstract ones, in any given chart. For example, by using $\mathtt{bimEX}$, the user can call only one function to write an abstract equation into components, without the need to explicitly call the $\mathtt{xAct}$ function $\mathtt{ToBasis}$, $\mathtt{TraceBasisDummy}$ and the others. The latter are used in the background by $\mathtt{bimEX}$. Also, another $\mathtt{bimEX}$ function takes the chosen ansatz on the primary variables—introduced in the next section—as its input and gives the bimetric decomposition as output, without the need for the user to explicitly write all the needed formulas defined in [16, 23].

All the commands and options in $\mathtt{bimEX}$ are documented through usage messages in $\mathtt{Mathematica}$, i.e., by typing $\mathtt{?\langle name\;of\;the\;function/option\rangle}$ after loading the package.

1.2 Bimetric relativity and its primary variables

In this section, we very briefly introduce the features of bimetric relativity needed to understand the structure of $\mathtt{bimEX}$, following the notation in [23].

In bimetric relativity, the two spin-2 fields are described by two metrics, {\color[rgb]{0,0,0}g} and {\color[rgb]{0,0,0}f}. Note, however, that the metrics do not coincide with the mass eigenstates of the theory [27]. It is possible to write each one of the two metrics in terms of dynamical and kinematical variables as [28, 29, 30] (see also [19, 31]),

[TABLE]

Here, the lapse functions {\color[rgb]{0,0,0}\alpha},{\color[rgb]{0,0,0}\mkern 1.0mu\widetilde{\mkern-1.0mu\alpha}} and the shift vectors {\color[rgb]{0,0,0}\beta},{\color[rgb]{0,0,0}\mkern 1.0mu\widetilde{\mkern-1.0mu\beta}} [32] are kinematical variables, as can be shown by performing the Hamiltonian analysis of the theory [1, 4]. The spatial metrics {\color[rgb]{0,0,0}\gamma}_{ij},{\color[rgb]{0,0,0}\varphi}_{ij} are instead dynamical, and are the induced metrics on a common spacelike hypersurface $\Sigma_{t}$ [19, 31, 1]. When writing the BFE in terms of these variables, one obtains a set of constraint equations—in the usual Hamiltonian sense—and a set of evolution equations. These equations can be written in different forms, one of them being the so-called standard $N+1$ decomposition [19, 31]. The bimetric standard $N+1$ decomposition was computed in [16].

In the standard $N+1$ decomposition, both in GR and in BR, the evolution equations are written in such a way that they contain only first-order time derivatives. In order to do that, one has to promote the time derivative of the metric to be a dynamical variable. This is done by introducing the extrinsic curvatures [19, 31],

[TABLE]

where $\mathcal{L}_{X}$ is the Lie derivative along the vector field $X$, {\color[rgb]{0,0,0}n}=(\partial_{t}-{\color[rgb]{0,0,0}\beta})/{\color[rgb]{0,0,0}\alpha} and {\color[rgb]{0,0,0}\mkern 1.0mu\widetilde{\mkern-1.0mun}}=(\partial_{t}-{\color[rgb]{0,0,0}\mkern 1.0mu\widetilde{\mkern-1.0mu\beta}})/{\color[rgb]{0,0,0}\mkern 1.0mu\widetilde{\mkern-1.0mu\alpha}} are the two normal vectors to the spacelike hypersurface $\Sigma_{t}$, with respect to {\color[rgb]{0,0,0}g} and {\color[rgb]{0,0,0}f}.

In the bimetric parametrization established in [3, 16], it is made clear that, given any two metrics in (1), the real square root ({\color[rgb]{0,0,0}g}^{-1}{\color[rgb]{0,0,0}f})^{1/2} does not necessarily exist, although necessary to be able to write down the theory, since the interaction potential between the metrics depends on it.

A necessary and sufficient condition for the real square root to exist is that it is possible to find a Lorentz transformation $L=\Lambda\,R$, with $\Lambda$ boost and $R$ spatial rotation, such that the geometric mean of the two metrics,

[TABLE]

exists [33, 3]. Here {\color[rgb]{0,0,0}E},{\color[rgb]{0,0,0}M_{\mathtt{o}}} are the vielbeins of {\color[rgb]{0,0,0}g},{\color[rgb]{0,0,0}f}, respectively, and {\color[rgb]{0,0,0}\boldsymbol{\mathrm{\eta}}} is the Minkowski metric.

The Lorentz transformation $L$ can be written [34, Sec. 1.6],

[TABLE]

with {\color[rgb]{0,0,0}\lambda}=\sqrt{1+{\color[rgb]{0,0,0}\boldsymbol{\mathrm{p}}}^{\mathtt{{\scriptscriptstyle\top}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{p}}}} being the Lorentz factor of the boost, and {\color[rgb]{0,0,0}\boldsymbol{\mathrm{p}}} is a real spatial vector in the Lorentz frame, called “separation parameter”, with components {\color[rgb]{0,0,0}\boldsymbol{\mathrm{p}}}^{\textbf{a}}=\sinh({\color[rgb]{0,0,0}\boldsymbol{\mathrm{w}}}^{\textbf{a}}), with {\color[rgb]{0,0,0}\boldsymbol{\mathrm{w}}}^{\textbf{a}} rapidities of the Lorentz boost. The condition (3), in the $N+1$ formalism, translates into a condition on the shifts of the two metrics, and another condition on the spatial part {\color[rgb]{0,0,0}\boldsymbol{\Lambda}_{\mathrm{s}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}} of $L$ [35]. The first condition tells us that the two shifts are related,

[TABLE]

where {\color[rgb]{0,0,0}q} is the shift vector of {\color[rgb]{0,0,0}h}. The second condition requires the spatial part {\color[rgb]{0,0,0}\chi} oh {\color[rgb]{0,0,0}h} to be symmetric,

[TABLE]

where {\color[rgb]{0,0,0}e},{\color[rgb]{0,0,0}m_{\mathtt{o}}},{\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}} are the spatial parts of {\color[rgb]{0,0,0}E},{\color[rgb]{0,0,0}M_{\mathtt{o}}},{\color[rgb]{0,0,0}\boldsymbol{\mathrm{\eta}}}, respectively. The Lorentz transformation $L$ such that the two conditions (5) and (6) hold, is found in two steps. First, one expresses {\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}} in terms of {\color[rgb]{0,0,0}\boldsymbol{\Lambda}_{\mathrm{s}}},{\color[rgb]{0,0,0}e},{\color[rgb]{0,0,0}m_{\mathtt{o}}}; second, one determines {\color[rgb]{0,0,0}\boldsymbol{\Lambda}_{\mathrm{s}}} by solving one of the constraints.111If an evolution equation for {\color[rgb]{0,0,0}\boldsymbol{\Lambda}_{\mathrm{s}}} is known, as in the case of spherical symmetry, one can also specify the value of {\color[rgb]{0,0,0}\boldsymbol{\Lambda}_{\mathrm{s}}} on the initial hypersurface and evolve it in time, without the need to solve it from one of the constraints. Here we focus on the first step, i.e., the determination of {\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}}. The expression of {\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}} is found to be [35]

[TABLE]

which is the polar decomposition of {\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}}^{-1}. Since the polar decomposition always exists, it is always possible to write {\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}} in these terms. In addition, {\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}}^{-1} being clearly invertible, it follows that {\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}}^{-1}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}}^{\mathtt{{\scriptscriptstyle\top}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}} is strictly positive definite [36, Sec. 2.5]. This means that the computation of {\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}} reduces to the computation of the square root of a $3\times 3$ symmetric positive definite matrix. This is the first step in computing the bimetric $3+1$ and BSSN decomposition, and all the successive steps depend on it.

The recasting of the standard $N+1$ equations in the covariant BSSN formalism, introduces new dynamical variables [17, 18, 19, 22, 23]. Following the notation in [23], they are the conformal metrics,

[TABLE]

the conformal extrinsic curvatures and traces,

[TABLE]

and the conformal connections,

[TABLE]

where the background connections {\color[rgb]{0,0,0}{\widebar{\Gamma}}_{\mathtt{B}}}^{i}_{jk},{\color[rgb]{0,0,0}{\widehat{\Gamma}}_{\mathtt{B}}}^{i}_{jk} are arbitrary but time-independent222The assumption of time-independency is due to the fact that the arbitrary connections do not have to fulfill any evolution equation. In bimetric relativity, there is the possibility to set {\color[rgb]{0,0,0}\chi} as the background geometry for {\color[rgb]{0,0,0}\gamma} and {\color[rgb]{0,0,0}\varphi}, hence the assumption of time-independency can be relaxed. See [24] for more details., possibly arising as the compatible connections for the background metrics. We also define the conformal mean metric {\color[rgb]{0,0,0}{\accentset{\circ}{{\chi}}}},

[TABLE]

which is not a dynamical variable. For the bimetric covariant BSSN equations, we refer the reader to [23].

We are finally in the position to state what are the primary variables needed to be specified in order to compute the bimetric covariant BSSN decomposition,

[TABLE]

where {\color[rgb]{0,0,0}{\widebar{e}}},{\color[rgb]{0,0,0}{\widehat{m}}} are the conformal vielbeins,

[TABLE]

These are the variables used to make an ansatz, in the covariant BSSN formalism. Once they are known, the entire decomposition defined in [16, 23] can be computed.333Note that choosing an ansatz for {\color[rgb]{0,0,0}\boldsymbol{\mathrm{p}}}^{\textbf{a}} and {\color[rgb]{0,0,0}q}^{i} is equivalent to choosing an ansatz for {\color[rgb]{0,0,0}\beta}^{i} and {\color[rgb]{0,0,0}\mkern 1.0mu\widetilde{\mkern-1.0mu\beta}}^{i}.

2 The abstract equations

2.1 Basic definitions in $\mathtt{xAct}$

First, $\mathtt{bimEX}$ loads the $\mathtt{xAct}$ bundle, after which the three dimensional spacelike hypersurface $\Sigma_{t}$ with abstract indices $\mathtt{\{i,j,k,q,r,s\}}$ is defined,

The time coordinate is not defined on this manifold, but we need our objects to depend on it. For this reason, $\mathtt{bimEX}$ defines a parameter $\mathtt{t}$,

All the abstract tensors defined in $\mathtt{bimEX}$ depend on this parameter. The time derivatives are performed using the $\mathtt{xAct}$ built-in function $\mathtt{ParamD[\langle parameter\rangle][\langle expression\rangle]}$. Since the evolution equations in the covariant BSSN formalism contain the differential operators \partial_{t}-{}_{\color[rgb]{0,0,0}\beta} and \partial_{t}-{}_{\color[rgb]{0,0,0}\mkern 1.0mu\widetilde{\mkern-1.0mu\beta}}, we also define the parameters $\mathtt{g\digamma}$ and $\mathtt{f\digamma}$. A third parameter $\mathtt{h\digamma}$ is also defined to represent the operator \partial_{t}-{}_{\color[rgb]{0,0,0}q}. All the defined abstract tensors depend on these three parameters as well. In this way, $\mathtt{bimEX}$ provides a way to represent the differential operators \partial_{t}-{}_{\color[rgb]{0,0,0}\beta},\partial_{t}-{}_{\color[rgb]{0,0,0}\mkern 1.0mu\widetilde{\mkern-1.0mu\beta}},\partial_{t}-{}_{\color[rgb]{0,0,0}q} in a simple way. Their explicit form must be provided by the user, when needed. These $\mathtt{xAct}$ representations should be considered as placeholders, but their utility in the abstract computations relies on the fact that the properties of a derivative operator, e.g., the Leibniz rule, are automatically implemented for them by $\mathtt{xAct}$.

Next, $\mathtt{bimEX}$ defines all the abstract tensors representing the geometrical objects defined in [16, 23]. We note that $\mathtt{bimEX}$ defines six metrics; the three spatial metrics {\color[rgb]{0,0,0}\gamma},{\color[rgb]{0,0,0}\varphi},{\color[rgb]{0,0,0}\chi} and their conformally related metrics {\color[rgb]{0,0,0}{\widebar{\gamma}}},{\color[rgb]{0,0,0}{\widehat{\varphi}}},{\color[rgb]{0,0,0}{\accentset{\circ}{{\chi}}}}. An important point to stress is that the abstract tensors representing the vielbeins of the metrics in $\mathtt{xAct}$, have only spatial indices, rather than a spatial and a Lorentz index. This can potentially lead to confusion, since one cannot distinguish between the vielbein {\color[rgb]{0,0,0}e}^{\textbf{a}}{}_{i} and its inverse {\color[rgb]{0,0,0}e}^{i}{}_{\textbf{a}}, with a Lorentz index and $i$ spatial index. For this reason, the vielbein and the inverse vielbein are two different abstract tensors in $\mathtt{bimEX}$,

This solution is effective and simple enough that we did not encounter any need to build more complex structures to describe the Lorentz frame, during the work that led to the results in [23].444We could introduce another manifold with another set of indices representing the Lorentz frame, and then connect the two manifolds with a suitable map. However, as we said, there was no need to do that, which would perhaps result in a more mathematically rigorous, but less intuitive structure of the code.

2.2 Raising and lowering indices

In bimetric relativity there are (at least) two metric sectors, therefore raising and lowering indices must be done accordingly. $\mathtt{xAct}$ allows the definition of one active metric only, which automatically raises and lowers all the indices. The active metric in $\mathtt{bimEX}$, at the present version, is {\color[rgb]{0,0,0}\gamma}, i.e., the physical spatial part of the metric {\color[rgb]{0,0,0}g}. This is clearly a problem in bimetric relativity, since we do not want one metric to raise or lower the indices in the other metric sectors.

This means that we need to take care of contractions explicitly. This is done by defining the tensors with some canonical indices and never write them with the indices in different positions. Suppose the user defines the extrinsic curvature in the {\color[rgb]{0,0,0}f}-sector {\color[rgb]{0,0,0}\mkern 1.0mu\widetilde{\mkern-1.0muK}}_{ij}, and wants to raise the index $i$. The user should not write it in upper position directly, as in {\color[rgb]{0,0,0}\mkern 1.0mu\widetilde{\mkern-1.0muK}}^{i}{}_{j}, because this expression in $\mathtt{xAct}$ would be equivalent to {\color[rgb]{0,0,0}\gamma}^{ik}{\color[rgb]{0,0,0}\mkern 1.0mu\widetilde{\mkern-1.0muK}}_{kj}. Rather, the user should explicitly write {\color[rgb]{0,0,0}\varphi}^{ik}{\color[rgb]{0,0,0}\mkern 1.0mu\widetilde{\mkern-1.0muK}}_{kj}. Also, $\mathtt{bimEX}$ includes a function, inspired by the built-in $\mathtt{xAct}$ function $\mathtt{Simplification}$, which simplifies the abstract expressions without raising and lowering indices. It is called $\mathtt{SafeSimplification}$ and it is defined as,

The explicit writing of every contraction and the use of $\mathtt{SafeSimplification}$ results in unambiguous computations.

2.3 The definitions of the covariant BSSN equations

The next step is the definition of the covariant BSSN equations as abstract tensors in $\mathtt{xAct}$. We take the Hamiltonian constraint in the {\color[rgb]{0,0,0}g}-sector as an example, but the same procedure is implemented for all the evolution and constraint equations in [23]. At the present stage, the package only includes the covariant BSSN equations, not the standard $N+1$ ones. However, $\mathtt{bimEX}$ easily allows to implement them, since it is enough to mimic the procedure described below. Hence, in order to be compatible with later versions that could include also other sets of equations, we prefix the objects referring to the covariant BSSN with the string \mathtt{cBSSN\}$. We can see this explicitly in the definition below

This tensor is now only a placeholder for the Hamiltonian constraint in the {\color[rgb]{0,0,0}g}-sector. We would like the component of this scalar in a given chart to be the Hamiltonian constraint in that chart. On the other hand, we also want to manipulate the abstract Hamiltonian constraint. Therefore, we define an $\mathtt{IndexRule}$ to replace \mathtt{cBSSN\gHamiltonianConstraint[]}$with its explicit tensorial expression in the covariant BSSN. The rule is called$\mathtt{Instantiate$gHC}$.

The command in Listing 6 prints the abstract Hamiltonian constraint in the covariant BSSN formulation. There are also other $\mathtt{IndexRules}$: $\mathtt{InstantiateConstraints}$, $\mathtt{InstantiateEvolution}$ and $\mathtt{InstantiatePDE}$ do the same job for the constraint equations only, for the evolution equations only and for both of them.

Once the equations are instantiated, one can use the $\mathtt{xAct}$ functions to manipulate them. All the bimetric interactions and sources are also defined as abstract tensors by $\mathtt{bimEX}$.

3 The concrete equations

3.1 The computation of the square root in (7a)

The computation of ({\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}}^{-1}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}}^{\mathtt{{\scriptscriptstyle\top}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}})^{1/2} in (7a) needs to be performed with care, since we are dealing with symbolic manipulation. Hence, the computation can be inefficient or result in a very complicated expression for ({\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}}^{-1}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}}^{\mathtt{{\scriptscriptstyle\top}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}})^{1/2}, depending on the chosen ansatz and on the algorithm used to compute it.

We have implemented three algorithms to compute this square root. Depending on the ansatz, one can perform better than the others. Which method is the best has to be inspected case by case. The user can choose in a quite simple way which method to use, as we will explain in the next subsection. The first method, which is the default and simplest one, is based on the built-in Mathematica function $\mathtt{MatrixPower}$. The implementation is the following,

The second method is an adapted version of the one reported in [37], and consists in computing the polar decomposition of {\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}}^{-1}, using the built-in Mathematica function $\mathtt{SingularValueDecomposition}$,

Note that, in these first two methods, we explicitly use the fact that we are dealing with real quantities. The third method is the implementation of the algorithm presented in [38]. This algorithm is very efficient when dealing with numbers. In the case of symbolic manipulation, it is not guaranteed that it can perform better than the other two. Its implementation is,

This function has the option $\mathtt{ApplyFunction}$, which is discussed in the next subsection.

To reproduce the results in [16] and obtain the results in [23, 24], we used the algorithm in Listing 7, since the chosen ansatz imply that {\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}}^{-1}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}}^{\mathtt{{\scriptscriptstyle\top}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}} is diagonal, and the function $\mathtt{MatrixPower}$ is efficient enough. In this case, the algorithm in Listing 9 gives a much more complicated expression, which has to be simplified in a second step to give the simpler expression obtained with the algorithm in Listing 7. Also, in this case the efficiency of the algorithm in Listing 8 is comparable with the one in Listing 7.

We stress the fact that, since we are dealing with symbolic manipulation, it is important to recognize repetitive patterns appearing in the equations. For this reason we defined the shifted elementary symmetric polynomials in [6], which considerably simplify the computations and increase the efficiency of the code. We verified this explicitly in the case of spherical symmetry. Our approach was to blindly compute a part of the decomposition as the first step, and then to recognize repetitive patterns which can be represented in Mathematica as single symbols. The substitution of the repetitive patterns with single symbols speeds up the symbolic manipulation tremendously.

3.2 Computation of the components in a given chart

Here we describe the functions that compute the bimetric BSSN decomposition in a given chart. First, one has to define a chart. This is made easier by the function $\mathtt{DefChartScalars}$, which is based on the built-in $\mathtt{xAct}$ function $\mathtt{DefChart}$ and has its same $\mathtt{Options}$. It defines a chart on the spacelike hypersurface, assigns components to the identity operator in that chart, and sets the independence of the parameters $\mathtt{t,g\digamma,f\digamma,h\digamma}$ from the coordinates in the chart, i.e.,

In addition, $\mathtt{DefChartScalars}$ sets the values of the variables $\mathtt{FirstChart}$ and $\mathtt{DefaultChart}$ to the name of the chart, henceforth considered as the default chart for all the functions. The variable $\mathtt{DefaultChart}$ can be modified, and in the case of multiple charts it is reset to the last defined chart. The value of $\mathtt{FirstChart}$ is $\mathtt{Protected}$ instead. In this way, $\mathtt{bimEX}$ can handle multiple charts simultaneously. In case the user is using only one chart, there is no need to specify it when dealing with the components of the tensors. $\mathtt{DefChartScalars}$ has the option $\mathtt{ChartAssumptions}$, which allows to specify a list containing the assumptions on the coordinates. As an explicit example, this commands defines the spherical polar chart on the spacelike hypersurface,

Once a chart is defined, the ansatz should be chosen. The ansatz is the input to the function $\mathtt{ComputeBSSNDecomposition}$. We must choose an ansatz for the primary variables in (12), plus a background metric for {\color[rgb]{0,0,0}{\accentset{\circ}{{\chi}}}}_{ij}:

The conformal factors {\color[rgb]{0,0,0}\phi}, {\color[rgb]{0,0,0}\psi} 2. 2.

The conformal vielbeins {\color[rgb]{0,0,0}{\widebar{e}}}^{\textbf{a}}{}_{i}, {\color[rgb]{0,0,0}{\widehat{m}}_{\mathrm{o}}}^{\textbf{a}}{}_{i} 3. 3.

The vector {\color[rgb]{0,0,0}\boldsymbol{\mathrm{p}}}^{\textbf{a}} which completely determines {\color[rgb]{0,0,0}\boldsymbol{\Lambda}_{\mathrm{s}}}^{\textbf{a}}{}_{\textbf{b}} and, together with {\color[rgb]{0,0,0}{\widebar{e}}}^{\textbf{a}}{}_{i}, {\color[rgb]{0,0,0}{\widehat{m}}_{\mathrm{o}}}^{\textbf{a}}{}_{i}, {\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}}^{\textbf{a}}{}_{\textbf{b}} (hence $L$) 4. 4.

The shift vector {\color[rgb]{0,0,0}q}^{i} of the geometric mean metric {\color[rgb]{0,0,0}h} 5. 5.

The conformal extrinsic curvatures {\color[rgb]{0,0,0}{\widebar{A}}}^{i}{}_{j}, {\color[rgb]{0,0,0}{\widehat{A}}}^{i}{}_{j} 6. 6.

The conformal connections {\color[rgb]{0,0,0}{\widebar{\Lambda}}}^{i}, {\color[rgb]{0,0,0}{\widehat{\Lambda}}}^{i} 7. 7.

The three background metrics for {\color[rgb]{0,0,0}{\widebar{\gamma}}}_{ij}, {\color[rgb]{0,0,0}{\widehat{\varphi}}}_{ij}, {\color[rgb]{0,0,0}{\accentset{\circ}{{\chi}}}}_{ij}

$\mathtt{ComputeBSSNDecomposition}$ uses background connections arising from background metrics. The background metric for {\color[rgb]{0,0,0}{\accentset{\circ}{{\chi}}}}_{ij} is an input variable needed to compute the dynamics of the geometric mean in the covariant BSSN formulation, included for forward compatibility.

After we choose the ansatz for these variables, we give it to $\mathtt{ComputeBSSNDecomposition}$ as input. This function computes all the components of the bimetric interactions and sources, storing them both as the components of the abstract tensors in $\mathtt{xAct}$, and as plain lists in $\mathtt{Mathematica}$. Hence, the user will be able to use the computed components within $\mathtt{bimEX}$ or not, as we will discuss in Sec. 3.5. Since the components are saved as plain lists, the user can extend $\mathtt{bimEX}$ to handle them in $\mathtt{CTensor}$ as well, in a straightforward way. All geometrical quantities—Ricci tensor, Christoffel symbols, etc.—in the metric sectors are also computed by $\mathtt{ComputeBSSNDecomposition}$, making use of the $\mathtt{xAct}$ built-in function $\mathtt{MetricCompute}$. $\mathtt{ComputeBSSNDecomposition}$ also performs many internal checks to verify some of the relations involving the bimetric interactions. If any of these checks fails, the evaluation aborts and an informative message about the error is printed.

The function $\mathtt{ComputeBSSNDecomposition}$ has five options that can make the computations more efficient, and allow the user to decide what to compute. The options are,

$\mathtt{ChartName}$. The default value for this option is $\mathtt{DefaultChart}$. If the user needs to use multiple charts, they should set $\mathtt{ChartName}$ to the name of the chart they want to compute the decomposition in. 2. 2.

$\mathtt{ComputeGeometryOf}$. This option should be set equal to a list containing the names of the metrics whose geometrical quantities we want to compute. The default value is \left\{{\color[rgb]{0,0,0}\gamma},{\color[rgb]{0,0,0}\varphi}\right\}, which means that the geometries of {\color[rgb]{0,0,0}\gamma},{\color[rgb]{0,0,0}{\widebar{\gamma}}} and {\color[rgb]{0,0,0}\varphi},{\color[rgb]{0,0,0}{\widehat{\varphi}}} will be computed, but the geometry of {\color[rgb]{0,0,0}\chi},{\color[rgb]{0,0,0}{\accentset{\circ}{{\chi}}}} will not, since it slows down the execution. The user can decide to compute the geometries of any of the metrics. 3. 3.

$\mathtt{ApplyFunction}$. This option allows the user to specify a function to be applied to the results of the computations performed in $\mathtt{ComputeBSSNDecomposition}$. The default value is the identity, i.e., nothing is applied to the computed expressions. The chosen function has a direct impact on the efficiency of $\mathtt{ComputeBSSNDecomposition}$. It is desirable to identify recurrent patterns in the equations to make the computation faster and more useful. The user can define a function $f$ that recognizes these patterns in the chosen ansatz, and set $f$ as the value of $\mathtt{ApplyFunction}$. This increased the efficiency tremendously for our computations in spherical symmetry. 4. 4.

$\mathtt{SqrtAlgorithm}$. This option allows the user to choose the algorithm to compute the square root matrix ({\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}}^{-1}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}}^{\mathtt{{\scriptscriptstyle\top}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}})^{1/2}. It can be set to three strings, $\mathtt{``MatSqrt"}$, $\mathtt{``PolDec"}$, $\mathtt{``PosDefSqrt"}$. They refer to the three algorithm discussed in Sec. 3.1, respectively to the algorithms in Listing 7, Listing 8 and Listing 9. 5. 5.

$\mathtt{IndVariables}$. This option allows the user to specify a list containing the independent variables to be set as arguments for the primary fields. The default value is the string $\mathtt{``AutoDetect"}$, which makes $\mathtt{ComputeBSSNDecomposition}$ to inspect the ansatz and find the independent variables. The following shortcut is then defined, for a generic function $f$,

[TABLE]

An example call of $\mathtt{ComputeBSSNDecomposition}$ is

Note that the order of the input variables in $\mathtt{ComputeBSSNDecomposition}$ is important, and that the two conformal vielbeins {\color[rgb]{0,0,0}{\widebar{e}}},{\color[rgb]{0,0,0}{\widehat{m}}_{\mathrm{o}}} have to be upper triangular.

We remark the following limitation. The $\mathtt{xAct}$ function $\mathtt{MetricCompute}$, used by $\mathtt{ComputeBSSNDecomposition}$ to compute the components of the geometrical objects associated with the metrics, does not store the computed components in all the charts. For example, suppose the user uses two charts $C_{1}$ and $C_{2}$. The user executes $\mathtt{MetricCompute}$ in the chart $C_{1}$ first, and immediately afterwards in the chart $C_{2}$. The second execution in the chart $C_{2}$ overwrites the previously computed components in the chart $C_{1}$. This has to be taken into account when using multiple charts in $\mathtt{bimEX}$.

3.3 How to access to the computed decomposition besides $\mathtt{xAct}$

There is a standard notation for the lists storing the components of the tensors computed by $\mathtt{ComputeBSSNDecomposition}$, which helps the user to remember their names. For definiteness, let’s consider the components of the metric {\color[rgb]{0,0,0}{\widehat{\varphi}}}_{ij} and the conformal extrinsic curvature {\color[rgb]{0,0,0}{\widebar{A}}}_{ij}. They are stored in the lists defined as,

where the “c” in $\mathtt{\varphi c}$ stands for “conformal” metric. If the user uses one chart only, there is no need to specify the name of the chart. The symbols in the names tell us what type of components are we looking at. The symbol $\blacktriangledown\blacktriangledown$ means that the list \mathtt{gA}$$\blacktriangledown\blacktriangledown$$\mathtt{[]} stores the components of {\color[rgb]{0,0,0}{\widebar{A}}}_{ij} with both lower indices. The components of {\color[rgb]{0,0,0}{\widebar{A}}}^{i}{}_{j} are stored in the list \mathtt{gA}$$\blacktriangleup\blacktriangledown$$\mathtt{[]}, and so on. This holds for all the rank 2 tensors defined in $\mathtt{bimEX}$. However, there are exceptions to this rule. One exception, as we see in Listing 13, concerns the metrics. The list of components of {\color[rgb]{0,0,0}{\widehat{\varphi}}} contains the symbol $\blacksquare$ in its name. The component of the inverse conformal metric are stored in the list named \overset{-1}{\varphi c}\blacksquare$$\mathtt{[]}. The mixed components of the metrics are those of the identity, and are stored in lists named with the $\blacktriangleup\blacktriangledown$ and $\blacktriangledown\blacktriangleup$ symbols. Another exception is given by the lists containing the components of Lorentz linear operator, e.g., {\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}},{\color[rgb]{0,0,0}\boldsymbol{\Lambda}_{\mathrm{s}}}, and the vielbeins. These lists have a rectangle in their name, given by the Mathematica code $\mathtt{\backslash[FilledRectangle]}$.

Now the user can access the components stored in plain lists. Suppose the user computes the decomposition and wants to see it. The function $\mathtt{PrintDecomposition}$ prints the components of the main variables in the decomposition. The user can call it as follows,

$\mathtt{PrintDecomposition}$ has no arguments and two options, $\mathtt{Of}$ and $\mathtt{ApplyFunction}$. The value of $\mathtt{Of}$ is a list containing the names of the sectors which we want to print the decomposition of. In Listing 14, the command will print the quantities associated with all three metric sectors.

3.4 From abstract to concrete equations in $\mathtt{xAct}$

Once $\mathtt{ComputeBSSNDecomposition}$ is executed, all the bimetric interactions, sources and the geometrical quantities of the metrics have assigned components in $\mathtt{xAct}$. The function $\mathtt{ToConcrete}$ takes an abstract equation, and writes it in components in the given chart. Its efficiency depends on the level of instantiation that the user can control via three boolean options, listed below and bespoke for the covariant BSSN equations,

$\mathtt{ConcreteSources}$. If $\mathtt{True}$, it will instantiate the bimetric sources. Setting this options to $\mathtt{False}$ can be very useful for the readability of the final equations, and their manipulation in Mathematica, since there are less terms to manipulate. These terms do not contain derivatives of the dynamical fields. 2. 2.

$\mathtt{ConcreteRicci}$. If $\mathtt{True}$, it will instantiate the Ricci tensors in

[TABLE]

appearing in the evolution equations for the conformal extrinsic curvatures {\color[rgb]{0,0,0}{\widebar{A}}}^{i}{}_{j},{\color[rgb]{0,0,0}{\widehat{A}}}^{i}{}_{j} [22, 23]. Turning this option to $\mathtt{False}$ speeds up the computation considerably, and helps to see the structure of the equations, which is needed to optimize them before proceeding to the numerical integration. 3. 3.

$\mathtt{ConcreteShift}$. If $\mathtt{True}$, it will instantiate the shift vector of {\color[rgb]{0,0,0}g} and {\color[rgb]{0,0,0}f} according to (5). This has the same advantages of $\mathtt{ConcreteSources}$.

The default value for the three options is $\mathtt{False}$. $\mathtt{ChartName}$ is also an option for $\mathtt{ToConcrete}$, and works in the same way. Setting the three options to $\mathtt{False}$ results in a very efficient instantiation algorithm.

Three more functions are defined to be able to instantiate the bimetric sources, the Ricci tensors and the shifts independently from $\mathtt{ToConcrete}$. They are called $\mathtt{ToConcreteSources}$, $\mathtt{ToConcreteRicci}$ and $\mathtt{ToConcreteShift}$. They can be applied to a concrete equation to instantiate the desired components.

Note that, even if we do not instantiate the components of any of these quantities by setting the three options of $\mathtt{ComputeBSSNDecomposition}$ to $\mathtt{False}$, only their non-zero components will be kept uninstantiated. The zero components will be set to zero not to keep unnecessary terms and simplify the output. The recognition of the zero components is done by $\mathtt{ComputeBSSNDecomposition}$. For definiteness, let’s consider the Ricci tensor as an example. $\mathtt{ComputeBSSNDecomposition}$ computes explicitly all the components of the Ricci tensor, and recognizes the zero ones. Then, it will set to zero all the respective components of the abstract Ricci tensor in $\mathtt{xAct}$. The non-zero components of the abstract Ricci tensor in $\mathtt{xAct}$ will be set equal to some (appropriately named) scalar functions, acting as placeholders. The functions $\mathtt{ToConcrete}$, $\mathtt{ToConcreteSources}$, $\mathtt{ToConcreteRicci}$ and $\mathtt{ToConcreteShift}$ only replace the placeholders scalar functions with the actual components computed only once by $\mathtt{ComputeBSSNDecomposition}$. This results in a fast execution for all the $\mathtt{ToConcrete}$ commands.

The function $\mathtt{ToConcrete}$ is a switch function which depends on the values of the three options of instantiation. It calls different functions depending on the values of its options. Suppose we choose to set the three options to False. Then, $\mathtt{ToConcrete}$ calls only the function \mathtt{ToConcrete\Basic}$. The function$\mathtt{ToConcrete$Basic}$is based on the built-in$\mathtt{xAct}$functions—$\mathtt{ToBasis}$,$\mathtt{TraceBasisDummy}$,$\mathtt{ComponentArray}$,$\mathtt{ToValues}$,$\mathtt{InChart}$—but it implements other algorithms as well, since those functions alone do not suffice to write everything in components.555As an example, the abstract coordinate derivative$\mathtt{PD}$, used by$\mathtt{xAct}$in the abstract equations, has to be replaced by hand by the partial derivative in the given basis$\mathtt{PDOfBasis}$. The Christoffel symbols of the metrics need to be taken into account separately as well. Essentially,$\mathtt{ToConcrete$Basic}$instantiates everything but the Ricci tensors, the bimetric sources and the shift vectors. Suppose that we want to instantiate an abstract equation including the shifts. We set the option$\mathtt{ConcreteShift}$in$\mathtt{ToConcrete}$to$\mathtt{True}$.$\mathtt{ToConcrete}$then calls$\mathtt{ToConcrete$Basic}$first, and$\mathtt{ToConcreteShift}$ later. The same applies for all the other options’ combinations.

As an example, the code to instantiate the Hamiltonian constraint in the {\color[rgb]{0,0,0}g}-sector in the $\mathtt{DefaultChart}$ is

3.5 Exportation of the bimetric decomposition and equations

At this point, we know how to compute the BSSN decomposition and how to write the abstract equations in components, with the desired level of instantiation. In the case of spherical symmetry, on a HP Z240 SFF Workstation, with an Intel(R) Core(TM) i7-6700 CPU 3.40GHz and 16.0 GB of RAM memory, running Windows 8 64-bit, $\mathtt{Mathematica}$ 11.0 and $\mathtt{xAct}$ 1.1.3, the execution of $\mathtt{ComputeBSSNDecomposition}$—computing the geometrical quantities for all the six metrics—takes about 2 minutes, and the instantiation and simplification of all the bimetric covariant BSSN equations takes about 10 minutes in total. Once the decomposition is computed and the equations are instantiated and simplified to the desired level, they can be exported into an .m file. This removes the need to make the same computations each time that the user needs to use, e.g., the constraint equations in spherical symmetry.

We already said that $\mathtt{ComputeBSSNDecomposition}$ saves the computed components of the tensors both as lists and as components stored in the $\mathtt{xAct}$ tensors. This is because we do not want to be necessarily constrained to the $\mathtt{xAct}$ bundle. The user should be able to export the decomposition in a file which does not have any memory of $\mathtt{xAct}$. In this case, the user could open a Mathematica notebook, load the .m file where the decomposition is saved in the form of plain lists, and use these lists to perform the desired computations using only the standard Mathematica functions. On the other hand, the user also needs to be able to export the decomposition within $\mathtt{xAct}$, i.e., to export all the definition and rules defined by $\mathtt{bimEX}$ in $\mathtt{xAct}$. In this second case, the user should be able open a Mathematica notebook, load $\mathtt{bimEX}$ and the .m file where the decomposition in $\mathtt{xAct}$ is exported, and use $\mathtt{bimEX}$ and the $\mathtt{xAct}$ bundle to perform the desired computations.

Before introducing the two functions which export the bimetric decomposition and equations in the two different ways outlined above, we need to hightlight another feature of $\mathtt{ComputeBSSNDecomposition}$. This function does not assign components to the abstract tensors representing the covariant BSSN constraint and evolution equations. For clarity, it does not assign components to \mathtt{cBSSN\gHamiltonianConstraint[]}in Listing [5](#LST5). This is the case for two reasons. First,\mathtt{ComputeBSSNDecomposition}$performs all computations$\mathtt{ToConcrete}$needs in order to write the equations in components. Therefore, for$\mathtt{ComputeBSSNDecomposition}$to assign components to the abstract equations,$\mathtt{ToConcrete}$would be needed to be part of it. On the other hand,$\mathtt{ToConcrete}$needs$\mathtt{ComputeBSSNDecomposition}$to be executed first, so this type of implementation would require a basic restructuring of the code. Second, after the execution of$\mathtt{ToConcrete}$, the equations are in concrete form, but they are not necessarily written in the simplest possible form. Hence, one may want to simplify or manipulate them, and then set the simplified equations as the components of the abstract tensors (which are the ones being exported).

The components are assigned to the tensors representing the covariant BSSN equations by the function $\mathtt{AssignComponents}$, which has $\mathtt{ChartName}$ and $\mathtt{ApplyFunction}$ as options. It takes as arguments two ordered lists, the first one containing the abstract tensors to which we want to assign the components, and the second one containing the components. It applies the value of $\mathtt{ApplyFunction}$ to the list of components, and set them as the components of the abstract tensors in the chart $\mathtt{ChartName}$. Again, the default value of $\mathtt{ChartName}$ is $\mathtt{DefaultChart}$, and the default value of $\mathtt{ApplyFunction}$ is the identity. At the present stage, this function works for tensors up to rank 2. After having assigned the components, the user can access the concrete equations with the following code, as an alternative to the one in Listing 15,

We are now ready to export the decomposition to an .m file. There are two functions to do that, $\mathtt{ExportDecomposition}$ and \mathtt{ExportDecomposition\xAct}$. The first one exports all the components stored into lists, in an .m file which is independent from the$\mathtt{xAct}$bundle and from$\mathtt{bimEX}$as well. Its arguments are the fields in the chosen ansatz, without the conformal factors. It has four options,$\mathtt{ChartName}$,$\mathtt{ApplyFunction}$,$\mathtt{NameFile}$and$\mathtt{OtherQuantities}$. The first two work in the same way as for the other functions, and$\mathtt{NameFile}$allows to specify the name of the exported .m file. The default value of$\mathtt{NameFile}$is$\mathtt{``Bimetric_decomposition_\langle time;and;date\rangle.m"}$. The option$\mathtt{OtherQuantities}$allows the user to add other quantities to be exported. It should be set to a list of two lists. The first list includes the names of the variables to be exported, and the second list the expressions to be saved into these variables. This makes it possible to export the concrete equations, which are not exported by$\mathtt{ExportDecomposition}$by default. If the option$\mathtt{OtherQuantitites}$is used, it is necessary to execute the command$\mathtt{MapThread[Set,{$$names,$$values}]}$after loading the .m file. In addition, the expressions included in the value of$\mathtt{OtherQuantities}$and specified by the user, must be totally instantiated, i.e., they should not contain any reference to xAct.$\mathtt{ExportDecomposition}$locally clears the$\mathtt{UpValues}$of the variables and uses$\mathtt{DumpSave}$to export them, since this function does not export the$\mathtt{FullDefinition}$of a variable, contrary to$\mathtt{Save}. The variables exported by default are named with the same convention described in [subsection 3.3](#S3.SS3), but the last two characters in the names, namely \mathtt{[]}$, are replaced by$\mathtt{$}$. For clarity,$\mathtt{gA}$\blacktriangleup\blacktriangledown$\mathtt{[]}$becomes$\mathtt{gA}\blacktriangleup\blacktriangledown\ in the exported .m file.

\mathtt{ExportDecomposition\xAct}$takes no arguments and exports all the information about the tensors defined in$\mathtt{xAct}$and their components into an .m file, using$\mathtt{Save}$. Before loading the .m file, the user needs to load$\mathtt{bimEX}$, since it defines the spacelike hypersurface and all the abstract tensors. The .m file exported by$\mathtt{ExportDecomposition$xAct}$will then *overwrite* all the information stored in it to the basic information defined by$\mathtt{bimEX}$. One could also think about this as *adding* the information concerning the components of all the tensors within$\mathtt{xAct}$, to the basic information contained in$\mathtt{bimEX}$.$\mathtt{ExportDecomposition$xAct}$only has the option$\mathtt{NameFile}$, whose default value is$\mathtt{``Bimetric_decomposition_xAct_\langle time;and;date\rangle.m"}$.

4 Conclusions

We presented $\mathtt{bimEX}$, a package to perform bimetric exact computations in the $3+1$ formalism in $\mathtt{Mathematica}$. It makes extensive use of the $\mathtt{xAct}$ bundle and can handle both abstract and concrete—i.e., with components—computations.

The package was written during the work to obtain the results in [23, 24]. Therefore, at the present stage, it includes the covariant BSSN equations written as abstract tensors, together with the bimetric $3+1$ and BSSN decompositions [16, 23]. The user can manipulate these equations abstractly by using the $\mathtt{xAct}$ built-in functions. However, there is a key factor that must be reminded: the user must have complete control on how the indices are raised and lowered, since the objects in one metric sector cannot be contracted with the other metric. This constitute a problem in $\mathtt{xAct}$, because it does not allow to define only frozen metrics, where by frozen metrics we mean metrics that do not raise and lower indices automatically. Our solution to this problem is to use the $\mathtt{xAct}$ function $\mathtt{ToCanonical}$ always with the option $\mathtt{UseMetricOnVBundle\rightarrow None}$; we defined the function $\mathtt{SafeSimplification}$ which automatically does this. The function $\mathtt{SafeSimplification}$ never raises and lowers indices during the simplifications. In addition, one has to write contractions explicitly.

The package allows to choose an ansatz on the primary variables of the theory, listed in (12), to give it as input to the function $\mathtt{ComputeBSSNDecomposition}$, and to obtain the bimetric decomposition as output. After the bimetric decomposition is computed, the function $\mathtt{ToConcrete}$, applied to an abstract quantity, writes it into components. This function allows the user to choose the level of instantiation of the output, in order to increase efficiency and readability of the equations written in components.

The package allows the user to print the decomposition via the function $\mathtt{PrintDecomposition}$, and export it, via the functions $\mathtt{ExportDecomposition}$ and \mathtt{ExportDecomposition\xAct}$. The first function exports the decomposition into a .m file which is independent from the$\mathtt{xAct}$bundle and$\mathtt{bimEX}$, since all the components of the tensors are saved into lists—i.e., matrices. The second function exports all the definitions of the$\mathtt{xAct}$tensors into an .m file, which requires$\mathtt{bimEX}$(hence$\mathtt{xAct}$) to be loaded first, in order to work properly. The user is then able to choose how to work with the decomposition.

In the computation of the bimetric decomposition, the square root of the positive definite $3\times 3$ matrix {\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}}^{-1}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}}^{\mathtt{{\scriptscriptstyle\top}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}}, which originates from the solution to the symmetrization condition (6), has to be computed. Three algorithms are implemented for this, one taken from [38], and the other two relying on Mathematica built-in functions. The user decides which algorithm to use, via an option of $\mathtt{ComputeBSSNDecomposition}$, which computes the square root. Since we deal with symbolic manipulation, there is no general rule about which method to use. The efficiency depends on the chosen ansatz and on the user-defined simplification functions. We suggest to define simplification functions which recognize repetitive patterns into the expressions and write them as single symbols.

We stress that this package can compute the covariant BSSN equations for any desired ansatz. It has been tested extensively for the spherically symmetric case, and should be useful to anyone working on numerical bimetric relativity.

Acknowledgments

We are thankful to Mikica Kocic for all the shared and fruitful Mathematica sessions, and to Edvard Mörtsell and Mikica Kocic for reading the paper and providing useful comments.