Covariant BSSN formulation in bimetric relativity

TL;DR

This paper develops a covariant BSSN formulation for bimetric relativity equations, highlighting differences from general relativity and discussing challenges for ensuring well-posedness in numerical simulations.

Contribution

It extends the covariant BSSN formulation to bimetric relativity, identifying key differences and challenges in achieving strong hyperbolicity for numerical solutions.

Findings

Shared features with general relativity's BSSN formulation

Differences depend on gauge choice affecting hyperbolic structure

Strong hyperbolicity not yet established for bimetric equations

Abstract

Numerical integration of the field equations in bimetric relativity is necessary to obtain solutions describing realistic systems. Thus, it is crucial to recast the equations as a well-posed problem. In general relativity, under certain assumptions, the covariant BSSN formulation is a strongly hyperbolic formulation of the Einstein equations, hence its Cauchy problem is well-posed. In this paper, we establish the covariant BSSN formulation of the bimetric field equations. It shares many features with the corresponding formulation in general relativity, but there are a few fundamental differences between them. Some of these differences depend on the gauge choice and alter the hyperbolic structure of the system of partial differential equations compared to general relativity. Accordingly, the strong hyperbolicity of the system cannot be claimed yet, under the same assumptions as in general relativity. In the paper, we stress the differences compared with general relativity and state the main issues that should be tackled next, to draw a road map towards numerical bimetric relativity.

Equations

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Delta}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\Delta}}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overset{\scriptscriptstyle{\#}}{\Delta}}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\Delta}}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\Delta}}}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\Delta}}}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\accentset{\circ}{\Delta}}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\accentset{\ast}{{\boldsymbol{\Delta}}}}}$$

$${\partial }_t{u}^i+A^k{}^i{}_j{\partial }_k{u}^j=S^i,$$

$$P^i{}_j:=n_k{A}^k{}^i{}_j.$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}G_{g}}^{\mu}{}_{\nu}$$

$$\dfrac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha}}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\alpha}}}=-\dfrac{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W_{f}}}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W_{g}}}\eqqcolon{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W}=\mbox{scalar field independent of the lapses and the shifts}.$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\gamma}}}_{ij}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\varphi}}}_{ij}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\phi}=\dfrac{1}{12}\log({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\Delta}),\qquad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\psi}=\dfrac{1}{12}\log({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\Delta}}).$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{A}}}_{ij}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{A}}}_{ij}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\Gamma}}}^{i}\coloneqq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\gamma}}}^{jk}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\Gamma}}}^{i}_{jk},\qquad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\Gamma}}}^{i}\coloneqq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\varphi}}}^{jk}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\Gamma}}}^{i}_{jk},$$

$$\displaystyle\partial_{t}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha}$$

$$\displaystyle\partial_{t}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta}^{i}$$

$$\displaystyle\partial_{t}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}B}^{i}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\gamma}={{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}e}}^{\mathsf{T}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}e},\qquad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varphi}={{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m}}^{\mathsf{T}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m},$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{R}}}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\Lambda}}=\begin{pmatrix}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\lambda}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}^{\mathsf{{\scriptscriptstyle T}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}}\\ {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}&{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\Lambda}_{\mathrm{s}}}\end{pmatrix}.$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta}\coloneqq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}q}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha}\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mathsfit{n}{\mkern 1.0mu}}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}q}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha}\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}e}^{-1}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\lambda}^{-1},\qquad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\beta}}\coloneqq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}q}-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\alpha}}\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\mathsfit{n}}{\mkern 1.0mu}}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}q}-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\alpha}}\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m}^{-1}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\lambda}^{-1}.$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\gamma}}}=\mathrm{e}^{-4{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\phi}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\gamma}\Longrightarrow{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{e}}}}^{\mathsf{T}}\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}}\,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{e}}}=\left(\mathrm{e}^{-2{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\phi}}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}e}}^{\mathsf{T}}\right){\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}}\left(\mathrm{e}^{-2{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\phi}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}e}\right),$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{e}}}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{e}}}^{\,-1}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\gamma}}}_{ij}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\varphi}}}_{ij}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{A}}}_{ij}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{A}}}_{ij}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{K}}}$$

$$\displaystyle{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\Lambda}}}^{i}\coloneqq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\gamma}^{jk}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bigtriangleup{\widebar{\Gamma}}}^{i}_{jk}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\gamma}^{jk}\Bigl{(}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\Gamma}}}^{i}_{jk}-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\Gamma}}_{\mathsf{B}}}^{i}_{jk}\Bigr{)},\qquad{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\Lambda}}}^{i}\coloneqq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varphi}^{jk}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\bigtriangleup{\widehat{\Gamma}}}^{i}_{jk}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varphi}^{jk}\Bigl{(}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\Gamma}}}^{i}_{jk}-{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\Gamma}}_{\mathsf{B}}}^{i}_{jk}\Bigr{)}.$$

$$\displaystyle\partial_{t}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha}$$

$$\displaystyle\partial_{t}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta}^{i}$$

$$\displaystyle\partial_{t}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}B}^{i}$$

Full text

1 Introduction

The Hassan–Rosen bimetric theory, or bimetric relativity (BR), is an extension of general relativity (GR) describing the nonlinear interaction of two metrics {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}_{\mu\nu} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}_{\mu\nu} Hassan:2011zd ; Hassan:2011ea . Its unambiguous definition and spacetime interpretation are established in HassanKocic2018 , and the complete Hamiltonian analysis of it is performed in HassanLundkvist2018 . The study of this theory is well-motivated. First, its cosmology is compatible with local gravity tests, and in particular with the recent observations of gravitational waves 10.1088/1361-6382/ab4f9b , since it describes the dynamics of both a massless and a massive spin-2 field, allowing for gravitational waves propagating with the speed of light. In addition, BR provides us with self-accelerating cosmological solutions Volkov:2011an ; vonStrauss:2011mq ; Comelli:2011zm ; Akrami:2012vf and a framework where the gravitational origin of dark matter can be studied Aoki:2014cla ; Blanchet:2015bia ; Enander:2015kda ; Babichev:2016bxi . This motivates further exploration of the theory and the work to obtain solutions to its field equations describing realistic physical systems. These are needed to compare the predictions of the theory against the observational data, possibly confirming or falsifying the theory.

At first, introducing a second metric in the same spacetime may sound exotic, but turns out to be necessary to make a spin-2 field massive.111There is no way to construct a dynamical mass term with only one metric, since {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}_{\mu\nu}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}^{\mu\nu}=4, the spacetime dimension. In GR, it is possible to couple a spin-2 field (the metric) to spin-0 and spin-1 fields. Hence, it is natural to try and see if and how a second spin-2 field can be coupled to the metric. Then, the question about the physical meaning of the second metric arises. In BR, independent matter sources are minimally coupled to only one of the two metrics PhysRevD.90.124042 ; Rham_2015 . This implies that a test particle coupled to the metric {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}_{\mu\nu} follows the geodesics determined by this metric, exactly as in GR. The interaction with {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}_{\mu\nu} is only experienced indirectly. The key difference is that the metric {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}_{\mu\nu} is not determined by the Einstein field equations (EFE), but rather by the bimetric field equations (BFE) which account for the interaction between the metrics. For more details about BR, we refer the reader to the review in Schmidt_May_2016 .

