Geometric mean of bimetric spacetimes

We use the geometric mean to parametrize metrics in the Hassan-Rosen ghost-free bimetric theory and pose the initial-value problem. The geometric mean of two positive definite symmetric matrices is a well-established mathematical notion which can be, under certain conditions, extended to quadratic forms having the Lorentzian signature, say metrics $g$ and $f$. In such a case, the null cone of the geometric mean metric $h$ is in the middle of the null cones of $g$ and $f$ appearing as a geometric average of a bimetric spacetime. The parametrization based on $h$ ensures the reality of the square root in the ghost-free bimetric interaction potential. Subsequently, we derive the standard $n+1$ decomposition in a frame adapted to the geometric mean and state the initial-value problem, that is, the evolution equations, the constraints, and the preservation of the constraints equation.


Introduction
We use the geometric mean to parametrize metrics in the Hassan-Rosen (HR) ghost-free bimetric theory and pose the initial-value problem. The HR bimetric theory [1][2][3][4] is a nonlinear theory of two interacting spin-2 classical fields which is free of instabilities such as the Boulware-Deser ghost [5]. The HR theory is related to de Rham-Gabadadze-Tolley massive gravity [6][7][8]. Recent reviews of these theories can be found in [9,10]. The geometric mean of two positive definite symmetric matrices can be written [11], (1.1) Here M 1/2 denotes the principal square root of a square matrix M with no eigenvalues on the negative real axis R − , that is, M 1/2 is the unique solution X of the matrix equation X 2 = M whose eigenvalues lie in the right half plane [12]. The notion of the geometric mean can be extended to symmetric bilinear forms with the Lorentzian signature under certain conditions. Given two spacetime metrics g and f , we can define their geometric mean as [3], provided that the principal square root S := (g −1 f ) 1/2 exists. According to the theorem from [3], the principal square root S exists if and only if there exists a common timelike direction and a common spacelike hypersurface element relative to both g and f . In this context, the geometric mean has a nice geometrical interpretation: The null cone of the mean metric h = g # f = gS will be in the middle of the null cones of g and f .
The shapes of g and f are given by lapses and spatial metrics.
The shifts of g and f are given by the separation p and the position relative to the shift q of h.
The lapse of h depends on p (a smaller p corresponds to a wider null cone of h).
t Σ g f h p q Consequently, if a solution to the HR bimetric field equations exists, there will always be a geometric mean metric h. This makes h useful when adapting the space-plus-time foliation. Moreover, it can be used to parametrize possible real square root realizations of metrics with respect to h. Let us consider a foliation adapted to h. We first specify the shift vector q of h, as illustrated in Figure 1. Next, we split the degrees of freedom of g and f across g, f and h. We parametrize the volume contributing parts of g and f by their lapse functions and the spatial metrics. Finally, in addition to the slice dependent shift q, we give the separation p between the null cones of g and f also relative to h. This gives all possible configurations of g and f which have the real principal square root.
After the parametrization, we decompose the bimetric field equations and state the initial-value problem. The space-plus-time decomposition is done for an arbitrary spacetime of dimension d = N +1 using the prescription from [13]. The projection yields the constraint equations and the evolution equations in standard N +1 form. Furthermore, the projection of the contracted Bianchi identities yields an additional equation that preserves the ghostfree bimetric potential when in N +1 form, ensuring the correct number of d(d − 2) − 1 propagating degrees of freedom.
The derived N +1 form of the bimetric field equations resembles two copies of General relativity (GR). For comparison, the 3+1 split is demonstrated for the bimetric spacetimes where the two sectors share a common spherically symmetry. The N +1 equations can be used as a starting point for further analysis, for instance, in numerical bimetric relativity.
Structure of the paper. In the rest of this section we review some technical properties of the field equations in GR and the HR bimetric theory. In section 2, we parametrize the metrics using the geometric mean, and decompose the effective stress-energy tensors of the ghost-free bimetric potential. After projecting the contracted Bianchi identities of the bimetric effective stress-energy tensor, we obtain the N +1 form of the bimetric conservation law (the so-called secondary constraint obtained using the Hamiltonian formalism [2,4]). In section 3, we give the complete set of the HR field equations in standard N +1 form, also summarized in Box 1 on page 18. In the case of spherically symmetry, the reduced N +1 equations are studied in section 4. The paper ends with a short summary and discussion.
Notation. The tilde indicates the most of the variables in the f -sector. The preprint of this paper with the color-highlighted variables can be found as an ancillary file on arXiv.

