Covariant BSSN formulation in bimetric relativity

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 roadmap towards numerical bimetric relativity.

We provide corrections to a typo and an error in the original paper.

Correction to the error in appendix A.1
This section replaces the first part of appendix A. 1, up to (A.25). Consider the spatial parts γ = e T δe and ϕ = m o T δm o of two Lorentzian metrics g, f . In [1], it is established that the existence of the real square root g −1 f 1/2 implies, orthogonal transformation R in terms of Λ s and the vielbeins e, m o . Such a solution always exists, as it is part of the polar decomposition of the invertible matrix R o [1,2] (see (3.7) in the main text). For the sake of simplicity, we define the new vielbein of ϕ to be m := Rm o ; we have the freedom to do that since m o T R T δRm o = m o T δm o , implying that ϕ is blind to this choice. The interaction terms are not affected as well, since they always contain both Λ s and R, irrespective of this choice. The matrix Λ s explicitly appears in them. On the contrary, R does not appear explicitly, but it is taken into account implicitly inside m.
We define the bimetric interactions as [3], See [3] for more details about the properties of e n (X) and Y n (X). Note that d is the dimension of the spacetime, that is, d = N + 1. Hence, some terms in the summations will be zero. The β (n) parameters are d + 1 real dimensionless constants appearing in the bimetric interaction potential, together with the energy scale m [4]. We define the bimetric sources (respectively, the bimetric energy densitites, the bimetric currents and the bimetric spatial stress-energy tensors) as [3], where the summation − m d d n=0 β (n) is understood in front of all the bimetric sources. Note the relation between the two bimetric currents j b i , j b i , which implies the relation (A. 35) between the momentum constraints in the main text.
Here we compute the expressions for the bimetric interaction and sources in the (c)BSSN formalism. We require that the symmetrization condition (2.1) holds for the BSSN variables as well. Since the shifts are the same in the BSSN formalism, we require conditions (2.1a) and (2.1b) to stay the same. The condition (2.1c) should instead lead to its analog in the BSSN formalism, where * Λ s , * R are the BSSN counterparts of the spatial part of the Lorentz boost (3.8) and the orthogonal transformation in (3.7) in the main text, whose expression is unknown yet.
We start by computing the conformal decomposition of the objects in the Lorentz frame. The requirement that (2.1a) stays the same implies, where the scalar ξ accounts for the conformal decomposition of pλ −1 . It follows that, We apply p T δ to (2.8) and obtain, 9) which is equivalent to (2.10) Hence, in general, we can rescale p as in (2.8) with a generic ξ when we recast the equations into the (c) BSSN formulation, as long as we satisfy (2.10). However, there is no need to rescale it since this is an unnecessary complication. Indeed, we can always satisfy (2.8) and (2.10) by choosing ξ = 1, which implies p = * p. It immediately follows, which implies (see (A.23) in the main text), (2.12) Using (3.11) in the main text and (2.11c), (2.12) and m o = e −2ψ m o (which follows from (3.11) in the main text), the spatial part of the symmetrization condition (2.1c) can be written as, χ = e T δΛ s Rm o = e 2(φ+ψ) e T δΛ s R m o =: e 2(φ+ψ) • χ = χ T = (e T δΛ s Rm o ) T = e 2(φ+ψ) (e T δΛ s R m o ) T =: e 2(φ+ψ) • χ T , (2.13) that is, if χ is symmetric, its BSSN counterpart • χ is also symmetric, as desired.

