Bimetric gravity can reproduce the accelerated expansion of the Universe, without a cosmological constant. However, the stability of these solutions to linear perturbations has been questioned, suggesting exponential growth of structure in this approximation. We present a simple model of structure formation, for which an analytical solution is derived. The solution is well-behaved, showing that there is no physical instability with respect to these perturbations. The model can yield a growth of structure exhibiting measurable differences from $\Lambda$CDM.

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.

The choice of gauge in numerical relativity is crucial in avoiding coordinate and curvature singularities. In addition, the gauge can affect the well-posedness of the system. In this work, we consider the mean gauges, established with respect to the geometric mean metric $h = g \, (g^{-1}f)^{1/2}$ in bimetric relativity. We consider three gauge conditions widely used in numerical relativity and compute them with respect to the geometric mean: The 1+log gauge condition and the maximal slicing for the lapse function of $h$, and the $\Gamma$-driver gauge condition for the shift vector of $h$. In addition, in the bimetric covariant BSSN formalism, there are other arbitrary choices to be made before evolving the system. We show that it is possible to make them by using the geometric mean metric, which is determined dynamically by the system, rather than using an arbitrary external metric, as in general relativity. These choices represent opportunities to recast the system in a well-posed form.

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

We study general properties of static and spherically symmetric bidiagonal black holes in Hassan-Rosen bimetric theory. In particular, we explore the behaviour of the black hole solutions both at the common Killing horizon and at the large radii. The former study leads to a new classification for black holes within the bidiagonal ansatz. The latter study shows that, among the great variety of the black hole solutions, the only solutions converging to Minkowski, Anti-de Sitter and de Sitter spacetimes at large radii are those of General Relativity, i.e., the Schwarzschild, Schwarzschild-Anti-de Sitter and Schwarzschild-de Sitter solutions.