The N +1 formalism in GR
The Einstein-Hilbert action for a metric g is given as, where R g denotes the Ricci scalar for g, and L m g denotes the Lagrangian density of matter fields minimally coupled to gravity. The parameter κ g is Einstein's gravitational constant. Varying (1.3) with respect to g gives the Einstein field equations (here in operator form), (1.4) where G g and T g denote the Einstein tensor and the stress-energy tensor, respectively, Both G g µ ν and T g µ ν are symmetric (self-adjoint) with respect to g. The Einstein tensor satisfies the contracted Bianchi identity ∇ µ G g µ ν = 0 which holds for any g. Due to the contracted Bianchi identity, the matter fields obey the conservation laws ∇ µ T g µ ν = 0. The kinematical and dynamical parts of the metric field g can be separated using the N +1 decomposition; here we follow the prescription from York [13]. One starts by foliating the spacetime into a family of spacelike hypersurfaces {Σ} with the future-pointing timelike unit normal n on {Σ} satisfying n µ n µ = −1 with respect to g. The geometrical objects are then projected onto the slices using the operator, The metric induced on the spatial slices reads, Adapting a suitable chart 8) where N and N i are respectively the standard lapse function and shift vector associated with the foliation. In this chart, the metric decomposes according to (1.7) as, The lapse N is a strictly positive scalar field (enforcing n 0 > 0), the shift N i is a purely spatial vector field, and γ ij is a spatial Riemannian metric (the shift vector N i is confined to the spatial hypersurface, and its indices are lowered using γ ij ). Also, the inverse of γ is obtained through γ ik γ kj = δ i j . In matrix notation, the metric (1.9) reads, (1.10) Shear operator. Given a spatial vector ξ, one can define the shear operator, (1.11) Note that the shear is unimodular, det Ξ[ξ] = 1, and it forms an Abelian group isomorphic Taking the shift vector N as a shear factor, the metric g can be expressed, (1.12) The inverse of g then simply follows from the properties of Ξ, (1.13) Moving frames. Two types of spacetime frames are used in this work. One is a Cauchy adapted frame [14] based on a coframe {θ µ } with a dual { θ µ }, where θ µ ( θ ν ) = δ µ ν and, Note that a Cauchy adapted frame is seemingly a 'coordinate frame' where the spacetime indices are abused to label the basis vectors. Recognizing θ µ ν = Ξ[ N ] µ ν from (1.14a), or equivalently θ = Ξ[ N ] in matrix notation, the expression (1.14c) is the same as (1.12).
The other type is an orthonormal frame or a vielbein, Here we used the beginning Latin letters to denote the indices in the local Lorentz frame.
Observe that e i a denotes e −1 in matrix notation, and the vector field E0 is the timelike unit normal n. The spatial metric is accordingly decomposed as γ ij = δ ab e a i e b j (or γ = e Tδ e) using a vielbein e a = e a i dy i tangential to Σ. Subsequently, the metric (1.12) becomes, 16) or g = E T ηE in terms of a vielbein E in the lower triangular form, (1.17) Beware that an arbitrary vielbein can be triangularized by a local Lorentz transformation if and only if the apparent lapse for the metric is real in a given coordinate chart (see Lemma 2 in appendix B). For a single metric, a general coordinate transformation can be used to ensure that the lapse is real. This might not be possible for two arbitrary metrics.
Projection meta-operators. Any symmetric tensor field X µν can be decomposed as, is the mixed projection, and (1.19b) The trace of X can be expressed by, . (1.20) For the stress-energy tensor T µν , the physical interpretation of the projections is following: Here L n denotes the Lie derivative along the vector field n. The perpendicular and the mixed projection of the Einstein tensor reads, where R = γ ij R ij is the trace of the Ricci tensor R ij defined by the spatial covariant derivatives D i compatible with γ ij . Combining (1.22) with the respective projections of the stress-energy tensor yields the constraint equations, Equation (1.23a) is called the scalar or Hamiltonian constraint while (1.23b) is called the vector or momentum constraint; they must be satisfied on the initial hypersurface.
To end up with a more robust form of the evolution equations, we decompose the field equations (1.4) based on the Ricci tensor (such a procedure is due to York [13]), The full spatial projection of the d-dimensional Ricci tensor gives, Combining (1.25) with the projection of the stress-energy tensor, then expanding the Lie derivatives in terms of n = N −1 ∂ t − N i ∂ i , we get the evolution equations, where the trace of T g was decomposed using (1.20). The equations (1.26) will be hereinafter referred to be in standard N +1 form.