The field equations governing the dynamics of the two metrics have been solved analytically in some important cases (see Schmidt_May_2016 and references therein). In the majority of those cases, they are reduced to a set of ordinary differential equations. This is done, e.g., by imposing spherical symmetry and staticity—the radial coordinate is the sole independent variable—or spatial homogeneity and spherical symmetry—the time coordinate is the sole independent variable. Other exact solutions are equivalent to those of GR, which one can always obtain in BR in vacuum under the conditions established in doi:10.1142/S0218271814430020 ; doi:10.1063/1.5100027 , and in the presence of external matter sources under other conditions doi:10.1142/S0218271814430020 .

We are interested in solving the bimetric field equations to obtain more realistic solutions, e.g., the spherical gravitational collapse of matter. In this case, very few results can be obtained analytically (see, e.g., Hogas:2019cpg ; Hogas:2019ywm ) and the numerical integration of the BFE becomes indispensable. In addition, the BFE are now a system of partial differential equations (PDEs)—both the time and radial coordinate are independent variables—and, as such, it is desirable to recast them as a well-posed problem before integrating them. Well-posedness can be naively thought of as the concept of “stability against small perturbations” for PDEs courant1962methods . The concept of well-posedness for first-order systems of PDEs is introduced in subsection 1.1.

In GR, the EFE can be recasted in a variety of well-posed forms (see, e.g., (alcubierre2008introduction, , Ch. 5)), one of them being the Baumgarte, Shapiro, Shibata, Nakamura, Oohara and Kojima (BSSNOK or, more commonly, BSSN) formulation 10.1143/PTPS.90.1 ; PhysRevD.52.5428 ; PhysRevD.59.024007 ; PhysRevD.66.064002 ; PhysRevD.70.104004 (see also baumgarte2010numerical for a pedagogical intorduction), with its covariant generalization PhysRevD.79.104029 . One begins by rewriting the EFE as a Cauchy problem. This can be accomplished doing a $N+1$ decomposition, where the field equations split into a set of constraint equations not involving time derivatives, and a set of evolution equations involving first-order time derivatives. Then, in the free evolution scheme, one finds the appropriate initial data by solving the constraint equations, and evolves them by solving the evolution equations only. In Frittelli:1996nj it is proven that the constraints stay satisfied in the free evolution scheme, and the same holds in BR Kocic2019con . If one adopts the free evolution scheme, as we do, only the well-posedness of the Cauchy problem arising from the evolution equations is relevant.

Following the path outlined by numerical relativity, we would like to recast the BFE as a well-posed problem. The first step, i.e., the $N+1$ decomposition of the BFE, was established in Kocic:2018ddp . In this paper, we present the covariant BSSN (cBSSN) formulation of the BFE. The choice of the BSSN formulation is motivated by the fact that it is one of the most widely used in numerical relativity, allowing for stable long-term numerical evolution—see, e.g., subsection 11.4.5 in gourgoulhon20123+1 . A future research goal is to obtain stable long-term bimetric simulations, hence, as the starting point, we followed closely what people are doing in numerical relativity. Examples of guiding lights are the Einstein Toolkit Lffler2012 and SENR/NRPy+ PhysRevD.97.064036 , which are both using the BSSN formulation. We show that, even though the cBSSN formulation of the EFE, together with the standard gauge and some other technical assumptions, is strongly hyperbolic and therefore its Cauchy problem is well-posed and suitable for the numerical integration PhysRevD.79.104029 , we cannot yet say if the cBSSN formulation of the BFE is also strongly hyperbolic. The reason is that the lapse functions of the two metrics are dependent Alexandrov2013 ; HassanLundkvist2018 ; Kocic_2019 . Their ratio, after a first-order reduction of the equations, involves the dynamical fields algebraically.222The shift vectors are dependent as well, and their relation contains the ratio between the lapses. This ratio appears in the equations in both the metric sectors and, depending on the gauge choice, its first and second spatial derivatives can appear in the equations of one or both sectors. These spatial derivatives affect the hyperbolic structure of the system of PDEs compared to GR, hence the result about the strong hyperbolicity of the cBSSN formulation of the EFE cannot be directly extended to the cBSSN formulation of the BFE. In the paper, we discuss in more detail these issues and stress the similarities and differences compared to GR. Observe that the bimetric features described in the paper, such as the discussion about the gauge fixing, the relations between lapses and shifts, and how these affect the well-posedness of the problem, do not depend on the chosen formulation and need to be accounted for in any formulation of the BFE.

We stress that recasting the BFE as a well-posed Cauchy problem, suitable for numerical integration, is a powerful strategy to obtain interesting physical solutions in BR. For example, since the Birkhoff theorem is not valid in BR Kocic:2017hve , a non-static spherically symmetric system, generically, emits gravitational radiation which is longitudinally polarized, coming from the helicity-0 mode of the massive spin-2 field. This could happen during a spherically symmetric gravitational collapse, or for non-static spherically symmetric black holes. This can provide us with predictions from BR that can be tested against the observational data in the future, which is one of the motivations of our work, inspired by the success of numerical relativity.

Structure of the paper.

In subsection 1.1, we concisely introduce the concepts of well-posedness, and strong hyperbolicity for a system of (first-order) PDEs. In section 2, we very briefly introduce bimetric relativity and review its $N+1$ formulation. In section 3, we introduce the BSSN formulation with its covariant extension and emphasize the differences between BR and GR. In section 4, we discuss gauge fixing in bimetric relativity and qualitatively describe how it affects the hyperbolic structure of the evolution equations. We summarize our results and state our view about what the next challenges in this field are in section 5. The appendix includes explicit equations and technical details.

Notation.

Consider the spatial metrics {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\gamma}, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varphi} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\chi} and their determinants. We shall denote the determinant of a spatial metric with $\Delta$, and define the following notation referring to the metric sectors,

[TABLE]

We denote tensors both with and without their indices, e.g., the metric {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g} or {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}_{\mu\nu}. Greek indices run from 0 to $d-1$, where $d$ is the dimension of spacetime; latin indices run from 1 to $d-1$; boldface indices are spatial Lorentz indices and run from 1 to $d-1$.

@fb@secFB

1.1 Well-posedness and strong hyperbolicity

