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.