Bimetric field equations
We first consider a general class of bimetric actions consisting of two Einstein-Hilbert terms S EH [g] and S EH [f ] each coupled to separate matter fields, and the interaction term L int g,f that depends on the scalar invariants of the operator g µρ f ρν (written g −1 f in matrix notation). 1 The bimetric action of such a class reads (in an arbitrary dimension d), where κ g and κ f are the gravitational constants for two sectors. By varying the action with respect to g and f , one obtains the field equations, where T g and T f are the stress-energy tensors of the matter fields, while V g and V f are the effective stress-energy tensors of the bimetric potential L int g,f (g −1 f ), The dependence of L int g,f only on g −1 f , combined with the definitions of V g and V f from (1.29), implies the following algebraic identities, and the differential identity, where ∇ µ and ∇ µ are the covariant derivatives compatible with g and f , respectively. The identity (1.30b) was first proved in [16]. The identity (1.30c) was shown in [15] assuming that the action (1.27) is invariant under the diagonal subgroup of diffeomorphisms. In particular, it can be proved that the differential identity (1.30c) is a consequence of the definitions (1.29) provided that L int g,f is a scalar function of g −1 f . Assuming that ∇ µ T g µ ν = 0 and ∇ µ T f µ ν = 0, the equations of motion (1.28) imply the Bianchi constraints, which are not independent according to (1.30c). The equations (1.31) will be hereinafter referred to as the bimetric conservation law.

The ghost-free bimetric potential
The absence of ghosts in the Hassan-Rosen bimetric theory is ensured by the interaction potential of the following form [1], where S is the principal square root of g −1 f , denoted as S = (g −1 f ) 1/2 . In particular, the principal branch of the square root provides an unambiguous definition of the HR theory that guaranties the existence of a spacetime interpretation [3]. The sign of the interaction term is parametrized by a label. All dependent signs will be indicated by the same label in the rest of this work. The negative sign convention for the interaction term places the bimetric stress-energy contributions V g and V f on the gravitational (left) side of the field equations, while the positive sign places V g and V f on the matter (right) side.
The bimetric potential is parametrized by a set of real constants {β n }. The mass scale is given by m, making the β-parameters dimensionless. The coefficients e n (S) in (1.32) are the elementary symmetric polynomials, which are the scalar invariants of S obtained through the generating function [17], (1.33) Note that e n (X) = 0 for all n above the rank of X due to the Cayley-Hamilton theorem. This implies that the summation in (1.32) stops at n = d for the nonsingular S.
Elementary symmetric polynomials. We summarize some important properties of the elementary symmetric polynomials that are used throughout this work. Let X a b be an invertible operator on a vector space of an arbitrary dimension D. From (1.33) we have, (1.34) A particularly useful expression for the elementary symmetric polynomials reads, an . The expansion of the antisymmetrizer gives, where we introduced, The function Y n satisfies the recursive relation (which can be used as a definition of Y n ), The Cayley-Hamilton theorem can be written Y D (X) = 0. From (1.38) and (1.34) follows, Using (1.36) one obtains the following identities involving a vector u and covector ω (these relations occur in the decompositions of Lorentz transformations and square roots), Moreover, the recursive relation (1.38) implies the following identities for traces, which are Newton's identities in disguise, Applying the derivative δ Tr(X n ) = n Tr(X n−1 δX) on (1.41) gives, Hence, Y n (X) can be equally defined as the derivative of e n+1 (X), The HR field equations. The HR bimetric field equations have the form (1.28) with the following effective stress-energy contributions of the ghost-free bimetric potential (1.32), The function Y n (S) is defined in (1.37), (1.38), or (1.44). The two contributions in (1.45) are not independent because they satisfy the identities (1.39) and (1.30).

Usage of the geometric mean
In this section we parametrize the metrics relative to the geometric mean and obtain the projections of the bimetric stress-energy tensors V g and V f as determined by the respective Eulerian observers. For this we utilize two properties of the bimetric potential.
(i) The class of bimetric actions (1.27) can exhibit a duality in vacuum under the interchange g ↔ f and κ g ↔ κ f . In the HR theory, the duality is explicit under the additional exchange β n ↔ β d−n as a consequence of (1.34). The practical use is that we need to derive only the relations for one sector and then simply apply the results to the other, provided that the symmetry g ↔ f is not algebraically broken (for example, when doing the N +1 decomposition). One way to achieve this is by using vielbeins.
(ii) The function space of the generating function (1.33) is spanned by the polynomial basis {t n }. In light of this, the bimetric potential (1.32) and its stress-energy contributions (1.45) can be seen as a linear combination of the span {β n }, where the elementary symmetric polynomials e n (S) and their derivatives Y n (S) act as the coefficients. Consequently, the apparent basis ±m d n β n can be neglected in the intermediate calculations; namely, we can assume without loss of generality that, and recover m d n β n (together with the chosen sign) at the end when stating the complete set of equations. Note also that the summation range is omitted from m d n β n in (1.32) and (1.45) since the vanishing e n (S) and Y n (S) will truncate the sum by killing all unnecessary terms depending on the rank of S; in particular, e n>d (S) = 0 and Y n>d−1 (S) = 0. The summation range will self-adapt to the lower dimension of the spatial hypersurface upon the projection of the field equations using the N +1 decomposition,