When dealing with PDEs, it is important to be able to write the mathematical boundary value problem arising from them in a well-posed way. A “mathematical boundary value problem” is a differential problem with some specified boundary or initial data, such as the Dirichlet problem or the Cauchy problem. The definition of a well-posed problem given in (courant1962methods, , p. 226) was introduced for the first time by Hadamard in hadamard1902 and reads: “A mathematical problem which is to correspond to physical reality should satisfy the following basic requirements: (1) The solution must exist. (2) The solution should be uniquely determined. (3) The solution should depend continuously on the data […]”. As stated in (tikhonov1977solutions, , p. 8), the first two requirements characterize the mathematical determinacy of the problem, whereas the third requirement characterizes its physical determinacy and the possibility to apply numerical methods to solve it. We stress that, from the numerical viewpoint, well-posedness is highly desirable, since it entails that small errors in the initial data imply controllable errors in the numerical solution.333Note that the requirement of well-posedness also applies to the case of ordinary differential equations (ODEs). In that case, the well-posedness of the Cauchy problem is established by the Picard–Lindelöf theorem picard ; lindelof ; lindelof2 (also called the Picard–Lipshitz theorem, or the fundamental theorem of ODEs), see also (grant2014theory, , Chapter. 1).

It is important to emphasize that, contrary to Hadamard’s opinion, ill-posed systems can in fact be physically relevant (see e.g. tikhonov1977solutions ). Examples of ill-posed problems can be found in Kabanikhin:2008 . However, in these cases the ill-posedness is acceptable, as it has to do with the physical description of the system. This is not the case in relativity (general or bimetric).

Following (baumgarte2010numerical, , Sec. 11.1), we now give the definition of a strongly hyperbolic system containing first-order time and spatial derivatives. The generalization to systems involving second-order spatial derivatives—as the (c)BSSN system in GR and BR—can be found in Gundlach_2006 . Consider the system of PDEs given by,

[TABLE]

where $u^{i}$ is the $n$-dimensional vector containing the unknown functions to solve for, $S^{i}$ is the source $n$-vector and each of the $A^{k}$ are $n\times n$ matrices of constant coefficients. Consider an arbitrary unit covector $n_{k}$, and define the “principal symbol” or “characteristic matrix” of the system (1) as,

[TABLE]

The system (1) is called “strongly hyperbolic” if, for all unit covectors $n_{k}$, the principal symbol $P^{i}{}_{j}$ has real eigenvalues and a complete set of eigenvectors. Strongly hyperbolic first-order systems of PDEs lead to well-posed Cauchy problems (Kreiss1989ivp, , Theorem 6.2.2).

2 The bimetric field equations and their standard $N+1$ formulation

The BFE can be written Hassan_2013 ,444We choose the same sign convention as in Kocic:2018ddp , for the sign in front of the bimetric interaction potential in the action.

[TABLE]

where {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}G_{g}} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}G_{f}} are the Einstein tensors for the two metrics {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T_{g}} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T_{f}} are two independent stress–energy tensors, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}V_{g}} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}V_{f}} are the bimetric stress–energy tensors that couple the metrics, and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\kappa_{g}} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\kappa_{f}} are the two Einstein gravitational constants. The bimetric stress–energy tensors may or may not satisfy the energy conditions (for the energy conditions, see hawking1973large ). In particular, note that in vacuum, if one of them satisfies the null energy condition, the other does not Baccetti2012 . When one of them does, it can be interpreted as a stress–energy tensor arising from the other spin-2 field. One example of this can be found in Kocic:2017hve .

A fundamental feature about the bimetric stress–energy tensors is that they are non-derivative, i.e., they are not defined in terms of the derivatives of the metrics, but are functions of the square root matrix {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}S}=\left({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}^{-1}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}\right)^{1/2} only, where the principal branch of the matrix square root function is assumed HassanKocic2018 . This allows us to rewrite the BFE (3) as,

[TABLE]

with {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T_{g}^{\mathrm{eff}}}{}^{\mu}{}_{\nu}\coloneqq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}V_{g}}^{\mu}{}_{\nu}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T_{g}}^{\mu}{}_{\nu} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T_{f}^{\mathrm{eff}}}{}^{\mu}{}_{\nu}\coloneqq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}V_{f}}^{\mu}{}_{\nu}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}T_{f}}^{\mu}{}_{\nu}.

The BFE are then formally equivalent to two sets of Einstein equations coupled via the effective stress–energy tensors only. Then, one can recast them in the standard $N+1$ decomposition by following the same steps as in GR. One defines a spacelike hypersurface embedded in the spacetime, where the initial data are specified, and projects the Einstein equations both on the hypersurface and on the direction orthogonal to it (see, e.g., baumgarte2010numerical ; gourgoulhon20123+1 ). The sources are now the sum of the external matter sources and the bimetric stress–energy tensors, whose decomposition has to be computed. This was done independently in Hassan:2011zd ; HassanLundkvist2018 and Kocic:2018ddp following different approaches. The result is a set of evolution equations, a set of constraints similar to those of GR, and the bimetric conservation law (BCL) {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{}_{\mathrm{b}}}=0 (the so-called secondary constraint), which is crucial in eliminating the Boulware–Deser ghost PhysRevD.6.3368 , as explained in HassanLundkvist2018 . Here, we follow the approach in Kocic:2018ddp , to which we refer the reader for the details. The $N+1$ BFE computed in Kocic:2018ddp are written explicitly in Appendix A.2.

The bimetric conservation law {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{}_{\mathrm{b}}}=0 must be preserved in time, therefore \partial_{t}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{}_{\mathrm{b}}}=0. This is called the “preservation of the bimetric conservation law”, and provides a relation between the lapse functions of the two metrics, of the form Alexandrov2013 ; HassanLundkvist2018 ; Kocic_2019

[TABLE]

Unfortunately, the explicit expression for {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W} is very complicated, even in spherical symmetry. See Kocic_2019 and the addendum to this paper in the ancillary files for its expression in the standard $3+1$ and cBSSN formalisms, respectively. Note that the existence of the relation between the lapses is consistent with the fact that bimetric relativity is invariant under the action of a diffeomorphism group acting in the same way on both sectors. In other words, we are free to choose one coordinate system for both metrics. In the $N+1$ decomposition, this translates in the freedom to choose one lapse function and one shift vector only, as explained in more detail in section 4.

3 The covariant BSSN formulation of the bimetric field equations

This section is devoted to the discussion of the BSSN and cBSSN formulation of the standard $3+1$ BFE. We will mainly focus on the differences compared with GR. The explicit equations are presented in Appendix A.1, Appendix A.3 and Appendix A.4.

@fb@secFB

3.1 The BSSN formulation

When rewriting the bimetric $3+1$ equations [eqs. (44, 45, 48)] in the BSSN formulation, the starting point is the definition of the conformal metrics (baumgarte2010numerical, , Sec. 11.5)

[TABLE]

which are assumed to have determinants equal to 1. This renders the conformal metrics tensor densities of weight $-2/3$, and the conformal factors scalar densities of weight $1/6$. The conformal structure implies

[TABLE]

We also decompose the traceless part of the extrinsic curvatures as follows

[TABLE]

Likewise, the conformal traceless extrinsic curvatures are tensor densitites of weight $-2/3$. The indices of the conformal tensors are raised and lowered by the conformal metrics of the corresponding metric sector. In the BSSN formulation, the “conformal connections” are defined as,