Introduction
The Hassan-Rosen bimetric theory, or bimetric relativity (BR), is an extension of general relativity (GR) describing the nonlinear interaction of two metrics g µν and f µν [1,2]. Its unambiguous definition and spacetime interpretation are established in [3], and the complete Hamiltonian analysis of it is performed in [4]. 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 [5], 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 [6][7][8][9] and a framework where the gravitational origin of dark matter can be studied [10][11][12][13]. 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 2 . 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 [14,15]. This implies that a test particle coupled to the metric g µν follows the geodesics determined by this metric, exactly as in GR. The interaction with f µν is only experienced indirectly. The key difference is that the metric g µν 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 [16].
The field equations governing the dynamics of the two metrics have been solved analytically in some important cases (see [16] 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 [17,18], and in the presence of external matter sources under other conditions [17].
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. [19]) 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 [20]. The concept of wellposedness for first-order systems of PDEs is introduced in section 1.1.
In GR, the EFE can be recasted in a variety of well-posed forms (see, e.g. [21, chapter 5]), one of them being the Baumgarte, Shapiro, Shibata, Nakamura, Oohara and Kojima (BSSNOK or, more commonly, BSSN) formulation [22][23][24][25][26] (see also [27] for a pedagogical introduction), with its covariant generalization [28]. One begins by rewriting the EFE as a Cauchy problem. This can be accomplished doing a N + 1 decomposition, where the EFE 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 [29] it is proven that the constraints stay satisfied in the free evolution scheme, and the same holds in BR [30]. 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 [31]. 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. section 11.4.5 in [32]. 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 [33] and SENR/NRPy + [34], 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 [28], 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 [4,35,36]. Their ratio, after a first-order reduction of the equations, involves the dynamical fields algebraically 3 . This ratio appears in the evolution 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 [37], 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 section 1.3, 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 γ, ϕ and χ and their determinants. We shall denote the determinant of a spatial metric with ∆, and define the following notation referring to the metric sectors, We denote tensors both with and without their indices, e.g. the metric g or g µν . 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.

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 [20, p 226] was introduced for the first time by Hadamard in [38] 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 [39, 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 implies controllable errors in the numerical solution 4 . It is important to emphasize that, contrary to Hadamard's opinion, ill-posed systems can in fact be physically relevant (see e.g. [39]). Examples of ill-posed problems can be found in [44]. 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 [27, section 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 [45]. Consider the system of PDEs given by, 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 × n matrices of constant coefficients. Consider an arbitrary unit covector n k , and define the 'principal symbol' or 'characteristic matrix' of the system (1.1) as, (1.2) The system (1.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 [46, theorem 6.2.2].

The bimetric field equations and their standard N + 1 formulation
The BFE can be written [47] where G g and G f are the Einstein tensors for the two metrics g and f , T g and T f are two independent stress-energy tensors, V g and V f are the bimetric stress-energy tensors that couple the metrics, and κ g and κ 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 [48]). In particular, note that in vacuum, if one of them satisfies the null energy condition, the other does not [49]. 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 [37].
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 S = g −1 f 1/2 only, where the principal branch of the matrix square root function is assumed [3]. This allows us to rewrite the BFE (2.1) as, 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. [27,32]). 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 [1,4] and [31] 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) C b = 0 (the so-called secondary constraint), which is crucial in eliminating the Boulware-Deser ghost [50], as explained in [4]. Here, we follow the approach in [31], to which we refer the reader for the details. The N + 1 BFE computed in [31] are written explicitly in appendix A.2.
The bimetric conservation law C b = 0 must be preserved in time, therefore ∂ t C b = 0. This is called the 'preservation of the bimetric conservation law' (PBCL), and provides a relation between the lapse functions of the two metrics, of the form [4,35,36]  Unfortunately, the explicit expression for W is very complicated, even in spherical symmetry. See [36] and the addendum to this paper in the supplementary material (available at stacks.iop. org/CQG/37/025013/mmedia) 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.

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 appendices A.1, A.3 and A.4.

The BSSN formulation
When rewriting the bimetric standard 3 + 1 equations (equations (A.32a), (A.33a) and (A.36a)) in the BSSN formulation, the starting point is the definition of the conformal metrics [27, section 11.5] 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 We also decompose the traceless part of the extrinsic curvatures as follows 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, where Γ i jk ,Γ i jk are the Christoffel symbols of the conformal metrics γ ij ,φ ij . The conformal connections transform as in equation (11.45) in [27]. The new dynamical variables in the BSSN formulation are the conformal metrics γ ij ,φ ij , the conformal factors φ, ψ, the traces of the extrinsic curvatures K,K, the conformal traceless parts of the extrinsic curvatures A ij ,Â ij and the conformal connections Γ i ,Γ i . With the appropriate transformation rules for these geometrical objects, the BSSN equations are covariant under spatial coordinate transformations not involving the time coordinate [28]. 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 A ij is traceless during the evolution [25,26,28]. The standard gauge consists in the 1 + log slicing for the lapse α [51], and the Γ-driver condition for the shift β [52], where B i is an auxiliary variable and η is a freely specifiable real constant. As explained in [28], the Γ-driver condition (3.5b) and (3.5c) 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 Γ-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 [28] (see also [53]) and summarized in section 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: 1. The presence of the bimetric conservation law C b = 0. 2. The fact that the lapse function and the shift vector of one of the metrics g and f are not freely specifiable. 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 e, m of the spatial metrics γ, ϕ transform under the conformal change (3.1). This is discussed in the next subsection.