Parametrization
We work in a Cauchy frame θ = Ξ[q] adapted to h that is parametrized by a shift vector q. The shape of the metric g is given by a lapse function N and a spatial vielbein e that defines the spatial metric γ = e Tδ e. The metric f is similarly parametrized by a lapse function M and a spatial vielbeinm that defines the spatial metric ϕ =m Tδm . The lapses and the spatial metrics are respectively perpendicular and tangential projections on the spacelike hypersurfaces {Σ}. The mixed projections are given by different shift vectors which define the respective Cauchy frames. The overview of the parametrization as a dependency graph is shown in Figure 2 (which also gives the evaluation order). The Cauchy frames adapted to h, g and f reads, The separation between the frames is be parametrized by a boost vector p a in the local Lorentz frame of h, giving the Lorentz factor, Moreover, given an arbitrary p andm, the vielbein m in (2.2) is obtained as, whereΛ is the spatial part of the Lorentz boost given by, The relation (2.4) ensures that the spatial metric χ of h is symmetric, An important remark is that ϕ = m Tδ m =m Tδm sinceR is an orthogonal transformation R TδR =δ. Note also that λ is an eigenvalue ofΛ with the eigenvector v := pλ −1 , where v Tδ v T < 1. All other eigenvalues ofΛ are 1, and detΛ = λ.
The resulting parametrization of the metrics reads, where the lapse of h is given by, Notice that whenever the Cauchy frames coincide (λ = 1,Λ =Î), we have, A detailed proof is given in appendix B. In brief, the effective shift vectors of g and f are, Subsequently, the parametrization is exhaustive based on the results from section 3 in [18], and valid based on the theorem from [3]. The resulting square root Notice how the algebraic form of S is naturally given in the Cauchy frame Ξ[q] of h. The square root can be further decomposed as, using the similarity transformation, This form is useful when calculating matrix functions of S since Z cancels out in traces and powers of S (for instance, when expanding the elementary symmetric polynomials e n or their derivatives Y n ).

Different views of the Lorentz frame.
The boost vector v can be given with respect to γ or ϕ by changing the frame (basis) using the spatial vielbeins, We can equally recastΛ 2 =Î + λ 2 v v Tδ by the corresponding similarity transformations, , the eigenvectors of Q and Q are (2.15) having the eigenvalue λ 2 , To decompose the square roots of γ −1 ϕ and ϕ −1 γ, we define, so that γ −1 ϕ and ϕ −1 γ can be symmetrically split into two components, Then, the spatial metric of h reads, Note that χ is symmetric, χ = χ T . Hence, any function of D or B is symmetric with respect to γ, and any function of D and B is symmetric with respect to ϕ. Finally, using (2.18), the boost vectors n and n can be related across the sectors,

Decomposition of the potential
The scalar density √ −g V can be written, Here we used (1.40) together with the following relations forΛ, where p := p Tδ and v := v Tδ denote the adjoint of of p and v, respectively. Since λ −1 e n−1 (B) = e d−n (D) det(me −1 ), we can rewrite (2.22) as,

The bimetric stress-energy tensors
Here we evaluate the projections of the bimetric stress-energy tensors V g and V f as seen by their Eulerian observers. The projection operators for the frames of h, g, and f are, We first determine the projections of V g = Y n (S) relative to the Eulerian observer of g.
The one-side projection of S can be evaluated as, where the double brackets are used to enclose matrix expressions when putting indices and, The composite symbols for the mixed terms (Q U) and ( QU) are defined as, where we used Q D = D Q from (2.16) and (2.18). Observe that γ ik U k j , γ ik (Q U) k j , ϕ ik U k j , and ϕ ik ( QU) k j are symmetric tensor fields. The projections of V f can be deduced from the duality between the sectors by using (1.34) and (1.39); we obtain, Having V = e n ( D), it holds, where U and U satisfy Un = U n and n T γ(Q U)n = λ 2 n T ϕ U n.
We summarize the projections of V g and V f , now given in a slightly different form, Recall that the projections are implicitly spanned by ±m d n β n . For completeness, the components of V g and V f in terms of (2.25) and (2.29) are given in appendix C.