[TABLE]

where {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\Gamma}}}^{i}_{jk},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\Gamma}}}^{i}_{jk} are the Christoffel symbols of the conformal metrics {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\gamma}}}_{ij},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\varphi}}}_{ij}. The conformal connections transform as in eq. (11.45) in baumgarte2010numerical . The new dynamical variables in the BSSN formulation are the conformal metrics {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\gamma}}}_{ij},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\varphi}}}_{ij}, the conformal factors {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\phi},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\psi}, the traces of the extrinsic curvatures {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}K},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0muK}}, the conformal traceless parts of the extrinsic curvatures {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{A}}}_{ij},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{A}}}_{\ij} and the conformal connections {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\Gamma}}}^{i},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\Gamma}}}^{i}. With the appropriate transformation rules for these geometrical objects, the BSSN equations are covariant under spatial coordinate transformations not involving the time coordinate PhysRevD.79.104029 .

In order to evolve the system in time, one needs to choose a gauge. The BSSN equations are strongly hyperbolic if one chooses the standard gauge, introduced below, and enforce that {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{A}}}_{ij} is traceless during the evolution PhysRevD.66.064002 ; PhysRevD.70.104004 ; PhysRevD.79.104029 . The standard gauge consists in the 1+$\hskip 1.00006pt\log$ slicing for the lapse {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha} PhysRevLett.75.600 , and the $\Gamma$-driver condition for the shift {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta} PhysRevD.67.084023 ,

[TABLE]

where {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}B}^{i} is an auxiliary variable and $\eta$ is a freely specifiable real constant. As explained in PhysRevD.79.104029 , the $\Gamma$-driver condition (10b, 10c) is not spatially covariant. Suppose we have some initial data on the spacelike hypersurface, written in Cartesian coordinates. We can rewrite them in spherical coordinates by using the transformation rules for the tensors densities and the conformal connections. We then evolve these initial data according to the BSSN equations with the standard gauge, up to some time $t_{f}$. Since the $\Gamma$-driver condition is not covariant, the dynamical variables at $t_{f}$ are not related by the same transformation rules for tensor densities and conformal connections. Therefore, the dynamical variables at $t_{f}$ in Cartesian coordinates and spherical coordinates differ geometrically. This problem can be solved by rewriting the BSSN equations and the standard gauge according to a procedure presented in PhysRevD.79.104029 (see also Alcubierre2011 ) and summarized in subsection 3.3.

In BR, the rewriting of the evolution equations (BEE) and the constraint equations (BCE) in terms of the BSSN dynamical fields mimics the analog computation in GR, whereas the rewriting of the BCL has no analog in GR. Actually, the parts of the equations not involving the bimetric interactions are exactly the same as in GR. The bimetric BSSN constraint and evolution equations are listed in Appendix A.3. The differences compared with GR are:

The presence of the bimetric conservation law {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{}_{\mathrm{b}}}=0. 2. 2.

The fact that the lapse function and the shift vector of one of the metrics {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f} are not freely specifiable. 3. 3.

The effective sources include both the contribution from the external matter sources and the bimetric sources.

The parts involving the bimetric interactions can be rewritten in the BSSN formulation by determining how the spatial vielbein {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}e},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m} of the spatial metrics {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\gamma},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varphi} transform under the conformal change (6). This is discussed in the next subsection.

@fb@secFB

3.2 The BSSN formulation of the bimetric interactions

The $N+1$ decomposition of the BFE as formulated in Kocic:2018ddp relies on the parametrization with respect to the geometric mean metric {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}h}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}\operatorname{\#}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}\coloneqq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}\left({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}^{-1}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}\right)^{1/2} of the metrics {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}.555In index notation we have {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}h}_{ij}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}_{ik}\left[\left({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}^{-1}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}\right)^{1/2}\right]{}^{k}{}_{j}. In this parametrization, the spatial metrics are written in terms of their vielbeins. In matrix notation,

[TABLE]

where the spatial vielbein {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}e} is freely specifiable and the spatial vielbein {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m} is defined as

[TABLE]

where the freely specifiable vielbein {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m_{\mathsf{o}}} of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varphi} is rotated into {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m} by the orthogonal transformation {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{R}}}, and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}} is the spatial part of the Minkowski metric, i.e., the Euclidean metric.666For the sake of clarity, we write {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}} explicitly in every equation where it appears, because it is needed to raise and lower the Lorentz indices. The transformation {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{R}}} is determined by the requirement that the geometric mean {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}h} exists MikicaB.794295 . The operator {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\Lambda}_{\mathrm{s}}} is the spatial part of a Lorentz boost with boost vector {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{v}}}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\lambda}^{-1}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}} and Lorentz factor {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\lambda}=(1+{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}}^{\mathsf{T}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}})^{1/2}. It can be written as {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\Lambda}_{\mathrm{s}}}=({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{I}}}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}}^{\mathsf{T}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}})^{1/2}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{I}}}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}}^{\mathsf{T}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}}/({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\lambda}+1), and the 4-dimensional Lorentz boost itself can be written as

[TABLE]

See Appendix A.1 for more details. In this framework, the real spatial vector {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}, called “separation parameter”, defines the shift vectors {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\beta}}, respectively of the metrics {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}, in terms of the shift vector {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}q} of the geometric mean metric {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}h},

[TABLE]

This is the most general parametrization of the bimetric $N+1$ decomposition Kocic:2018ddp .

Since all the spatial bimetric interactions terms, defined in Kocic:2018ddp and reported in Appendix A.1, depend on the spatial vielbeins, on {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\Lambda}_{\mathrm{s}}} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{R}}}, their rewriting in the BSSN formalism relies on the conformal decomposition of these variables. The conformal decomposition of the vielbeins read,

[TABLE]

which tells us that

[TABLE]

From (16) and the conformal decomposition of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\Lambda}_{\mathrm{s}}} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{R}}}, reported in Appendix A.1, we can derive the BSSN formulation of all the spatial bimetric interactions and sources. The explicit derivations are presented in Appendix A.1.

@fb@secFB

3.3 The covariant extension of the BSSN formulation

As we outlined in subsection 3.1, the BSSN formulation with the standard gauge is not spatially covariant. As described in PhysRevD.79.104029 , this is a problem when comparing the same physical system in different coordinates. Therein, the BSSN formulation was generalized making it spatially covariant, obtaining the cBSSN formulation. Since the computations in PhysRevD.79.104029 do not alter the expressions of the matter sources in the evolution equations, the covariant generalization applies to both metric sectors in bimetric relativity.777The computations in PhysRevD.79.104029 are made in vacuum, but it is straightforward to add an external source and generalize them. As a consequence, the bimetric sources have the same expression as in the BSSN formulation, given by (38, 40, 43).888As a general rule in this context, every computation that does not concern the matter sources in the evolution and constraint equations in GR can be directly translated in BR without modification.