The BSSN formulation of the bimetric interactions
The N + 1 decomposition of the BFE as formulated in [31] relies on the parametrization with respect to the geometric mean metric h = g#f := g g −1 f 1/2 of the metrics g and f 6 . In this parametrization, the spatial metrics are written in terms of their vielbeins. In matrix notation, where the spatial vielbein e is freely specifiable and the spatial vielbein m is defined as where the freely specifiable vielbein m o of ϕ is rotated into m by the orthogonal transformation R, and δ is the spatial part of the Minkowski metric, i.e. the Euclidean metric 7 . The transformation R is determined by the requirement that the geometric mean h exists [54]. The operator Λ s is the spatial part of a Lorentz boost with boost vector v = λ −1 p and Lorentz factor λ = (1 + p T δp) 1/2 . It can be written as Λ s = (I + pp T δ) 1/2 = I + pp T δ/(λ + 1), and the 4-dimensional Lorentz boost itself can be written as See appendix A.1 for more details. In this framework, the real spatial vector p, called 'separation param eter', defines the shift vectors β and β , respectively of the metrics g and f , in terms of the shift vector q of the geometric mean metric h, This is the most general parametrization of the bimetric N + 1 decomposition [31]. 6 In index notation we have h ij = g ik g −1 f 1/2 k j . 7 For the sake of clarity, we write δ explicitly in every equation where it appears, because it is needed to raise and lower the Lorentz indices.
Since all the spatial bimetric interactions terms, defined in [31] and reported in appendix A.1, depend on the spatial vielbeins, on Λ s and R, their rewriting in the BSSN formalism relies on the conformal decomposition of these variables. The conformal decomposition of the vielbeins read, which tells us that From (3.11) and the conformal decomposition of Λ s and 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.

The covariant extension of the BSSN formulation
As we outlined in section 3.1, the BSSN formulation with the standard gauge is not spatially covariant. As described in [28], 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 [28] do not alter the expressions of the matter sources in the evolution equations, the covariant generalization applies to both metric sectors in BR 8 . As a consequence, the bimetric sources have the same expression as in the BSSN formulation, given by (A.26a), (A.28a) and (A.31a) 9 .
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. [37], 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 W appearing in (2.3) and the separation parameter p appearing in (3.9) in the difference between the shifts (see [31,36] and appendix A.5 for more details). Note that p is also known in the most general β (1) -model [4,55], where the explicit expression for ñ = m −1 pλ −1 in (3.9) is computed 10 .
In the BSSN formulation, the determinant of the conformal metric γ 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, which is the same as before, without the restriction on the determinants ∆ and ∆ , and The tensors A and Â are not traceless, as can be seen by comparing (3.13) with (3.3). The conformal connections in (3.4) are made covariant by introducing two background connections [56]. 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, (3.14) Since the difference between two Christoffel symbols is a tensor, our set of dynamical variables includes tensors only, i.e. γ ij , ϕ ij , φ, ψ, A ij ,Â ij , K,K, Λ,Λ. The standard gauge (3.5) in the cBSSN formalism reads [28], 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 BR, a third metric h is defined 11 . Hence, we can choose the Levi-Civita connection of the conformal spatial metric • χ to define the covariant conformal connection in (3.14). The consequences of this choice are described in more detail in [58].
In the cBSSN formulation, the determinants ∆,∆ of the conformal metrics and the traces A,Â are left undetermined, making the cBSSN evolution equations in (A.45a) and (A.47a) incomplete [28]. We must choose how these quantities evolve, in order to be able to evolve the full system. Following [28], we choose ∂ t A = ∂ tÂ = 0. Regarding the determinants, there are two natural choices, referred to as 'Lagrangian' and 'Eulerian' [28,59,60], given by 11 Actually, in the space of pseudo-Riemannian metrics built on our manifold, we have a path of metrics connecting g and f , corresponding to a geodesic of the trace metric and parameterized by h α = g g −1 f α , α ∈ R (see [31,57]).
with analog expressions for the 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 [58]. The expression for ∂ t Λ i in (3.15c) is explicitly substituted with the Lagrangian formulation of (A.45e).