The bimetric conservation law
In GR, the conservation law ∇ µ T g µ ν = 0 needs to be specifically imposed upon the N +1 decomposition, otherwise the constraint equations will fail to be preserved during the evolution [19]. For a class of bimetric theories having the field equations (1.28), the effective stress-energy tensor contributions (1.29) are included in the Bianchi constraints (1.31). Therefore, after putting the bimetric field equations into the N +1 form, we also need to assume that the bimetric conservation law ∇ µ V g µ ν = 0 is specifically satisfied. In the following, we derive the projection of the conservation law for the ghost-free bimetric potential V = e n (S). Let V g be given in terms of the energy density ρ, the energy flux covector j i , and the stress tensor J i j , all of them measured as (2.28) by the Eulerian observer relative to g. A straightforward manipulation yields the following projections of ∇ µ V g µ ν = 0 (summarized in appendix A), where the perpendicular projection of ∇ µ V g µ ν with respect to h can be identified as, The energy density and the energy flux of the bimetric potential satisfy, Moreover, for the Lie flow along an arbitrary vector field ξ i , it holds, The proof of the lemma is given in the first part of appendix D. Useful choices for ξ i are the shifts N i and M i , or the boost parameters n i and n i , where L n n = L n n = 0.
Applying Lemma 1 on (2.33) gives the N +1 form of the conservation law for V g (the derivations are in the second part of appendix D), A more symmetric form that reflects the duality g ↔ f is given by (D.17), copied here, Note that the covariant derivatives in (2.38) can be easily converted to the partial derivatives of the respective metric since γU and ϕ U are symmetric.

The HR bimetric field equations in standard N +1 form
The HR field equations (1.28) are two copies of (1.4) with the additional V g and V f that couple two sectors through the identities (1.30). Therefore, to state the HR field equations in the N +1 form, we only need to include the effective stress-energy contributions when projecting the particular sector (using an appropriate timelike unit normal): where T g = 0 and T f = 0 in the "bimetric vacuum." The substitution is straightforward (it is mostly a "copy & paste" of the GR results). The end result is the complete set of the HR equations in standard N +1 form, given below. The evolution equations for the dynamical pair (γ ij , K ij ) in the g sector reads, while the evolution equations for the dynamical pair (ϕ ij , K ij ) in the f sector reads, The full spatial projections of V f + T f and V g + T g are given by, The traces of V g + T g and V f + T f can be obtained from (1.42), The constraint equations for the two sectors are given by, The constraints for the mixed projections (3.5c)-(3.5d) are coupled by (3.6c) yielding, The last equation can be used as a replacement for either (3.5c) or (3.5d), where the other becomes the equation of motion for p (i.e., for n and n). Observe that the lapses N , M , and the shift vector q of h are absent from the constraints, as expected. In addition to (3.2)-(3.7), we also have to assume specifically that (i) the conservation laws ∇ µ T g µ ν = 0 and ∇ µ T f µ ν = 0 hold, and (ii) the conservation law for the bimetric potential ∇ µ V g µ ν = 0 holds. Failing to do so, the constraints will not be preserved during the dynamical evolution by (3.2). Note that we do not need to regard ∇ µ V f µ ν = 0 because of the differential identity (1.30c). Consequently, the propagation of the constraints requires that the equation (2.37) is satisfied, The equation (3.8) is equivalent to (and algebraically the same as) the so-called secondary constraint which is obtained using the Hamiltonian formalism [2,4]. In the context of the N +1 decomposition employed here, it is rather seen as the conservation law for the ghost-free bimetric potential. This concludes the HR bimetric field equations in standard N +1 form. An overview of the used variables is given in Table 1. A handy overview of the N +1 equations is given in Box 1. The counting of the degrees of freedom is given in Table 2. The counting begins with 2d(d − 1) equations for the dynamical fields (γ ij , K ij ) and (ϕ ij , K ij ). Besides, there are also d + 1 constraint equations, the bimetric conservation law equation, and an overall gauge freedom which allows us to remove additionally d of the dynamical fields. As a result, we are left with d(d − 2) − 1 truly dynamical conjugate pairs, which is the number of the propagating degrees of freedom in the HR theory.