Having a covariant version of the BSSN formulation is important in BR, since it allows us to safely use spherical coordinates. Since the Birkhoff theorem is not valid in BR, see, e.g, Kocic:2017hve , a spherically symmetric solution of the BFE does not need to be static. From one side, this may be interpreted as an undesired feature of the theory; on the other hand, it makes the study of spherically symmetric systems much more interesting in BR than in GR, as discussed in the Introduction. Specifically, spherically symmetric systems in BR emit longitudinally polarized gravitational radiation, which can be tested against observational data.

In addition, using spherical coordinates made it possible to compute both the ratio between the lapses {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W} appearing in (5) and the spatial Lorentz vector {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}} appearing in (14) in the difference between the shifts (see Kocic:2018ddp ; Kocic_2019 and Appendix A.5 for more details). Note that {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}} is also known in the most general $\beta_{(1)}$-model Hassan:2011tf ; HassanLundkvist2018 , where the explicit expression for {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\mathsfit{n}}{\mkern 1.0mu}}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m}^{-1}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\lambda}^{-1} in (14) is computed.999Knowing {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\mathsfit{n}}{\mkern 1.0mu}}, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}=\big{[}1-({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\mathsfit{n}}{\mkern 1.0mu}})^{\mathsf{{\scriptscriptstyle T}}}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\mathsfit{n}}{\mkern 1.0mu}})\big{]}^{-1/2}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\mathsfit{n}}{\mkern 1.0mu}}.

In the BSSN formulation, the determinant of the conformal metric {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\gamma}}} is taken to be $1$, making it a scalar rather than a scalar density. This also alters the transformation properties of the metric and the extrinsic curvature, making them tensor densities. In the cBSSN formulation, nothing is assumed on the transformation properties of the determinant of the conformal metric. The new conformal decomposition of the metrics and the extrinsic curvatures becomes,

[TABLE]

which is the same as before, without the restriction on the determinants {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\Delta}}} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\Delta}}}, and

[TABLE]

The tensors {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{A}}} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{A}}} are not traceless, as can be seen by comparing (18) with (8). The conformal connections in (9) are made covariant by introducing two background connections {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\Gamma}}_{\mathsf{B}}}^{i}_{jk},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\Gamma}}_{\mathsf{B}}}^{i}_{jk} Garfinkle_2008 . It is possible, but not necessary, to introduce two background metrics whose Levi–Civita connections serve as the background connections. We define the new dynamical variables,

[TABLE]

Since the difference between two Christoffel symbols is a tensor, our set of dynamical variables includes tensors only, i.e., {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\gamma}_{ij},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varphi}_{ij},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\phi},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\psi},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{A}}}_{ij},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{A}}}_{ij},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{K}}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{K}}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\Lambda}}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\Lambda}}}. The standard gauge (10) in the cBSSN formalism reads,

[TABLE]

which is manifestly covariant. The explicit bimetric cBSSN constraint and evolution equations are written in Appendix A.4.

We emphasize that the background connections are completely arbitrary, and in GR, there is no preferred connection to be chosen. In bimetric relativity, a third metric {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}h} is defined.101010Actually, in the space of pseudo-Riemannian metrics built on our manifold, we have a path of metrics connecting {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}, corresponding to a geodesic of the trace metric and parameterized by h_{\alpha}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}\left({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}^{-1}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}\right)^{\alpha},\;\alpha\in\mathbb{R} (see PhysRevD.97.084022 ; Kocic:2018ddp ). Hence, we can choose the Levi–Civita connection of the conformal spatial metric {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\accentset{\circ}{{\chi}}}} to define the covariant conformal connection in (19). The consequences of this choice are described in more detail in Torsello_2019 .

In the cBSSN formulation, the determinants {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\Delta}}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\Delta}}} of the conformal metrics and the traces {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{A}}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{A}}} are left undetermined, making the cBSSN evolution equations in (57, 59) incomplete PhysRevD.79.104029 . We must choose how these quantities evolve, in order to be able to evolve the full system. Following PhysRevD.79.104029 , we choose \partial_{t}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{A}}}=\partial_{t}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{A}}}=0. Regarding the determinants, there are two natural choices, referred to as “Lagrangian” and “Eulerian” PhysRevD.71.104011 ; Brown_2008 ; PhysRevD.79.104029 , given by

[TABLE]

with analog expressions for the {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}-sector. In BR, we need to specify the evolution of the determinants and the traces in both sectors. More details on this can be found in Torsello_2019 . The expression for \partial_{t}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\Lambda}}}^{i} in (20c) is explicitly substituted with the Lagrangian formulation of (57e).

4 Well-posedness and gauge choices

As established in Kocic:2018ddp , the three lapses and shifts of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}h} are related by,

[TABLE]

where {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\lambda}=\sqrt{1+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}^{\mathsf{{\scriptscriptstyle T}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}H} is the lapse of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}h}. We are free to choose one lapse function and one shift vector, exactly as in GR; the other four quantities are determined by (22).

When imposing a gauge in BR, it can be written in terms of any of the three lapse functions or the shift vectors. Suppose, for example, that we choose the standard gauge with respect to {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g} in (20). This gauge can be rewritten in terms of the lapse functions {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\alpha}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}H} and the shift vectors {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\beta}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}q} by using (22). Hence, we can impose the standard gauge with respect to {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g} by gauge fixing, say, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}H} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\beta}}. We say that we “choose a gauge condition with respect to a metric”, to emphasize that the geometry of the slicing is determined by that metric. In addition, we say that we “gauge fix” one of the lapses and shifts. The same gauge choice (or gauge condition) can be expressed via different, but equivalent, gauge fixings. It follows that, in BR one can have “mixed” gauges, i.e., one can choose the 1+$\hskip 1.00006pt\log$ slicing with respect to {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\chi}, and the $\Gamma$-driver condition with respect to {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varphi}. In this case, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}h} would determine the time slicing, whereas {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f} would determine the spatial gauge. If these gauges are singularity avoiding or horizon penetrating for any of the metrics remains an open question. See Torsello_2019 for a study of the “mean gauges”, i.e., the gauge choices with respect to the mean metric {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}h}.

In GR, the cBSSN formulation is strongly hyperbolic if one chooses the standard gauge and fulfills some other technical conditions PhysRevD.79.104029 . In BR, the well-posedness of the evolution equations involves both of the metric sectors. Suppose that we fix the lapse and shift of one metric to be determined by the standard gauge. Now, the bimetric source {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}J^{\mathrm{b}}}^{i}{}_{j} appearing in the evolution equation for the conformal extrinsic curvature (28) contains the ratio of the lapses {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W}. The general explicit expression of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W} is not known, but it can be computed explicitly in spherical symmetry (Kocic_2019 and the addendum to this paper). In that case, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W} is a lengthy expression—roughly, it fills two pages—which depends on the radial derivatives of the dynamical fields. If the radial derivatives can not be replaced by algebraic expressions, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W} will affect the characteristic matrix in (2) and alter the hyperbolic structure of the equations in the {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}-sector compared with GR. Following the procedure described in ALCUBIERRE200576 , which promotes the logarithmic radial derivatives of the dynamical fields to be new dynamical variables, thus achieving a first-order reduction of the system, one ends up with an expression for {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W} which only depends on three radial derivatives, namely,