Well-posedness and gauge choices
As established in [31], the three lapses and shifts of g, f and h are related by, where λ = 1 + p T δp and H is the lapse of h. We are free to choose one lapse function and one shift vector, exactly as in GR; the other four quantities are determined by (4.1). 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 g in (3.15). This gauge can be rewritten in terms of the lapse functions α, H and the shift vectors β , q by using (4.1). Hence, we can impose the standard gauge with respect to g by gauge fixing, say, H and β . 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 + log slicing with respect to χ, and the Γ-driver condition with respect to ϕ. In this case, h would determine the time slicing, whereas 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 [58] for a study of the 'mean gauges', i.e. the gauge choices with respect to the mean metric h.
In GR, the cBSSN formulation is strongly hyperbolic if one chooses the standard gauge and fulfills some other technical conditions [28]. 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 J bi j appearing in the evolution equation for the conformal extrinsic curvature (A.5a) contains the ratio of the lapses W . The general explicit expression of W is not known, but it can be computed explicitly in spherical symmetry ( [36] and the supplementary material to this paper). In that case, 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, W will affect the characteristic matrix in (1.2) and alter the hyperbolic structure of the equations in the g-sector compared with GR. Following the procedure described in [61], 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 W which only depends on three radial derivatives, namely, Here, p 1 is the only nonzero component of p and A 2 ,Â 2 are the A θ θ ,Â θ θ components of the conformal extrinsic curvatures (see appendix A.5 for more details). By using the two momentum constraints (A.54b), (A.54e) and the BCL (A.57) rewritten according to the procedure in [61], these three derivatives can be substituted with expressions involving the dynamical fields only algebraically. In more detail, the BCL can be solved for ∂ r p 1 , and the momentum constraints can be solved for ∂ r (A 2 + K/3) and ∂ r (Â 2 +K/3) 12 . This means that the ratio of the lapses W is a purely algebraic expression in the dynamical fields, and, as such, does not enter the characteristic matrix in (1.2) and does not alter the hyperbolicity of the equations in the g-sector compared with GR. Therefore, if we choose the standard gauge (3.15) for α and β, the equations in the g-sector are strongly hyperbolic.
Note that the algebraicity of W in the dynamical fields is compatible with the fact that, in the case of static spherically symmetric black hole solutions in BR (see [62][63][64] and reference therein) the function describing the ratio of the lapses, τ in [64], is determined by an algebraic relation. This relation corresponds precisely to the PBCL ∂ t C b = 0. Also, this confirms that the bimetric stress-energy tensors (of which the spatial bimetric interactions are the projections [31]) 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 f -sector as well. The equations in the f -sector formally appear the same as in the g-sector, but now the lapse and shift of f are determined by (4.1). As a consequence, the ratio of the lapses appears wherever the lapse and shift of 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 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 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 (4.1) to replace the lapse of f , we could equivalently use it to determine its value on the initial hypersurface, and impose the preservation in time of the PBCL (2.3) by setting its time derivative to 0. This gives an evolution equation for the lapse of f , which becomes a dynamical variable. 13 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 f involves the time derivative of the ratio of the lapses W which nontrivially alters the hyperbolic structure compared with GR.
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, 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 W . Therefore, the computation of 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 g and 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 h. In this case, we see from (4.1) that W appears in the spatial derivatives of the lapses and shifts of both g and f , thus modifying the hyperbolic structure in both sectors. This is investigated in more detail in [58], where both the standard gauge and the maximal slicing for H and q are computed. 12 Note that, in this way, we can freely specify p on the initial hypersurface, and the value of these three derivatives will depend on this choice. 13 More in general, this is an evolution equation for the lapse that we do not gauge fix.

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 field equations, summarized below.
1. In addition to the Hamiltonian and momentum constraints for both metrics, there is an additional bimetric conservation law (the so-called secondary constraint) C b = 0 that, in the free evolution scheme, has to be solved for the initial data on the spacelike hypersurface. 2. The sources in the equations include both the external matter sources and the bimetric sources in (A.5). 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. 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, γ, ϕ or their geometric mean χ. The remaining lapses and shifts are determined by (4.1). [4,35] by imposing the preservation of the bimetric conservation law in time, ∂ t C b = 0, and it is computed explicitly in [36] for spherically symmetric spacetimes in the standard 3 + 1 formulation. The expression in the covariant BSSN formalism is given in the supplementary material of this paper. The ratio between the lapses, 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 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.