Constraint equations
Box 1. The HR vacuum field equations in standard N +1 form.

Dictionary
In conclusion to this section, we address one important condition that the projected bimetric conservation law (3.8) itself stays preserved under the dynamical evolution (i.e., that its time derivative vanishes identically). Such a requirement imposes a relation between the two lapses as shown in [4], which gives the correct number of kinematical variables related to the gauge freedom. To illustrate this, consider a β 1 -model where all the β-parameters vanish except β 1 ; then, U = λ −1Î , U = D, and the preservation of (3.8) requires, The time derivative in (3.9) can be resolved using the evolution equations. For example, we can determine N in terms of the other fields after setting q = 0 and M = 1. By a more general analysis, it can be shown that necessarily M/N = W , where W is a spatial scalar field determined from the other fields [4]. 2 Explicit examples of tying up M/N in spherical symmetry will be tackled in section 4. Note that the preservation of the equation (3.8) is sometimes referred to as the stability condition for the secondary constraint in the context of the Hamiltonian formalism [4,20,21]. However, the preservation and the stable propagation of the constraints are usually considered as two distinct concepts [22], where one assumes that a quantity propagates in a stable manner if it depends continuously on its initial values [23]. Such a treatment will be studied elsewhere.

Spherically symmetric spacetimes
We now investigate the special case where the two metric sectors share the same spherical symmetry. One reason is that the bimetric field equations in spherical symmetry reduce to 1+1 dimensions and are simpler to pose and solve. For instance, all the shift vectors become scalars and we do not need to explicitly determine the spatial rotationR from (2.4). Also, the reduced 1+1 equations can serve as a toy model for experimenting with different gauge choices and slicings. As a physical motivation, many important aspects of gravitational collapse can be studied in spherical symmetry, for example, black hole or structure formation and growth. Finally, when treated numerically, the spherical equations can be solved faster, and with higher accuracy [24]. Hence, treating problems in spherical symmetry seems as a good starting point for the bimetric numerical relativity. Let d = 4. The general form of the metrics in spherical polar coordinates reads [25], where α and α are the lapse functions, and (A, B, A, B) denote the nontrivial components of the spatial vielbeins such that γ rr := A 2 , γ θθ := B 2 , ϕ rr := A 2 , and ϕ θθ := B 2 . The radial components of the shifts are parametrized by the mean shift q and the radial separation p, In addition, we have the components of the extrinsic curvature, All these variables are functions of (t, r) to be solved for in general. For the following, we introduce R := B/B and define a partial span by {β n }, The nonzero components of the projections of the bimetric stress-energy tensor V g are,  Table 1 are given for the case of spherical symmetry in appendix F. The scalar constraints (3.5a)-(3.5b) are, The vector constraints (3.5c)-(3.5d) are, The last two equations can be recombined using (3.6c), The radial shift separation can be determined from (4.7d) as, The projection (3.8) of the bimetric conservation law reads, where p can be eliminated using (4.8). The evolution equations for the spatial metrics are, The evolution equations for the extrinsic curvatures are, The algebraic symmetry between the two sectors is obvious. The set of 3+1 equations in each bimetric sector is directly comparable to GR [24,25]. As a sanity check, let us decouple g and f by setting m 4 β 0 = Λ g , m 4 β 4 = Λ f , and β 1 = β 2 = β 3 = 0. We get the stress-energy components for the cosmological constants, Here we also have U = 0 and U = 0, so the bimetric conservation law becomes automatically satisfied. Notice the implications of choosing the V -sign convention (inside R n k ) on the bimetric stress-energy projections (4.5), (4.6), and (4.12).
As another sanity check, we can compare the 3+1 split to the existing cosmological and black hole solutions exhibiting spherical symmetry in the HR bimetric theory. We start with the derivation of cosmological solutions [26] based on the Friedmann-Lemaître-Robertson-Walker ansatz. Therein, the lapse of f is denoted as X, while the lapse of g is fixed to one, so the ratio of the two lapses reads W = α/α = X. The constraint equations reduce to the equation (2.12) in [26]. Moreover, the equation (2.13) fixes the ratio between the lapses and corresponds to the preservation of the bimetric conservation law condition. On the other hand, for a comparison with the black hole solutions, let us consider a choice of variables in the "radial gauge" B = r where we take R = B/B as a primary field. Furthermore, we express the equations in terms of two additional scalar functions τ = α/α and Σ = A/A (which appear in the eigenvalues of the square root). These names are complaint with the spherically symmetric ansatz from [27]. The constraints reduce to the fifth equation (23e) in [27]. The algebraic condition for τ (25a) in [27] fixes τ relative to the other fields, and roughly corresponds to the preservation of the bimetric conservation law condition imposing the ratio between the two lapses W = τ = α/α. This shows the usefulness of the space-plus-time dissection of the field equations in analyzing the role of the fields chosen in a specific ansatz.
Note that the variables can be differently reparametrized to simplify the equations. For example, one can use the exponential functions instead of A, B, A, B, or introduce the densitized lapses α = αA −1 and α = α A −1 which nicely appear as the factors in the shift separations. Also, one option is to do a spatial gauge fixing and work in a zero shift q = 0. The next step would be to study the required conditions for the gauge choice in general. For example, the lapse function can be determined from the maximal slicing which is done with respect to the geometric mean metric, functioning as a common "singularity avoiding" slicing condition for both sectors. However, there are many other choices for gauge fixing that may give the condition for the lapse function (the time slicing condition), and the condition for the shift vector (the spatial gauge condition). Posing such conditions and then solving the spherical equations will be treated elsewhere.