[TABLE]

Here, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}^{\textbf{1}}} is the only nonzero component of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{A}}_{2}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{A}}_{2}} are the {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{A}}}^{\theta}{}_{\theta},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{A}}}^{\theta}{}_{\theta} components of the conformal extrinsic curvatures (see Appendix A.5 for more details). By using the two momentum constraints (66b, 66e) and the bimetric conservation law (A.5) rewritten according to the procedure in ALCUBIERRE200576 , these three derivatives can be substituted with expressions involving the dynamical fields only algebraically. In more detail, the BCL can be solved for \partial_{r}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}^{\textbf{1}}}, and the momentum constraints can be solved for \partial_{r}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{A}}_{2}}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{K}}}/3) and \partial_{r}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{A}}_{2}}+{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{K}}}/3).111111Note that, in this way, we can freely specify {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}} on the initial hypersurface, and the value of these three derivatives will depend on this choice. This means that the ratio of the lapses {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W} is a purely algebraic expression in the dynamical fields, and, as such, does not enter the characteristic matrix in (2) and does not alter the hyperbolicity of the equations in the {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}-sector compared with GR. Therefore, if we choose the standard gauge (20) for {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\alpha} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta}, the equations in the {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}-sector are strongly hyperbolic.

Note that the algebraicity of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W} in the dynamical fields is compatible with the fact that, in the case of static spherically symmetric black hole solutions in BR (see PhysRevD.85.124043 ; Brito:2013xaa ; PhysRevD.96.064003 and reference therein) the function describing the ratio of the lapses, $\tau$ in PhysRevD.96.064003 , is determined by an algebraic relation. This relation corresponds precisely to the PBCL \partial_{t}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{}_{\mathrm{b}}}=0. Also, this confirms that the bimetric stress–energy tensors (of which the spatial bimetric interactions are the projections Kocic:2018ddp ) are non-derivative and cannot spoil the hyperbolic structure of the equations compared with GR (in spherical symmetry).

On the other hand, we need to consider the {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}-sector as well. The equations in the {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}-sector formally appear the same as in the {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}-sector, but now the lapse and shift of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f} are determined by (22). As a consequence, the ratio of the lapses appears wherever the lapse and shift of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f} appear. Since there are terms involving first and second spatial derivatives of the lapse and shift, they contain first and second spatial derivatives of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W}. Hence, they contain the spatial derivatives of all the dynamical fields, and we cannot eliminate all of them by using the constraints. Therefore, the hyperbolic structure of the PDEs in the {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}-sector is drastically different compared with GR. This means that one cannot carry over the GR results to BR, and additional steps are needed towards a definite answer regarding the strong hyperbolicity of the equations.

Instead of using (22) to replace the lapse of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}, we could equivalently use it to determine its value on the initial hypersurface, and impose the preservation in time of the PBCL (5) by setting its time derivative to 0. This gives an evolution equation for the lapse of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}, which becomes a dynamical variable.121212More in general, this is an evolution equation for the lapse that we do not gauge fix. Hence, with this choice, the principal symbol of the system of PDEs is largely different than in GR. For example, the evolution equation for the lapse of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f} involves the time derivative of the ratio of the lapses {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W} which nontrivially alters the hyperbolic structure compared with GR. See (Torsello1390334, , Sec. 6.2, p. 103) for more details.

The analysis above holds in the case of spherical symmetry, which is presented in Appendix A.5. This study suggests that the ratio of the lapses, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W}, always contains some spatial derivatives of the dynamical fields which can be eliminated by using the constraints, but the hyperbolicity is altered compared to GR by the spatial derivatives of the lapses and shifts involving {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W}. Therefore, the computation of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W} in the general case is a prerequisite for the study of the well-posedness of any formulation of the BFE.

A possible gauge choice, which preserves the symmetry of the equations between {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f} and modifies the hyperbolic structures of both sectors in a more symmetric—and hopefully better-behaved—way, is to fix the lapse and the shift of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}h}. In this case, we see from (22) that {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W} appears in the spatial derivatives of the lapses and shifts of both {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}, thus modifying the hyperbolic structure in both sectors. This is investigated in more detail in Torsello_2019 , where both the standard gauge and the maximal slicing for {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}H} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}q} are computed.

5 Summary and outlook

In this paper, we presented the covariant BSSN formalism of the bimetric field equations. We emphasized why this formulation is important in bimetric relativity and we stressed the differences with the analogous formulation of the Einstein equations, summarized below.

In addition to the Hamiltonian and momentum constraints for both metrics, there is an additional bimetric conservation law (the so-called secondary constraint) {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{}_{\mathrm{b}}}=0 that, in the free evolution scheme, has to be solved for the initial data on the spacelike hypersurface. 2. 2.

The sources in the equations include both the external matter sources and the bimetric sources in (28). After a first-order reduction of the PDEs, the bimetric sources do not contain the derivatives of the dynamical fields. Hence, they do not alter the hyperbolic structure of the equations compared to GR. 3. 3.

Bimetric relativity is diffeomorphism invariant. This provides us with the possibility to choose one lapse function and one shift vector of any of the metrics, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\gamma}, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varphi} or their geometric mean {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\chi}. The remaining lapses and shifts are determined by (22). 4. 4.

The relation between the lapses is established in Alexandrov2013 ; HassanLundkvist2018 by imposing the preservation of the bimetric conservation law in time, \partial_{t}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{}_{\mathrm{b}}}=0, and it is computed explicitly in Kocic_2019 for spherically symmetric spacetimes in the standard $3+1$ formulation. The expression in the covariant BSSN formalism is given in the addendum to this paper. The ratio between the lapses, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W}, is a lengthy algebraic expression in terms of the dynamical fields and their spatial derivatives. Since the evolution equations involve the spatial derivatives of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W}, the hyperbolic structure of the system of PDE is different compared with the corresponding equations in GR. The hyperbolic structure is changed in either one of the two metric sectors, or in both, depending on the gauge choice.

Other than these four differences, the system is analogous to the covariant BSSN formulation of the Einstein field equations presented in PhysRevD.79.104029 . In particular, from the viewpoint of numerical relativity, the bimetric fields equations can be tackled numerically in the same way as the Einstein field equations. However, since we showed that the results in GR cannot be carried over to BR in a straightforward way, and the well-posedness of the problem is not proved yet, we do not know how successful this can be. Nonetheless, the equations do offer some stimulating challenges:

The computation of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}} is necessary to be able to solve the bimetric equations in any formulation. At present, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W} is only computed under the assumption of spherical symmetry, whereas {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}} is computed in spherical symmetry and in the most general $\beta_{(1)}$-model; we refer the reader to Kocic_2019 ; Hassan:2011tf ; HassanLundkvist2018 for more details. 2. 2.