The relation between the lapses is established in
Other than these four differences, the system is analogous to the covariant BSSN formulation of the Einstein field equations presented in [28]. 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: 1. The computation of W and p is necessary to be able to solve the bimetric equations in any formulation. At present, W is only computed under the assumption of spherical symmetry, whereas p is computed in spherical symmetry and in the most general β (1) -model; we refer the reader to [4,36,55] for more details. 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 W . Since the latter is known in spherical symmetry, one can study the hyperbolicity of the evolution equations in (A.58). 3. The choice of the gauge is essential in bimetric relativity (as it is in GR as well). In [58], we study some possible gauge choices which alter the hyperbolic structure of the evolution equations in both sectors. In particular, we investigate the gauge conditions with respect to the geometric mean metric h. 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 W and p significantly reduces the accuracy of the integra-tion. We have written a Mathematica/C + + code to perform the simulations, see [65] 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. [37]), 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.

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.

A.1. The BSSN decomposition of the bimetric interactions and sources
Consider the spatial parts γ = e T δe and ϕ = m o T δm o of two Lorentzian metrics g, f . In [54], it is established that the existence of the real square root To be more precise, the freely specifiable spatial vielbein m o is used to compute the vielbein Λ s Rm o such that the spatial part χ of the geometric mean metric h = g#f is given by χ = e T δ(Λ s Rm o ). This is obtained by imposing (A.1c) and solving it for the Euclidean orthogonal transformation R in terms of Λ s and the triangular vielbeins e, m o . Such a solution always exists, as it is part of the polar decomposition of the invertible matrix R o [54,66] (see (3.7)). For the sake of simplicity, we define the new vielbein of ϕ to be m := Rm o rather than Λ s Rm o ; we have the freedom to do that since implying that ϕ is blind to this choice. The interaction terms are not affected as well, since they always contain both Λ s and R, irrespective of this choice. The matrix Λ s explicitly appears in them. On the contrary, R does not appear explicitly, but it is taken into account implicitly inside m. We define the bimetric interactions as [31], where e n (X) are the elementary symmetric polynomials of the linear operator X, and Y n (X) is defined as, See [31] 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 β (n) parameters are d + 1 real dimensionless constants appearing in the bimetric interaction potential, together with the energy scale m [1]. We define the bimetric sources (respectively, the bimetric energy densitites, the bimetric currents and the bimetric spatial stress-energy tensors) as [31], where the summation − m d d n=0 β n is understood in front of all the bimetric sources. Note the relation between the two bimetric currents j b i ,j b i , which implies the relation (A.35) between the momentum constraints.
Here we compute the expressions for the bimetric interactions and sources in the (c)BSSN formalism. We require that the symmetrization condition (A.1) holds for the BSSN variables as well. Since the shifts are the same in the BSSN formalism, we require conditions (A.1a) and (A.1b) to stay the same. The condition (A.1c) should instead lead to its analog in the BSSN formalism, where * Λs, * 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. This naturally leads us to define the vielbeins of φ in the BSSN formalism as, We start by computing the conformal decomposition of the objects in the Lorentz frame. We have, implying, The tetrad m is a set of three Lorentz vectors m a 1 , m a 2 , m a 3 . The operator Λ s a b acts in the Lorentz frame and therefore it is blind to the spacetime index. This means that (A.13) is a set of three equations. Let us denote any of the three Lorentz vectors m a i with u. Then (A.13) becomes, The coefficient in front of p in (A.14) is a Lorentz (and spacetime) scalar, i.e. p and * p are collinear, p = ξ * p . Next, we apply Λ s to (A.13) from the left, Two boosts with collinear boost vectors commute [67, p 50]. This implies that their spatial parts commute as well. Indeed, (A.17) and since p = ξ * p and * p are collinear, the spatial part of (A.17) reads, 18) i.e. the spatial parts commute. Hence (A.15) becomes, and using Λ s u = * Λsu (from (A.14)) in (A.19) we get, It follows that, The case ξ = −1 changes (A.1a) and (A.1b) by introducing a minus sign in front of p. Therefore, we choose ξ = 1 and get 14 , Let us now compute the BSSN version of R o , introduced in (3.7). From (3.11), (A.10) and (A.22c) it follows that, where v is the boost vector and λ the Lorentz factor of the boost in (3.8). Therefore, changing the sign of p is equivalent to consider the boost Λ −1 rather than Λ in the parametrization (3.8), which is an arbitrary choice that one can make. This choice would introduce a minus sign in the relation between the shifts (A.1a) and (A.1b), the latter remaining consistent. Since the two choices are equivalent, and we do not want to change the parametrization in recasting the bimetric standard 3 + 1 decomposition into the bimetric (c)BSSN decomposition, we choose ξ = 1 without losing generality.
We note the following property of the elementary symmetric polynomials, a2 · · · X an] an = f n e n (X), where f is a scalar field. We use this property to compute the bimetric interactions in terms of the BSSN variables, 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 × N identity. Note that both the bimetric interactions and sources can be computed starting with the physical vielbeins and metrics e, m, γ, ϕ, or with the conformal vielbeins and metrics e,m, γ,φ and the conformal factors φ, ψ.