Degrees of freedom in spherical symmetry.
Here we determine the number of truly dynamical variables in the case of spherical symmetry. In GR, the spherical reduction of the number of degrees of freedom is reflected in the fact that we can obtain a fully constrained system with no evolution equations, which is compliant with Birkhoff's theorem. On the other hand, the constraint equations in the HR theory cannot remove all of the dynamical degrees of freedom after imposing the spherical symmetry. In particular, the spherical symmetry reduction does not suppress the propagating degrees of freedom associated with the monopole radiation (contrary to Birkhoff's theorem). To see this, let us start with the field inventory (p, q, α, A, B, K 1 , K 2 , β, A, B, K 1 , K 2 ). We also have eight evolution equations in (4.10) and (4.11), four constraint equations in (4.7), one conservation of the bimetric potential equation (4.9), and one equation relating the two lapses. The gauge freedom can be used to fix α and q eliminating one conjugate pair. One of the constraint equations can be used to determine p (4.8). At the end, we are left with one dynamical conjugate pair governing radial fluctuations in time.

Summary and outlook
Various averaging operations on positive definite matrices have been of interest to operator theorists, physicists, and statisticians for a long time [28]. Particularly intriguing has been the notion of the geometric mean [28][29][30]. The geometric mean on the space of Lorentzian matrices was earlier not addressed in the mathematical literature. In the context of bimetric theory, such a geometric mean first appeared in [3]. Now, let us consider a space P(u, Σ) of all Lorentzian metrics for which a vector field u is future-pointing timelike and a family of nonintersecting hypersurfaces {Σ} is spacelike, as illustrated in Figure 3a. We assume that the spacelike slices {Σ} locally arise as the level surfaces of a smooth scalar field τ called the time function. In particular, the foliation is set out by a closed one-form Ω, and since Ω is closed, there locally exists a scalar field τ such that Ω = dτ , where we also assume u µ Ω µ = 0 so that u is never tangent to {Σ}.
The space P(u, Σ) can be endowed with a Riemannian metric ∆ measuring the distance between two metrics g 1 and g 2 (which are points in P) defined as, where σ[X] denotes the set of the eigenvalues (spectra) of X. Strictly, the metric ∆ can be defined only on a space of positive definite symmetric matrices [31]. 3 Here we extend such a definition to P(u, Σ) since the conditions for the existence of the principal logarithm are the same as for the existence of the principal square root [12]. Consequently, the validity of the definition (5.1) comes from the theorem from [3]. This makes applicable the considerations from [29][30][31][32]. The set of spaces P(u, Σ) can be identified as a natural habitat of the solutions in the HR bimetric theory. An open problem is still a definition of the geometric mean of several Lorentzian metrics. An important observation is that any two points g, f in P(u, Σ) can be joined by a unique geodesic for which a natural parametrization is given by, The geodesic can be illustrated by a null cone sweep, as shown in Figure 3b. In [33] it was shown that, if any two metrics on a geodesic (5.2) share the same Killing vector field, then all the metrics on the geodesic share the same isometry. The geometric mean (1.2) is obviously the midpoint of the geodesic (5.2), which justifies the expression that the null cone of h is 'in the middle' of the null cones of g and f . The midpoint is unique since The geometric mean metric can be used, in principle, to measure a mean curvature when adapting a singularity avoiding slicing. Namely, having GR as the guideline, a simple choice of setting the time and the spatial slicing conditions will not work. The selection of the gauge conditions needs to be physically and geometrically motivated so that the lapse and shift are adapted to be suitable for the corresponding spacetime. In this sense, the parametrization based on the geometric mean is usable in gauge fixing. Note that any dynamical choice of the lapse and shift as functions of the metric or first-order derivatives of the metric is allowed, provided that the evolution of the bimetric constraints is well-posed. Ensuring the stable propagation of the constraints is important since, if the constraint evolution equations are not well-posed, the unphysical modes (which are normally suppressed by the constraint equations) will not be bounded but propagated as amplified by the free evolution. The causal propagation of the bimetric constraints will be studied elsewhere.
The N +1 form of bimetric equations derived here is compatible with York's standard version that is used in GR. The applied procedure is complementary to the one using the Hamiltonian formalism [4]. The variables in the Hamiltonian formalism are evaluated by varying the N +1 form of the HR action, while the variables in this work are derived by N +1 projecting of the bimetric field equations obtained from the HR action varied in a general form. As shown, the N +1 projection was straightforward because of the duality between the two metrics in the HR theory. Both procedures yield the same form of the secondary constraint (for a detailed comparison of the variables see appendix E).
Finally, the 3+1 form of the equations for the case of the spherical symmetric reduction was given in section 4. One of the applications would be in a systematic study of gravitational collapse in spherical symmetry (e.g., featuring black hole or structure formation and growth). In such a case, however, further adapting of the evolution equations to a scheme like the BSSN formulation [34,35] might be needed for achieving numerical stability. . The operator X can be, for instance, T g , G g − κ g T g , or V g . We summarize the derivations [19] inspired by the paper of Frittelli [22] (similar expressions can also be found in [25]). A direct expansion yields, The projection of ∇ µ X µ ν = 0 along n ν reads, which can be simplified to, Here we used ∇ µ n ν = + K µν + n µ D ν log N and, Also, note that for an arbitrary tangential vector ξ holds Similarly, the projection ∇ µ X µ ν = 0 using ⊥ ν α reads, which can be simplified to,