Investigating if the bimetric covariant BSSN evolution equations, together with a suitable gauge, are strongly hyperbolic is of great importance and depends on the computation of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W}. Since the latter is known in spherical symmetry, one can study the hyperbolicity of the evolution equations in (70). 3. 3.

The choice of the gauge is essential in bimetric relativity (as it is in GR as well). In Torsello_2019 , we study some possible gauge choices which alter the hyperbolic structure of the evolution equations in both sectors. In particular, we investigate the gauge fixing on the geometric mean metric {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}h}. 4. 4.

The challenge is to integrate the bimetric BSSN equations numerically in spherical symmetry, e.g., for a gravitational collapse of matter or a non-static black hole solution. The numerical computation of both {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}W} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}} significantly reduces the accuracy of the integration. We have written a Mathematica/C++ code to perform the simulations, see Kocic:2019gxl for results obtained using the standard $3+1$ equations. The results concerning the cBSSN formulation will be the subject of another work. We remind the reader that, since the Birkhoff theorem is not valid in bimetric relativity (see, e.g., Kocic:2017hve ), the spherically symmetric case is very interesting to study. One can look for both vacuum and non-vacuum spherical solutions with nontrivial dynamics, and perhaps gravitational wave emission. This can potentially lead to results that could directly be compared against observational data.

@fb@secFB

Acknowledgments

It is a pleasure to thank Anders Lundkvist and Fawad Hassan for many fruitful discussions. We thank Anders Lundkvist and Giovanni Camelio for providing valuable remarks after a careful reading of the paper.

Appendix A Explicit equations and computations

@fb@secFB

A.1 The BSSN decomposition of the bimetric interactions and sources

Consider the spatial parts {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\gamma}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}e}^{\mathsf{{\scriptscriptstyle T}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}e} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varphi}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m_{\mathsf{o}}}^{\mathsf{{\scriptscriptstyle T}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m_{\mathsf{o}}} of two Lorentzian metrics {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}. In MikicaB.794295 , it is established that the existence of the real square root \left({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}^{-1}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}\right)^{1/2} implies,

[TABLE]

To be more precise, the freely specifiable spatial vielbein {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m_{\mathsf{o}}} is used to compute the vielbein {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{R}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m_{\mathsf{o}}} such that the spatial part {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\chi} of the geometric mean metric {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}h}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}\operatorname{\#}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f} is given by {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\chi}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}e}^{\mathsf{{\scriptscriptstyle T}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\Lambda}_{\mathrm{s}}}({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{R}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m_{\mathsf{o}}}). This is obtained by imposing (24c) and solving it for the Euclidean orthogonal transformation {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{R}}} in terms of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\Lambda}_{\mathrm{s}}} and the vielbeins {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}e},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m_{\mathsf{o}}}. Such a solution always exists, as it is part of the polar decomposition of the invertible matrix {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}} MikicaB.794295 ; Hassan:2014gta [see (3.7) ]. For the sake of simplicity, we define the new vielbein of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varphi} to be {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m}\coloneqq{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{R}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m_{\mathsf{o}}}; we have the freedom to do that since {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m_{\mathsf{o}}}^{\mathsf{{\scriptscriptstyle T}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{R}}}^{\mathsf{{\scriptscriptstyle T}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{R}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m_{\mathsf{o}}}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m_{\mathsf{o}}}^{\mathsf{{\scriptscriptstyle T}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m_{\mathsf{o}}}, implying that {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varphi} is blind to this choice. The interaction terms are not affected as well, since they always contain both {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\Lambda}_{\mathrm{s}}} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{R}}}, irrespective of this choice. The matrix {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\Lambda}_{\mathrm{s}}} explicitly appears in them. On the contrary, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{R}}} does not appear explicitly, but it is taken into account implicitly inside {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m}.

We define the bimetric interactions as Kocic:2018ddp ,

[TABLE]

where $e_{n}(X)$ are the elementary symmetric polynomials of the linear operator $X$,

[TABLE]

and $Y_{n}(X)$ is defined as,

[TABLE]

See Kocic:2018ddp for more details about the properties of $e_{n}(X)$ and $Y_{n}(X)$. Note that $d$ is the dimension of the spacetime, i.e., $d=$ $N+1$. Hence, some terms in the summations will be zero. The $\beta_{(n)}$ parameters are $d+1$ real dimensionless constants appearing in the bimetric interaction potential, together with the energy scale $m$ Hassan:2011zd . We define the bimetric sources (respectively, the bimetric energy densitites, the bimetric currents and the bimetric spatial stress–energy tensors) as Kocic:2018ddp ,

[TABLE]

where the summation $-\,{m^{d}{\sum_{n=0}^{d}}\beta_{n}}$ is understood in front of all the bimetric sources. Note the relation between the two bimetric currents {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}j^{\mathrm{b}}}_{i},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0muj}{}^{\mathrm{b}}}_{i}, which implies the relation (A.35) between the momentum constraints .

Here we compute the expressions for the bimetric interaction and sources in the (c)BSSN formalism. We require that the symmetrization condition (24) holds for the BSSN variables as well. Since the shifts are the same in the BSSN formalism, we require conditions (24a, 24b) to stay the same. The condition (24c) should instead lead to its analog in the BSSN formalism,

[TABLE]

where {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{{\accentset{\ast}{{\Lambda}}}}_{\mathrm{s}}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{{\accentset{\ast}{{\mathrm{R}}}}}} are the BSSN counterparts of the spatial part of the Lorentz boost (3.8) and the orthogonal transformation in (3.7) , whose expression is unknown yet.

We start by computing the conformal decomposition of the objects in the Lorentz frame. The requirement that (24a) stays the same implies,

[TABLE]

where the scalar $\xi$ accounts for the conformal decomposition of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\lambda}^{-1}. It follows that,

[TABLE]

We apply {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}^{\mathsf{{\scriptscriptstyle T}}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}} to (31) and obtain,

[TABLE]

which is equivalent to

[TABLE]

Hence, in general, we can rescale {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}} as in (31) with a generic $\xi$ when we recast the equations into the (c)BSSN formulation, as long as we satisfy (33). However, there is no need to rescale it since this is an unnecessary complication. Indeed, we can always satisfy (31) and (33) by choosing $\xi=1$, which implies {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\accentset{\ast}{{\boldsymbol{\mathrm{p}}}}}}. It immediately follows,

[TABLE]

Let us now compute the BSSN version of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{R}}_{\mathrm{o}}}, introduced in (3.2). From (16, 34c) and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{m}}_{\mathrm{o}}}=\mathrm{e}^{-2{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\psi}}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m_{\mathsf{o}}} [which follows from (16)] it follows that,

[TABLE]

which implies,

[TABLE]

Using (16, 34c, 36), the spatial part of the symmetrization condition (24c) can be written as,

[TABLE]