A.2. The bimetric standard N + 1 equations
The bimetric standard N + 1 evolution equations read [31], Lβ, β and β being the shift vectors of the two metrics. The bimetric standard N + 1 constraint equations (BCE) are, The effective sources are the sum of the bimetric sources given by (A.5), and the external matter sources, Note that the relation between the two bimetric currents in (A.5), implies that the two momentum constraints are related to each other 15 , The bimetric conservation law (BCL), in its asymmetric and symmetric form, reads (A.36b)

A.3. The bimetric BSSN equations
We now write down the bimetric BSSN equations, applying the procedure in [27, section 11.5] to (A.32), (A.33) and (A.36), using (A.31) for the bimetric sources. The bimetric BSSN 3 + 1 evolution equations for the g-sector are, and for the f -sector The bimetric BSSN constraint equations are Now we compute the BSSN decomposition of the bimetric conservation law (A.36). We note that, due to (3.1), the relation between the two conformally related covariant derivatives of a vector field X i is, with analogous formulas for the f -sector. Using the definition of the conformal extrinsic curvatures (3.13) in the cBSSN formalism, the bimetric conservation law becomes Finally, we insert the expressions computed in appendix A.1 for the quantities appearing in (A.44), which are written in terms of the BSSN variables.

A.4. The bimetric covariant BSSN equations
At this point, we apply the method described in [28] 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, where D B and R B are the covariant derivative and the Riemann tensor of the background geometry, and Note that all Lie derivatives in (A.46) are Lie derivatives of tensors, not tensor densities as in (A.38). Equation (A.46e) is an identity for the Ricci tensor of γ ij in terms of the background geometry [28,68] (see appendix A of [68] for the proof). The superscript TF means 'tracefree'. Note that all the traces are with respect to γ except the trace of J eff , which is with respect to γ. For the f -sector we have, where the same conventions and notations as for the g-sector are implied. It holds, The bimetric cBSSN constraint equations are The cBSSN BCL is the same as in one in (A.44), since the bimetric interactions (A.26) and (A.28) are not changed, and the bimetric sources are the same as in (A.31).

A.5. The bimetric covariant BSSN equations in spherical symmetry
The equations in this appendix were computed with the Mathematica package bimEX [69], which can compute the bimetric cBSSN equations for any desired ansatz. We write down the cBSSN equations in spherical symmetry, assuming the following ansatz, The background geometries for both γ and ϕ are chosen to be the spatial part of the flat metric in spherical coordinates 16 , From now on, we will assume the time and radial dependence of all the fields. We define the function, 52) and the 'shifted elementary symmetric polynomials' [18], where the β (n) are five real dimensionless constants appearing in the bimetric interaction potential [1]. 16 Note that this is not the Minkowski metric in the Lorentz frame, whose spatial part remains δ = diag [1, 1, 1].
Consistency conditions between the various equations imply that q θ = p 2 = 0, and q φ , p 3 do not appear explicitly into the equations, so we can set them to zero without losing generality. Also, R is the identity and m =m o , which simplifies the computations considerably.
The constraint equations (A.49a) read, ã´− κ f j mr +j br = 0, where the bimetric energy density and currents are given by We can solve p 1 from both (A.54b) and (A.54e), obtaining respectively (A.56b) The asymmetric-and simpler-version of the BCL (A.44a) is, (A.57) The evolution equations (A.45) and (A.47), modified to get the evolution of the components A i j ,Â i j rather than A ij ,Â ij , reduce to 2b∂ rb a 2 r +b 2Λr r , (A.59d) and J br r =¨Û R The evolution equation for p 1 is obtained from the cBSSN version of (A.8) in [31], (A. 62) and appropriately using the Hamiltonian constraints in the evolution equations for the traces K and K , and the momentum constraints in the evolution equations for the conformal factors φ, ψ.