Appendix B Proof of Proposition 1
Here we combine the results from the sections 3 and 5 of [18], the theorem from [3], and one additional lemma (given below). If S is an arbitrary congruence, f = S T gS, then the symmetrization condition h = gS = h T is equivalent to the requirement that the congruence S is symmetric (self-adjoint), S = S where S := g −1 S T g, which is further equivalent to the fact that S is a square root, f = gS 2 . In terms of vielbeins, where g = E(g) T ηE(g) and f = E(f ) T Λ T ηΛE(f ), the symmetrization condition reads [36], where Λ is a residual overall local Lorentz transformation (LLT). The relation between the square root S and the Lorentz transformation Λ is just the change of basis S = E(g) −1 ηΛE(f ). Before proceeding, we address one important issue.

Lemma 2. An arbitrary vielbein can be triangularized by a local Lorentz transformation if and only if the apparent lapse of the associated metric is real in a given coordinate chart.
Proof. Consider an arbitrary vielbein E that we want to triangularize by a LLT Λ, The equation for the upper right component reads e T 1 + p R e 3 = 0, which becomes v R e 3 = −e T 1 for v = pλ −1 . The vector v can be determined, However, the Lorentz factor must satisfy λ ≥ 1, which gives the condition on the vielbein, Now, assuming a common spacelike hypersurface for both the metrics (which exists if and only if the principal square root exists [3]), we can triangularize the vielbeins in (B.1) using the above lemma. Then, starting from the vielbeins in the lower triangular form, which can be expanded as (noteΛpλ −1 = p),

Appendix C Components of V g and V f
The components of V g can be determined using (1.18), Note that (1.30b) can be expanded, Using (2.36b) from the lemma, then using ρ + V = j k n k where j i = −(Q U) ik n k , we have, Expanding the covariant derivatives in (D.12) then canceling terms gives, The factor of N in (D.13) vanishes identically; this is again a consequence of the lemma noting that j i L n n i = 0. Therefore, (D.14) To rewrite the above equation in a more symmetric form (exhibiting the g ↔ f duality) we note that for any X µ ν and ε µ holds the identity, A similar expression can be stated on the spatial slice, yielding, 16) Substituting (D.16) into (D.14) gives, All the covariant derivatives in (D.17) can be easily converted to partial derivatives of the respective metric since γU and ϕ U are symmetric.   The HR/HL variables are derived by varying the action in 3+1 form, while the variables in this work are obtained by projecting the field equations. Based on Table 3, we conclude,