The ghost-free bimetric theory describes interactions of gravity with another spin-2 field in terms of two Lorentzian metrics. However, if the two metrics do not admit compatible notions of space and time, the formulation of the initial value problem becomes problematic. Furthermore, the interaction potential is given in terms of the square root of a matrix which is in general nonunique and possibly nonreal. In this paper we show that both these issues are evaded by requiring reality and general covariance of the equations. First we prove that the reality of the square root matrix leads to a classification of the allowed metrics in terms of the intersections of their null cones. Then, the requirement of general covariance further restricts the allowed metrics to geometries that admit compatible notions of space and time. It also selects a unique definition of the square root matrix. The restrictions are compatible with the equations of motion. These results ensure that the ghost-free bimetric theory can be defined unambiguously and that the two metrics always admit compatible 3+1 decompositions, at least locally. In particular, these considerations rule out certain solutions of massive gravity with locally Closed Causal Curves, which have been used to argue that the theory is acausal.

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.

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.