i.e., if {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\chi} is symmetric, its BSSN counterpart {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\accentset{\circ}{{\chi}}}} is also symmetric, as desired.

In light of (16, 34, 36), we compute the BSSN decomposition of bimetric interactions,

[TABLE]

We note the following property of the elementary symmetric polynomials,

[TABLE]

where $f$ is a scalar field. We use this property to compute the bimetric interactions in terms of the BSSN variables,

[TABLE]

From (25g, 38b, 40c, 40d) it follows,

[TABLE]

We compute the bimetric sources in (28) in terms of the BSSN variables,

[TABLE]

We remind the reader that $N=d-1$ is the dimension of the spacelike hypersurface, in our case $N=3$, and $I$ is the $N\times N$ identity. Note that both the bimetric interactions and sources can be computed starting with the physical vielbeins and metrics {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}e},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}m},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\gamma},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varphi}, or with the conformal vielbeins and metrics {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{e}}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{m}}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\gamma}}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{\varphi}}} and the conformal factors {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\phi},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\psi}.

@fb@secFB

A.2 The bimetric standard $N+1$ equations

The bimetric standard $N+1$ evolution equations read Kocic:2018ddp ,

[TABLE]

with {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\partial_{\perp}}\coloneqq\partial_{t}-\mathscr{L}_{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta}, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\widetilde{\partial}_{\perp}}\coloneqq\partial_{t}-\mathscr{L}_{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\beta}}, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\beta} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\mkern 1.0mu\widetilde{\mkern-1.0mu\beta}} being the shift vectors of the two metrics. The bimetric standard $N+1$ constraint equations are,

[TABLE]

The effective sources are the sum of the bimetric sources given by (28), and the external matter sources,

[TABLE]

Note that the relation between the two bimetric currents in (28), implies that the two momentum constraints are related to each other,131313Note, however, that solving one momentum constraint does not imply that the other is automatically satisfied. Both of them need to be solved independently.

[TABLE]

The bimetric conservation law (BCL), in its asymmetric and symmetric form, reads

[TABLE]

@fb@secFB

A.3 The bimetric BSSN equations

We now write down the bimetric BSSN equations, applying the procedure in (baumgarte2010numerical, , Sec. 11.5) to (44, 45, 48), using (43) for the bimetric sources.

The bimetric BSSN $3+1$ evolution equations for the $g$-sector are,

[TABLE]

where

[TABLE]

and for the $f$-sector

[TABLE]

where

[TABLE]

The bimetric BSSN constraint equations are

[TABLE]

Now we compute the BSSN decomposition of the bimetric conservation law (48). We note that, due to (6), the relation between the two conformally related covariant derivatives of a vector field $X^{i}$ is,

[TABLE]

which implies

[TABLE]

with analogous formulas for the {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}-sector. Using the definition of the conformal extrinsic curvatures (18) in the cBSSN formalism, the bimetric conservation law becomes

[TABLE]

Finally, we insert the expressions computed in Appendix A.1 for the quantities appearing in (56), which are written in terms of the BSSN variables.

@fb@secFB

A.4 The bimetric covariant BSSN equations

At this point, we apply the method described in PhysRevD.79.104029 to the bimetric BSSN equations in Appendix A.3. The bimetric interactions and sources are not changed.

These are the cBSSN $3+1$ bimetric evolution equations for the $g$-sector, assuming that the background connections do not depend on time,

[TABLE]

where {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{D}}_{\mathsf{B}}} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{R}}_{\mathsf{B}}} are the covariant derivative and the Riemann tensor of the background geometry, and

[TABLE]

Note that all Lie derivatives in (58) are Lie derivatives of tensors, not tensor densities as in (50). Equation (58e) is an identity for the Ricci tensor of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\gamma}}}_{ij} in terms of the background geometry PhysRevD.79.104029 ; PhysRevD.85.084004 (see Appendix A of PhysRevD.85.084004 for the proof). The superscript $\mathrm{TF}$ means “trace-free”. Note that all the traces are with respect to {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{\gamma}}} except the trace of {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}J^{\mathrm{eff}}}, which is with respect to {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\gamma}. For the {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}f}-sector we have,

[TABLE]

where the same conventions and notations as for the {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}g}-sector are implied. It holds,

[TABLE]

The bimetric cBSSN constraint equations are

[TABLE]

The cBSSN BCL is the same as in one in (56), since the bimetric interactions (38, 40) are not changed, and the bimetric sources are the same as in (43).

@fb@secFB

A.5 The bimetric covariant BSSN equations in spherical symmetry

The equations in this Appendix were computed with the Mathematica package $\mathtt{bimEX}$ TORSELLO2020106948 , which can compute the bimetric cBSSN equations for any desired ansatz.

We write down the cBSSN equations in spherical symmetry, assuming the following ansatz,

[TABLE]

The background geometries for both {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\gamma} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\varphi} are chosen to be the spatial part of the flat metric in spherical coordinates,141414Note that this is not the Minkowski metric in the Lorentz frame, whose spatial part remains {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{\delta}}}=\operatorname{diag}\left[1,1,1\right].

[TABLE]

From now on, we will assume the time and radial dependence of all the fields. We define the function,

[TABLE]

and the “shifted elementary symmetric polynomials” doi:10.1063/1.5100027 ,

[TABLE]

where the $\beta_{(n)}$ are five real dimensionless constants appearing in the bimetric interaction potential Hassan:2011zd .

Consistency conditions between the various equations imply that {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}q}^{\theta}}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}^{\textbf{2}}}=0, and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}q}^{\phi}},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}^{\textbf{3}}} do not appear explicitly into the equations, so we can set them to zero without losing generality. Also, {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{R}}} is the identity and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{m}}}={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{m}}_{\mathrm{o}}}, which simplifies the computations considerably.

The constraint equations (61) read,

[TABLE]

where the bimetric energy density and currents are given by

[TABLE]

We can solve {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}^{\textbf{1}}} from both (66b) and (66e), obtaining respectively

[TABLE]

The asymmetric—and simpler—version of the BCL (56a) is,

[TABLE]

The evolution equations (57, 59), modified to get the evolution of the components {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{A}}}^{i}{}_{j},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{A}}}^{i}{}_{j} rather than {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{A}}}_{ij},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{A}}}_{ij}, reduce to

[TABLE]

where

[TABLE]

and

[TABLE]

The evolution equation for {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\boldsymbol{\mathrm{p}}}^{\textbf{1}}} is obtained from the cBSSN version of (A.8) in Kocic:2018ddp ,

[TABLE]

A direct comparison between (A.5) and (2.10) in Kocic_2019 , reveals their equivalence.

Equations (66, 67, 68, A.5, 70, 71, 72, A.5) reduce to the standard $3+1$ equations in spherical symmetry presented in Kocic:2018ddp after imposing,

[TABLE]

and appropriately using the Hamiltonian constraints in the evolution equations for the traces {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widebar{K}}} and {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\widehat{K}}}, and the momentum constraints in the evolution equations for the conformal factors {\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\phi},{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\psi}.