Gravitational collapse and odd-parity black hole perturbations in Minimal Theory of Bigravity

Masato Minamitsuji Centro de Astrofísica e Gravitação - CENTRA, Departamento de Física, Instituto Superior Técnico - IST, Universidade de Lisboa - UL, Avenida Rovisco Pais 1, 1049-001 Lisboa, Portugal Antonio De Felice Center for Gravitational Physics and Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto University, 606-8502, Kyoto, Japan Shinji Mukohyama Center for Gravitational Physics and Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto University, 606-8502, Kyoto, Japan Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study, The University of Tokyo, Kashiwa, Chiba 277-8583, Japan Michele Oliosi Center for Gravitational Physics and Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto University, 606-8502, Kyoto, Japan

Abstract

We investigate dynamical properties of static and spherically symmetric systems in the self-accelerating branch of the Minimal Theory of Bigravity (MTBG). In the former part, we study the gravitational collapse of pressureless dust and find special solutions, where, in both the physical and fiducial sectors, the exterior and interior spacetime geometries are given by the Schwarzschild spacetimes and the Friedmann-Lemaître-Robertson-Walker universes dominated by pressureless dust, respectively, with specific time slicings. In the case where the Lagrange multipliers are trivial and have no jump across the matter interfaces in both the physical and fiducial sectors, the junction conditions across them remain the same as those in general relativity (GR). For simplicity, we foliate the interior geometry by homogeneous and isotropic spacetimes. For a spatially flat interior universe, we foliate the exterior geometry by a time-independent flat space, while for a spatiallycurved interior universe, we foliate the exterior geometry by a time-independent space with deficit solid angle. Despite the rather restrictive choice of foliations, we find interesting classes of exact solutions that represent gravitational collapse in MTBG. In the spatially flat case, under a certain tuning of the initial condition, we find exact solutions of matter collapse in which the two sectors evolve independently. In the spatially closed case, once the matter energy densities and the Schwarzschild radii are tuned between the two sectors, we find exact solutions that correspond to the Oppenheimer-Snyder model in GR. In the latter part, we study odd-parity perturbations of the Schwarzschild $-$ de Sitter solutions written in the spatially flat coordinates. For the higher-multipole modes $\ell\geq 2$ , we find that, in general, the system reduces to that of four physical modes, where two of them are dynamical and the remaining two are shadowy, i.e., satisfying only elliptic equations. In the case where the ratio of the lapse functions between the physical and fiducial sectors is equal to a constant determined by the parameters of the theory, the two dynamical modes are decoupled from each other but sourced by one of the shadowy modes. Otherwise, the two dynamical modes are coupled to each other and sourced by the two shadowy modes. At least for the cases of collapse described in this paper, we find that the ratio of the lapse functions is determined by the properties of the collapse itself. On giving appropriate boundary conditions to the shadowy modes so as to not strongly backreact/influence the dynamics of the master variables, in the high frequency and short wavelength limits, we show that the two dynamical modes do not suffer from ghost or gradient instabilities. For the dipolar mode $\ell=1$ , the two copies of the slow-rotation limit of the Kerr $-$ de Sitter metrics cannot be a solution in the self-accelerating branch, unless the mass and spin of black holes and effective cosmological constants are tuned to be the same. Therefore, deviation from GR is expected for rotating black holes in the self-accelerating branch of MTBG.

^†^†preprint: YITP-22-157, IPMU22-0068

I Introduction

Modified gravity theories have been proposed from theoretical and observational points of view [1, 2]. While general relativity (GR) has passed all the experimental tests so far in the weak-field regime [3], a new window for testing gravitational theories has opened with the dawn of gravitational-wave astronomy [4, 5]. Massive and bigravity theories are promising candidates to elucidate the origin of the present day’s cosmic acceleration. The first model of massive gravity that is free from the ghost instability was formulated in the context of the linearized gravity by Fierz and Pauli [6]. The attempts to extend the Fierz-Pauli theory to the nonlinear case have not been successful for a long time, because of the reintroduction of a ghost degree of freedom (DOF) through nonlinear higher-derivative interactions in the scalar sector, which is known as the Boulware–Deser (BD) ghost [7]. The first construction of nonlinear massive gravity theory free from the BD ghost was eventually provided by de Rham, Gabadadze, and Tolley (dRGT) [8]. dRGT massive gravity was then extended to a bigravity theory by Hassan and Rosen (HR) [9] by promoting the fiducial metric to be a dynamical field. ¹¹1HR bigravity, however, still suffers from the BD ghost, when matter is coupled to both the physical and fiducial metrics [10, 11]. See also [12].

The Minimal Theory of Massive Gravity (MTMG) [13, 14] is a different extension of dRGT massive gravity. In MTMG, there are two tensorial DOFs as in GR. While the four-dimensional diffeomorphism invariance is explicitly broken, as the fiducial metric is given as a function of the coordinates, the absence of otherwise problematic scalar and vector DOFs makes it easy for the theory to be consistent with experimental and observational data. For instance, while MTMG shares the same homogeneous and isotropic background cosmology with dRGT massive gravity, including the existence of two separate branches, it does not suffer from the instabilities which appear in dRGT massive gravity [15, 16, 17], thanks indeed to the absence of the extra scalar and vector DOFs. The cosmological solutions in MTMG are, as in dRGT, classified into the normal and self-accelerating branches [18, 19, 20]. Any solution in GR with or without matter written in the spatially flat coordinates has been shown to be a solution in the self-accelerating branch of MTMG [21].

Using a construction similar to MTMG, HR bigravity was then extended to the Minimal Theory of Bigravity (MTBG) [22], where the four-dimensional diffeomorphism invariance is broken down to the three-dimensional one and time-reparametrization invariance. While, by construction, MTBG shares the same background cosmological dynamics with HR bigravity, the number of propagating DOFs are down to four, two tensorial DOFs in the physical sector and the other two in the fiducial sector. The absence of the extra scalar and vector DOFs in MTBG implies the absence of ghost or gradient instabilities associated with them [23, 24, 25]. In the case where the matter sectors are coupled to the two metrics separately and independently, cosmology in MTBG has been investigated in Ref. [22], which showed that in the self-accelerating branch both background cosmology and dynamics of the scalar and vector cosmological linear perturbations behave in the same way as in two copies of GR. On the other hand, two of the tensor modes acquire a nonzero mass. In the normal branch, another nontrivial difference arises in addition to the presence of two massive tensor modes: deviation in the dynamics for both background and the scalar sector from GR could be already observed and possibly lead to nontrivial interesting phenomenology. The absence of extra scalar and vector modes in MTBG also helps in developing a new production scenario of spin-2 dark matter based on the transition from an anisotropic fixed point solution to an isotropic one within axisymmetric Bianchi type-I universe [26].

MTMG and MTBG are theories with constraints. By construction, they possess constraints out of which the unwanted modes can be removed nonlinearly from the theory. This nice feature is motivated by the fact that whatever new theory of gravity we introduce should reproduce the whole range of gravitational phenomenology we already know. However, the introduction of such constraints is a nontrivial task. On one hand, the choice of constraints is not unique. Other MTMGs and/or MTBGs can be, in principle, introduced, and a priori it is not an easy task to discriminate one choice of constraints from another. On the other hand, each set of constraints we pick up for MTMG, including the one we will study in this paper, may restrict the space of solutions of the theory in a different way. While this is to be expected, by the same nature of a constraint, once more, it is not clear, a priori, how to determine which configurations are not allowed to be solutions of the specific choice of theory.

Despite the fact that this fact by itself would not be already introducing new and/or unexpected phenomenology, there is also another feature shared by these theories: the presence of shadowy modes [27]. In fact, as we will also see in this paper, the constraints lead to the presence of modes that only exist on a spacelike three-dimensional hypersurface. These modes are present not only in MTMG and MTBG, but also in other theories: for instance Refs. [28, 29], which also allow for the existence of shadowy modes. They all share the property that they satisfy elliptic equations for which some appropriate boundary conditions need to be imposed. As a consequence, the very existence of such modes introduces a preferred frame for the theory, i.e., the frame on which their equations of motion are manifestly elliptic. This leads to the consequence that, even for a given GR solution, different choices of frames may or may not be compatible with the existence of the shadowy modes; i.e., some slicings of GR solutions may be allowed as solutions by the theory, but others may not. This fact leads, for instance, to the absence of the Birkhoff theorem for these theories. Therefore, it is of primary importance to try to find all the solutions for a given symmetry or physical configuration.

Some of these foliation-and-branch specific solutions of MTBG have already been found. For instance, in a previous work [30], we have investigated the static and spherically symmetric solutions in the self-accelerating and normal branches of MTBG. We have shown that a pair of Schwarzschild $-$ de Sitter spacetimes with different cosmological constants and black hole masses written in the spatially flat [Gullstrand–Painlevé (GP)] coordinates is a solution in the self-accelerating branch of MTBG. In the case where the two matter sectors are coupled to the two metrics separately and independently, the self-accelerating branch also admits static and spherically symmetric copies of two GR solutions written in the spatially flat coordinates, which include neutron stars with arbitrary equations of state of matter. In addition to these solutions we also find nontrivial solutions that correspond to Schwarzschild $-$ de Sitter solutions, written in nonstandard coordinates. It should be added that these solutions have been found by imposing trivial configurations/profiles for the Lagrangian multipliers, introduced in the theory so as to impose the same above-mentioned constraints that define the theory itself. The same approach will be taken in this paper. It is not clear what kinds of new solutions are to be expected once we move away from this simple ansatz for the Lagrange multipliers.

On the other hand, in the normal branch, while the spatially flat coordinates of the Schwarzschild $-$ de Sitter metrics cannot be solutions, those written in the coordinates with $D^{i}D_{i}K=0$ , where $D_{i}$ is the covariant derivative associated with the spatial metric and $K$ is the trace of the extrinsic curvature tensor on the constant-time hypersurfaces, could be solutions, provided that the two metrics are parallel, showing also for this case a severe fine-tuning of the possible configurations. An example for such solutions were given, i.e., the ones with $K={\rm constant}$ , which turn out to be also solutions for the case of a single physical metric, in the context of the VCDM models that replace the cosmological constant $\Lambda$ in the $\Lambda$ CDM ( $\Lambda$ cold dark matter) model with a function of an auxiliary field $\phi$ , $V(\phi)$ [28, 31, 32, 33].

As the next step, in this paper we will study the dynamical processes only in the self-accelerating branch of MTBG and investigate whether any deviation from the GR predictions can be observed with the inclusion of time dependence at the level of the background solutions and perturbations, respectively.

Our analysis will be composed of two parts. In the former, we will study spherical gravitational collapse of pressureless dust in the self-accelerating branch of MTBG. In the case where the two species of matter are independently coupled to the two metrics separately and independently, we will derive special solutions describing spherical collapse of pressureless dust in both the physical and fiducial sectors, where the interior regions of the matter distribution are described by the spatially flat or closed Friedmann-Lemaître-Robertson-Walker (FLRW) universes and the exterior geometry is the Schwarzschild spacetime with specific time slicings. For simplicity, we foliate the interior geometry by homogeneous and isotropic spacetimes. For a spatially flat interior universe, we foliate the exterior geometry by a time-independent flat space, while for a spatiallycurved interior universe, we foliate the exterior geometry by a time-independent space with deficit solid angle. Despite the rather restrictive choice of foliations, we find interesting classes of exact solutions that represent gravitational collapse in MTBG. The spatially flat solutions are obtained under certain tuning of the initial conditions of the collapse. The spatially closed case corresponds to an extension of the Oppenheimer-Snyder model in GR [34, 35, 36] to MTBG. We will show that, in this case, gravitational collapse happens in the physical and fiducial sectors in the same manner as in two copies of GR independently, under certain tuning of the matter energy densities and Schwarzschild radii between the two sectors. Needless to say, these tunings reflect the fact that we restrict our considerations to rather specific choice of time slicings.

In the latter part, we will investigate the odd-parity perturbations of the Schwarzschild $-$ de Sitter solutions written in the spatially flat coordinates in the self-accelerating branch of MTBG. For the higher-multipolar modes $\ell\geq 2$ , where $\ell$ represents the angular multipole moment, we derive the master equations governing the dynamics of the odd-parity perturbations in both physical and fiducial sectors in MTBG. We will find that there are four physical modes, where two of them are dynamical and the remaining two are shadowy, which obey elliptic equations and are fixed by the spatial boundary conditions on each step of the time evolution [27, 29]. In the case where the ratio of the lapse functions in the physical and fiducial sectors are equal to a constant determined by the free parameters of theory, the two dynamical modes are decoupled and sourced by one of the shadowy modes. Otherwise, the two dynamical modes are coupled to each other and sourced by the two shadowy modes. In the high frequency and short wavelength limits, we also verify the modes are not suffering from ghost or gradient instabilities, provided we can neglect backreactions from the shadowy modes, by giving appropriate boundary conditions to them. For the dipolar mode $\ell=1$ , we will show that the two copies of the slow-rotation limit of the Kerr $-$ de Sitter metrics, in general, cannot be a solution in MTBG, unless the mass and spin of the black holes and the effective cosmological constants in the two sectors are tuned to be the same. Therefore, deviation from GR is expected for rotating black holes in the self-accelerating branch of MTBG. Readers should also refer to the last paragraph in Sec.IV of [37] for a similar statement in the context of MTMG.

The structure of this paper is as follows: In Sec. II, we review the MTBG theory. In Sec. III, we derive the exact solutions that describe gravitational collapse of pressureless dust where the metrics in the interior regions filled with matter are written in terms of the spatially flat and spatially closed FLRW universes, respectively. In Sec. IV, we investigate the odd-parity gravitational perturbations of the Schwarzschild $-$ de Sitter solutions written in the spatially flat coordinates in the self-accelerating branch of MTBG. The last section is devoted to giving a brief summary and conclusion.

II The Minimal Theory of Bigravity

In this section, we review MTBG following Refs. [22, 30].

II.1 Metrics

As in HR bigravity, MTBG is composed of two metric sectors, the physical and fiducial metrics denoted by $g_{\mu\nu}$ and $f_{\mu\nu}$ , respectively. Choosing the unitary gauge, the two metrics $g_{\mu\nu}$ and $f_{\mu\nu}$ are expressed in the Arnowitt-Deser-Misner form, respectively, as

	$\displaystyle g_{\mu\nu}dx^{\mu}dx^{\nu}=-N^{2}dt^{2}+\gamma_{ij}\left(dx^{i}+N^{i}dt\right)\left(dx^{j}+N^{j}dt\right),$
	$\displaystyle f_{\mu\nu}dx^{\mu}dx^{\nu}=-M^{2}dt^{2}+\phi_{ij}\left(dx^{i}+M^{i}dt\right)\left(dx^{j}+M^{j}dt\right),$		(1)

where $t$ and $x^{i}$ ( $i=1,2,3$ ) are the temporal and spatial coordinates, and $(N,N^{i},\gamma_{ij})$ and $(M,M^{i},\phi_{ij})$ are the sets of the lapse function, shift vector, and spatial metric on the constant $t$ hypersurfaces for the metrics $g_{\mu\nu}$ and $f_{\mu\nu}$ , respectively. The extrinsic curvature tensors on constant-time hypersurfaces are defined by

	$\displaystyle K_{ij}:=\frac{1}{2N}\left(\partial_{t}\gamma_{ij}-D_{i}N_{j}-D_{j}N_{i}\right),$		(2)
	$\displaystyle\Phi_{ij}:=\frac{1}{2M}\left(\partial_{t}\phi_{ij}-{\tilde{D}}_{i}M_{j}-{\tilde{D}}_{j}M_{i}\right),$		(3)

where $D_{i}$ and ${\tilde{D}}_{i}$ are covariant derivatives associated with the spatial metrics $\gamma_{ij}$ and $\phi_{ij}$ , respectively.

II.2 Action and equations of motion

In the unitary gauge, the action of MTBG [22] is then given by

\displaystyle S

\displaystyle=

\displaystyle\frac{1}{2\kappa^{2}}\int d^{4}x\left({\cal L}_{g}\left[N,N^{i},\gamma_{ij};M,M^{i},\phi_{ij};\lambda,{\bar{\lambda}},\lambda^{i}\right]+{\cal L}_{m}\left[N,N^{i},\gamma_{ij};M,M^{i},\phi_{ij};\Psi\right]\right),

(4)

where $\kappa^{2}$ represents the gravitational constant in the physical sector, ${\cal L}_{g}$ and ${\cal L}_{m}$ represent the gravitational and matter parts of the Lagrangian, respectively, $\lambda$ , ${\bar{\lambda}}$ , and $\lambda^{i}$ are two scalar and one spatial-vector Lagrange multipliers, and $\Psi$ represents the matter sector.

The gravitational Lagrangian ${\cal L}_{g}$ of MTBG is further decomposed into the precursor and constraint parts as

\displaystyle{\cal L}_{g}

\displaystyle:=

\displaystyle{\cal L}_{\rm pre}\left[N,N^{i},\gamma_{ij};M,M^{i},\phi_{ij}\right]+{\cal L}_{\rm con}\left[N,N^{i},\gamma_{ij};M,M^{i},\phi_{ij};\lambda,{\bar{\lambda}},\lambda^{i}\right],

(5)

with

$\displaystyle{\cal L}_{\rm pre}$	$\displaystyle:=$	$\displaystyle\sqrt{-g}R[g]+\tilde{\alpha}^{2}\sqrt{-f}R[f]-m^{2}\left(N\sqrt{\gamma}{\cal H}_{0}+M\sqrt{\phi}\tilde{\cal H}_{0}\right),$
$\displaystyle{\cal L}_{\rm con}$	$\displaystyle:=$	$\displaystyle\sqrt{\gamma}\alpha_{1\gamma}\left(\lambda+\Delta_{\gamma}{\bar{\lambda}}\right)+\sqrt{\phi}\alpha_{1\phi}\left(\lambda-\Delta_{\phi}{\bar{\lambda}}\right)+\sqrt{\gamma}\alpha_{2\gamma}\left(\lambda+\Delta_{\gamma}{\bar{\lambda}}\right)^{2}+\sqrt{\phi}\alpha_{2\phi}\left(\lambda-\Delta_{\phi}{\bar{\lambda}}\right)^{2}$	(6)
	$\displaystyle-$	$\displaystyle m^{2}\left[\sqrt{\gamma}U^{i}{}_{k}D_{i}\lambda^{k}-\beta\sqrt{\phi}{\tilde{U}}_{k}{}^{i}{\tilde{D}}_{i}\lambda^{k}\right],$	(6)

and

	$\displaystyle\alpha_{1\gamma}$	$\displaystyle:=$	$\displaystyle-m^{2}U^{p}{}_{q}K^{q}{}_{p},\quad\alpha_{1\phi}:=m^{2}{\tilde{U}}^{p}{}_{q}\Phi^{q}{}_{p},$
	$\displaystyle\alpha_{2\gamma}$	$\displaystyle:=$	$\displaystyle\frac{m^{4}}{4N}\left(U^{p}{}_{q}-\frac{1}{2}U^{k}{}_{k}\delta^{p}{}_{q}\right)U^{q}{}_{p},\quad\alpha_{2\phi}:=\frac{m^{4}}{4M\tilde{\alpha}^{2}}\left({\tilde{U}}_{q}{}^{p}-\frac{1}{2}{\tilde{U}}_{k}{}^{k}\delta_{q}{}^{p}\right){\tilde{U}}_{p}{}^{q},$		(7)

where the constant $\tilde{\alpha}$ represents the ratio of the two gravitational constants, $m$ is a parameter with dimensions of mass that can be related with the graviton mass, $\beta$ is a constant, and $\gamma:={\rm det}(\gamma_{ij})$ and $\phi:={\rm det}(\phi_{ij})$ are the determinants of the three-dimensional spatial metrics $\gamma_{ij}$ and $\phi_{ij}$ respectively. We also note that $K^{q}{}_{p}=\gamma^{qr}K_{rp}$ and $\Phi^{q}{}_{p}=\phi^{qr}\Phi_{rp}$ . Furthermore, ${\cal H}_{0}$ and $\tilde{\cal H}_{0}$ are defined by ${\cal H}_{0}:=\sum_{n=0}^{3}c_{4-n}e_{n}({\cal K})$ and $\tilde{\cal H}_{0}:=\sum_{n=0}^{3}c_{n}e_{n}(\tilde{\cal K})$ with

\displaystyle e_{0}({\cal K})=1,\quad e_{1}({\cal K})=\left[{\cal K}\right],\quad e_{2}({\cal K})=\frac{1}{2}\left(\left[{\cal K}\right]^{2}-\left[{\cal K}^{2}\right]\right),\quad e_{3}({\cal K})={\rm det}({\cal K}),

(8)

and similar for $e_{n}(\tilde{\cal K})$ with ${\cal K}^{i}{}_{k}$ and $\tilde{\cal K}_{k}{}^{i}$ characterized by ${\cal K}^{i}{}_{k}{\cal K}^{k}{}_{j}={\gamma}^{ik}\phi_{kj}$ and ${\cal{\tilde{K}}}_{j}{}^{k}{\cal{\tilde{K}}}_{k}{}^{i}=\gamma_{jk}{\phi}^{ki}$ ; $\Delta_{\gamma}:=\gamma^{ij}D_{i}D_{j}$ and $\Delta_{\phi}:=\phi^{ij}{\tilde{D}}_{i}{\tilde{D}}_{j}$ are the Laplacian operators in the physical and fiducial sectors, respectively, and the spatial tensors $U^{i}{}_{j}$ and ${\tilde{U}}{}_{j}{}^{i}$ are defined by

	$\displaystyle U^{i}{}_{j}$	$\displaystyle:=$	$\displaystyle\frac{1}{2}\sum_{n=1}^{3}c_{4-n}\left(U_{(n)}{}^{i}{}_{j}+\gamma^{ik}\gamma_{j\ell}U_{(n)}{}^{\ell}{}_{k}\right),$
	$\displaystyle{\tilde{U}}_{j}{}^{i}$	$\displaystyle:=$	$\displaystyle\frac{1}{2}\sum_{n=1}^{3}c_{n}\left({\tilde{U}}_{(n)j}{}^{i}+\phi^{ik}\phi_{j\ell}{\tilde{U}}_{(n)k}{}^{\ell}\right),$		(9)

with $U_{(n)}{}^{i}{}_{k}:=\frac{\partial e_{n}({\cal K})}{\partial{\cal K}^{k}{}_{i}}$ and ${\tilde{U}}_{(n)k}{}^{i}:=\frac{\partial e_{n}({\cal{\tilde{K}}})}{\partial\tilde{\cal K}^{k}{}_{i}}$ , and $c_{j}$ ( $j=0,1,2,3,4$ ) being dimensionless coupling constants. Although one of $c_{j}$ ( $j=0,1,2,3,4$ ) may be set to unity by a redefinition of the mass parameter $m$ , in order to discuss the various limiting cases of each coefficient, we keep them as independent parameters.

Variation of the action (4) with respect to $N$ , $N^{i}$ , $\gamma_{ij}$ , $M$ , $M^{i}$ , and $\phi_{ij}$ provides the equations of motion for the two metrics, and variation with respect to $\Psi$ provides the equations of motion for the matter field. On the other hand, variation with respect to the Lagrange multipliers $\lambda$ , ${\bar{\lambda}}$ , and $\lambda^{i}$ gives the constraint equations

	$\displaystyle\sqrt{\gamma}\alpha_{1\gamma}+\sqrt{\phi}\alpha_{1\phi}+2\sqrt{\gamma}\alpha_{2\gamma}\left(\lambda+\Delta_{\gamma}{\bar{\lambda}}\right)+2\sqrt{\phi}\alpha_{2\phi}\left(\lambda-\Delta_{\phi}{\bar{\lambda}}\right)=0,$		(10)
	$\displaystyle\sqrt{\gamma}\Delta_{\gamma}\alpha_{1\gamma}-\sqrt{\phi}\Delta_{\phi}\alpha_{1\phi}+2\sqrt{\gamma}\Delta_{\gamma}\left[\alpha_{2\gamma}\left(\lambda+\Delta_{\gamma}{\bar{\lambda}}\right)\right]-2\sqrt{\phi}\Delta_{\phi}\left[\alpha_{2\phi}\left(\lambda-\Delta_{\phi}{\bar{\lambda}}\right)\right]=0,$		(11)
	$\displaystyle\sqrt{\gamma}D_{p}U^{p}{}_{q}-\beta\sqrt{\phi}{\tilde{D}}_{p}{\tilde{U}}_{q}{}^{p}=0.$		(12)

II.3 Matter

We also assume that the matter sector $\Psi$ is decomposed into those in the physical and fiducial sectors $\Psi_{g}$ and $\Psi_{f}$ and the matter Lagrangian ${\cal L}_{m}$ is given by

\displaystyle{\cal L}_{m}\left[N,N^{i},\gamma_{ij};M,M^{i},\phi_{ij};\Psi\right]={\cal L}_{m,g}\left[N,N^{i},\gamma_{ij};\Psi_{g}\right]+{\cal L}_{m,f}\left[M,M^{i},\phi_{ij};\Psi_{f}\right].

(13)

The two matter sectors ${\cal L}_{m,g}$ and ${\cal L}_{m,f}$ are individually coupled to the metrics $g_{\mu\nu}$ and $f_{\mu\nu}$ , respectively, and they are not directly coupled to each other. In other words, the two matter sectors can interact only through gravitation.

III Gravitational collapse of pressureless dust

In this section, we consider two models of collapsing matter in the self-accelerating branch of MTBG. We assume the existence of matter surfaces that move inward in both the sectors separately, and the spacetime metrics inside the matter surfaces are described by either spatially flat or closed FLRW universes composed of pressureless dust. The spatially closed case corresponds to an extension of the Oppenheimer-Snyder model in GR [34, 35, 36] to MTBG.

III.1 Collapse in the spherically symmetric spacetimes in the spatially flat coordinates

First, we consider the case in which the spherically symmetric spacetimes can be expressed in the spatially flat coordinates

	$\displaystyle g_{\mu\nu}dx^{\mu}dx^{\nu}$	$\displaystyle=$	$\displaystyle-dt^{2}+\left(dr+N^{r}(t,r)dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b},$		(14)
	$\displaystyle f_{\mu\nu}dx^{\mu}dx^{\nu}$	$\displaystyle=$	$\displaystyle C_{0}^{2}\left[-b^{2}dt^{2}+\left(dr+N^{r}_{(f)}(t,r)dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b}\right],$		(15)

where $b>0$ and $C_{0}>0$ are constants, $t$ , $r$ , and $\theta^{a}=(\theta,\varphi)$ are the temporal, radial, and angular coordinates, respectively, and $\theta_{ab}$ represents the metric of the unit two-sphere. In vacuum GR with the vanishing cosmological constant, since $N^{r}_{(f)}=\sqrt{r_{(f)}/r}$ , by a redefinition of ${\tilde{r}}=\frac{r}{b}$ , $C_{0}\to C_{0}/b$ , and $r_{(f)}\to r_{(f)}b^{3}$ the metric (15) can be brought to that with $b=1$ . In MTBG, the constant $bC_{0}$ measures the ratio of the lapse functions between the fiducial and physical sectors.

In the spherically symmetric systems, the general ansatz for the Lagrange multipliers is given by

\displaystyle\lambda=\lambda(t,r),\qquad\bar{\lambda}=\bar{\lambda}(t,r),\qquad\lambda^{r}=\lambda^{r}(t,r),\qquad\lambda^{a}=0.

(16)

We assume that in each of the physical and fiducial sectors the spacetime is divided into the interior and exterior regions, where the interior regions are filled with pressureless dust. The interfaces between the two regions follow the trajectories on the $(t,r)$ plane

\displaystyle(t,r)=\left(T_{(g)}(\tau_{(g)}),R_{(g)}(\tau_{(g)})\right),\qquad\left(T_{(f)}(\tau_{(f)}),R_{(f)}(\tau_{(f)})\right),

(17)

where $\tau_{(g)}$ and $\tau_{(f)}$ represent the proper times on the interfaces in the physical and fiducial sectors, and $T_{(g,f)}(\tau_{(g,f)})<0$ and $R_{(g,f)}(\tau_{(g,f)}>0$ are smooth functions of $\tau_{(g,f)}$ representing the temporal and radial positions of the interfaces in both the sectors, respectively. The four velocities of the interfaces are, respectively, given by

\displaystyle u_{(g)}^{\mu}=\left(\dot{T}_{(g)},\dot{R}_{(g)},0,0\right),\qquad u_{(f)}^{\mu}=\left(\dot{T}_{(f)},\dot{R}_{(f)},0,0\right),

(18)

where the dot represents the derivative with respect to the proper time $\tau_{(g)}$ or $\tau_{(f)}$ in each sector, respectively. Since matter is collapsing, the interfaces are moving inward in the radial directions, respectively, $\dot{R}_{(g)}<0$ and $\dot{R}_{(f)}<0$ . The trajectories of the interfaces are timelike $g_{\mu\nu}u_{(g)}^{\mu}u_{(g)}^{\nu}=f_{\mu\nu}u_{(f)}^{\mu}u_{(f)}^{\nu}=-1$ , which can be solved as

	$\displaystyle\dot{T}_{(g)}$	$\displaystyle=$	$\displaystyle\frac{1}{1-\left(N^{r}{}\right)^{2}}\left[N^{r}\dot{R}_{(g)}+\sqrt{1-\left(N^{r}{}\right)^{2}+\dot{R}_{(g)}^{2}}\right],$		(19)
	$\displaystyle\dot{T}_{(f)}$	$\displaystyle=$	$\displaystyle\frac{1}{C_{0}\left[b^{2}-(N_{(f)}^{r}{})^{2}\right]}\left[C_{0}N_{(f)}^{r}\dot{R}_{(f)}+\sqrt{b^{2}-(N_{(f)}^{r}{})^{2}+C_{0}^{2}b^{2}\dot{R}_{(f)}^{2}}\right].$		(20)

For the unique definition of the time derivatives (19) and (20), $N^{r}$ and $N^{r}_{(f)}$ defined in the interior and exterior regions have to coincide at the interfaces.

The unit normals to the interfaces satisfying the normalization and orthogonality conditions $g_{\mu\nu}n_{(g)}^{\mu}n_{(g)}^{\nu}=f_{\mu\nu}n^{\mu}_{(f)}n_{(f)}^{\nu}=1$ and $g_{\mu\nu}u_{(g)}^{\mu}n_{(g)}^{\nu}=f_{\mu\nu}u_{(f)}^{\mu}n_{(f)}^{\nu}=0$ point outward in the radial direction. The induced metrics on the interfaces are then given by

	$\displaystyle h_{ij}dy^{i}dy^{j}=g_{\mu\nu}e^{(g)\mu}{}_{i}e^{(g)\nu}{}_{j}dy^{i}dy^{j}=-d\tau^{2}+R_{(g)}^{2}\theta_{ab}d\theta^{a}d\theta^{b},$		(21)
	$\displaystyle k_{ij}dy^{i}dy^{j}=f_{\mu\nu}e^{(f)\mu}{}_{i}e^{(f)\nu}{}_{j}dy^{i}dy^{j}=-d\tau_{(f)}^{2}+C_{0}^{2}R_{(f)}^{2}\theta_{ab}d\theta^{a}d\theta^{b},$		(22)

where the nonzero components of the projection tensors are given by $e^{(g)t}{}_{\tau}=\dot{T}_{(g)}$ , $e^{(g)r}{}_{\tau}=\dot{R}_{(g)}$ , $e^{(g)a}{}_{b}=\delta^{a}{}_{b}$ , $e^{(f)t}{}_{\tau}=\dot{T}_{(f)}$ , $e^{(f)r}{}_{\tau}=\dot{R}_{(f)}$ , and $e^{(f)a}{}_{b}=\delta^{a}{}_{b}$ with the indices $a,b$ representing the angular directions.

The extrinsic curvature tensors on the interfaces in the physical and fiducial sectors are given by

\displaystyle\tilde{K}^{(g)}_{ij}=e^{(g)\mu}{}_{i}e^{(g)\nu}{}_{j}\nabla^{(g)}_{\mu}n^{(g)}_{\nu},\quad\tilde{K}^{(f)}_{ij}=e^{(f)\mu}{}_{i}e^{(f)\nu}{}_{j}\nabla^{(f)}_{\mu}n^{(f)}_{\nu},

(23)

whose nonzero components are given by

$\displaystyle\tilde{K}^{(g)}_{\tau\tau}$	$\displaystyle=$	$\displaystyle\frac{1}{\sqrt{1-\left(N^{r}{}\right)^{2}+\dot{R}_{(g)}^{2}}}\left\{-\ddot{R}_{(g)}+\frac{\left(1-(N^{r})^{2}\right)^{2}N^{r}N^{r}_{,r}-\left(N^{r}\dot{R}_{(g)}+\sqrt{1+\dot{R}_{(g)}^{2}-(N^{r})^{2}}\right)^{2}N^{r}_{,t}}{\left(1-(N^{r})^{2}\right)^{2}}\right\},$
$\displaystyle\tilde{K}^{(g)}_{ab}$	$\displaystyle=$	$\displaystyle R_{(g)}\sqrt{1-\left(N^{r}{}\right)^{2}+\dot{R}_{(g)}^{2}}\theta_{ab},$
$\displaystyle\tilde{K}^{(f)}_{\tau\tau}$	$\displaystyle=$	$\displaystyle\frac{1}{\sqrt{b^{2}-(N_{(f)}^{r}{})^{2}+b^{2}C_{0}^{2}\dot{R}_{(f)}^{2}}}$
	$\displaystyle\times$	$\displaystyle\left\{-bC_{0}\ddot{R}_{(f)}+\frac{\left(b^{2}-(N_{(f)}^{r})^{2}\right)^{2}N_{(f)}^{r}N^{r}_{(f),r}-b^{2}\left(C_{0}N_{(f)}^{r}\dot{R}_{(f)}+\sqrt{b^{2}+b^{2}C_{0}^{2}\dot{R}_{(f)}^{2}-(N_{(f)}^{r})^{2}}\right)^{2}N^{r}_{(f),t}}{bC_{0}\left(b^{2}-(N_{(f)}^{r})^{2}\right)^{2}}\right\},$
$\displaystyle\tilde{K}^{(f)}_{ab}$	$\displaystyle=$	$\displaystyle\frac{C_{0}R_{(f)}}{b}\sqrt{b^{2}-(N_{(f)}^{r}{})^{2}+b^{2}C_{0}^{2}\dot{R}_{(f)}^{2}}\theta_{ab}.$	(24)

We consider the background solutions where Lagrange multipliers are trivial in both the interior and exterior regions and continuous across the interfaces

\displaystyle\lambda=\lambda^{r}=0,\qquad{\bar{\lambda}}={\bar{\lambda}}_{0}={\rm constant}.

(25)

Since the contributions of the constraints to the metric equations of motion do not contain second-order derivatives of the metric functions, in the case where the Lagrange multipliers are trivial and continuous across the interfaces, the junction conditions in both the physical and fiducial sectors are identical to those in the two copies of GR.

The first junction conditions correspond to the continuity of the induced metrics on the interfaces

\displaystyle\left[h_{ij}\right]:=h^{(+)}_{ij}-h^{(-)}_{ij}=0,\qquad\left[k_{ij}\right]:=k^{(+)}_{ij}-k^{(-)}_{ij}=0,

(26)

where below $(+)$ and $(-)$ denote the regions outside and inside the interfaces, which means the continuity of $N^{r}$ and $N^{r}_{(f)}$ across the interfaces

\displaystyle N^{r(+)}\left(T_{(g)},R_{(g)}\right)=N^{r(-)}\left(T_{(g)},R_{(g)}\right),\qquad N_{(f)}^{r(+)}\left(T_{(f)},R_{(f)}\right)=N_{(f)}^{r(-)}\left(T_{(f)},R_{(f)}\right).

(27)

The second junction conditions are the jump of the extrinsic curvature tensors across the interfaces

	$\displaystyle\left[{\tilde{K}}^{(g)}_{ij}\right]:={\tilde{K}}^{(g)(+)}_{ij}-{\tilde{K}}^{(g)(-)}_{ij}=\kappa^{2}\left(S^{(g)}_{ij}-\frac{1}{2}S^{(g)}h_{ij}\right),$
	$\displaystyle\left[{\tilde{K}}^{(f)}_{ij}\right]:={\tilde{K}}^{(f)(+)}_{ij}-{\tilde{K}}^{(f)(-)}_{ij}=\frac{\kappa^{2}}{{\tilde{\alpha}}^{2}}\left(S^{(f)}_{ij}-\frac{1}{2}S^{(f)}k_{ij}\right),$		(28)

where $S^{(g)}_{ij}$ and $S^{(f)}_{ij}$ represent the surface energy-momentum tensors of matter localized on the interfaces.

III.1.1 The exterior and interior solutions

For the spatially flat coordinates (14) and (15), the solutions of the constraint equations (10)-(12) can be divided into the self-accelerating and normal branches [30]. We focus on the self-accelerating branch given by the condition

\displaystyle C_{0}^{2}c_{1}+2C_{0}c_{2}+c_{3}=0.

(29)

Note that the Schwarzschild solutions written in the Schwarzschild coordinates cannot describe black hole solutions in MTBG, as they are singular at the positions of the event horizon [30].

Outside the interfaces, $r>R_{(g)}$ and $r>R_{(f)}$ , the metric solution is described by the Schwarzschild solutions [30] written in the spatially flat coordinates (14) and (15),

\displaystyle N^{r(+)}(r)=\sqrt{\frac{r_{(g)}}{r}},\qquad N_{(f)}^{r(+)}(r)=\sqrt{\frac{r_{(f)}}{r}},

(30)

where $r_{(g)}$ and $r_{(f)}$ represent the Schwarzschild radii in both the physical and fiducial sectors, respectively, which exist under the conditions for asymptotically Minkowski spacetimes

	$\displaystyle C_{0}^{3}c_{1}+3C_{0}^{2}c_{2}+3C_{0}c_{3}+c_{4}=0,$		(31)
	$\displaystyle C_{0}^{3}c_{0}+3C_{0}^{2}c_{1}+3C_{0}c_{2}+c_{3}=0,$		(32)

together with Eq. (29). Note that the Schwarzschild solutions exist for an arbitrary $b>0$ .

Inside the interfaces ( $r<R_{(g)}$ and $r<R_{(f)}$ ), we assume that the metric solution is described by the spatially flat coordinates (14) and (15) with

\displaystyle N^{r(-)}(t,r)=-r\frac{a_{,t}}{a(t)},\qquad N_{(f)}^{r(-)}(t,r)=-r\frac{a_{,t}}{a(t)},

(33)

where $a_{,t}\equiv da/dt$ . Introducing the comoving radial coordinates $t=t_{c}$ , $r=r_{c}a(t)$ , the interior solutions can be expressed by the spatially flat FLRW metrics

	$\displaystyle g_{\mu\nu}dx^{\mu}dx^{\nu}=-dt_{c}^{2}+a(t_{c})^{2}\left(dr_{c}^{2}+r_{c}^{2}\theta_{ab}d\theta^{a}d\theta^{b}\right),$		(34)
	$\displaystyle f_{\mu\nu}dx^{\mu}dx^{\nu}=C_{0}^{2}\left[-b^{2}dt_{c}^{2}+a(t_{c})^{2}\left(dr_{c}^{2}+r_{c}^{2}\theta_{ab}d\theta^{a}d\theta^{b}\right)\right].$		(35)

The fact that for the two metrics the ratio of the two effective scale factors is a constant $C_{0}$ , satisfying (29), is a consequence of the constraints in MTBG for the self-accelerating branch. We assume that the spacetimes are filled with pressureless dust fluids whose energy-momentum tensors are given by

\displaystyle T_{(g)}^{\mu\nu}=\rho_{(g)}v_{(g)}^{\mu}v_{(g)}^{\nu},\qquad T_{(f)}^{\mu\nu}=\rho_{(f)}v_{(f)}^{\mu}v_{(f)}^{\nu},

(36)

where $\rho_{(g)}$ and $\rho_{(f)}$ represent the matter energy densities that behave as

\displaystyle\rho_{(g)}=\frac{\rho_{(g,0)}}{a(t)^{3}},\qquad\rho_{(f)}=\frac{\rho_{(f,0)}}{a(t)^{3}},

(37)

with $\rho_{(g,0)}>0$ and $\rho_{(f,0)}>0$ being constants, and the four velocities of matter $v_{(g)}^{\mu}$ and $v_{(f)}^{\mu}$ are comoving with the Hubble flows of the spacetimes²²2In comoving coordinates, along the radial geodesics of matter, $r_{c}$ is constant. In this case, the normalized four-velocity is given by $v^{\mu}_{(g)}\partial_{\mu}=\partial_{t_{c}}$ . For the change of variables between $(t_{c},r_{c})$ to $(t,r)$ one can show that $\partial_{t_{c}}=\partial_{t}+ra_{,t}/a\partial_{r}$ , from which we obtain the first of Eq. (38). Along the same lines, one finds the second one of Eq. (38).

\displaystyle v_{(g)}^{\mu}=\left(1,r\frac{a_{,t}}{a},0,0\right),\qquad v_{(f)}^{\mu}=\frac{1}{bC_{0}}\left(1,r\frac{a_{,t}}{a}0,0\right).

(38)

The only nontrivial components of the gravitational equations in both the sectors of MTBG set the following constraints:

\displaystyle\frac{3a_{,t}^{2}}{a^{2}}=\kappa^{2}\frac{\rho_{(g,0)}}{a^{3}},\qquad b^{2}=\frac{\tilde{\alpha}^{2}}{C_{0}^{2}}\,\frac{\rho_{(g,0)}}{\rho_{(f,0)}}.

(39)

This implies that the second Friedmann equation effectively sets the value of $b$ in terms of the matter content on both the sectors. Furthermore, this result states that the value of $b$ for the exterior metric, which was seen as a free parameter, is actually determined by the details of collapsing matter. If matter is present only in one of the two sectors, the value of the parameter $b$ is still fixed through the relation (39), although its value becomes either zero or divergent. We note that in such a case zero or divergent $b$ is obtained only if the effective cosmological constants $\Lambda_{(g)}$ and $\Lambda_{(f)}$ [see Eq. (99)] are tuned to zero. On the other hand, if $\Lambda_{(g)}$ and $\Lambda_{(f)}$ take positive values then the value of $b$ is fixed to be a nonzero finite value even in the case where one or both of the matter energy densities vanishes [12]. On the other hand, if both $\Lambda_{(g)}$ and $\Lambda_{(f)}$ are exactly zero, then the value of $b$ in the limit of vanishing energy densities in both sectors depends on the way how the limit is taken. The solution to the Friedmann equation is given by

\displaystyle a(t)=c_{(g)}(-t)^{\frac{2}{3}},

(40)

where we have supposed that $t<0$ (and $a_{,t}<0$ ) during collapse, which ends when $t=0$ , and made the identification

\displaystyle\rho_{(g,0)}=\frac{4c_{(g)}^{3}}{3\kappa^{2}}.

(41)

The solutions for the shift vectors (33) are then given by [35]

\displaystyle N^{r(-)}(t,r)=N_{(f)}^{r(-)}(t,r)=\frac{2r}{3(-t)}.

(42)

The first junction condition (26) for (30) and (33) determines the conditions on the trajectories of the interfaces

\displaystyle r_{(g)}=\frac{4R_{(g)}^{3}}{9(-T_{(g)})^{2}},\qquad r_{(f)}=\frac{4R_{(f)}^{3}}{9(-T_{(f)})^{2}},

(43)

and relate trajectories of the matter interfaces with the external Schwarzschild radii. From Eq. (18), the four velocities of the matter interfaces are given by

\displaystyle u_{(g)}^{\mu}=\left(1,-\sqrt{\frac{r_{(g)}}{R_{(g)}}},0,0\right),\qquad u_{(f)}^{\mu}=\frac{1}{bC_{0}}\left(1,-\sqrt{\frac{r_{(f)}}{R_{(f)}}},0,0\right).

(44)

From Eqs. (38) (40), and (43), we find that $u_{(g)}^{\mu}=v_{(g)}^{\mu}$ and $u_{(f)}^{\mu}=v_{(f)}^{\mu}$ at the interfaces, showing that the interfaces are comoving with geodesics of pressureless dust in both the sectors.

This result also shows that $\dot{R}_{(g)}=-\sqrt{r_{(g)}/R_{(g)}}$ , which implies that collapse starts from infinity with zero velocity. This seems to imply that we are dealing with a star with an infinite initial radius. However, we can interpret this result for a finite-radius star as follows: the model here describes that collapse has already started before the time $t_{i}\leq t(<0)$ , after which we consider the model to be an approximate description of real stellar collapse. In this case, our approximate model describes collapse in the time region $t_{i}\leq t<0$ . We should notice that, at the beginning of collapse, we have the following condition hold true: $\dot{R}_{(g)}(\tau_{(g)}=\tau_{(g)i})=-\sqrt{r_{(g)}/R_{(g)}(\tau_{(g)i})}$ . This corresponds to collapse which satisfies some particular initial conditions. This will imply that the description in this section will qualitatively/approximately describe collapse, which will have $\dot{R}_{(g)}(\tau_{(g)}=\tau_{(g)i})\approx-\sqrt{r_{(g)}/R_{(g)}(\tau_{(g)i})}$ , at least as long as we are not too close to the final singularity.

On the other hand, when $\dot{R}_{(g)}(\tau_{(g)}=\tau_{(g)i})$ will not be well approximated by the previous condition, then we should expect deviations from the description of collapse in this section. It is not clear what we should expect in the latter case, especially because the other known solutions for collapse, as we will see later on, will also be (even more) fine-tuned. We defer discussion on possible yet-unknown collapse solutions to the discussion of the spatially closed case (Sec. III.2). Furthermore, an analogous tuned initial condition also holds in the fiducial sector for $\dot{R}_{(f)}(\tau_{(f)}=\tau_{(f)i})$ , see Eq. (44), since $C_{0}$ is defined by the parameters in the MTBG theory, whereas $b$ is set by Eq. (39). Nonetheless, fine-tuned initial conditions for collapse discussed here do not impose constraints on the matter content in them. For instance, no fine-tuned condition has to be imposed on masses of black holes in both the physical and fiducial sectors, namely, $r_{(g)}\neq r_{(f)}$ in general, as will happen, in contrast, for the spatially closed case. Because of this, we believe these spatially flat collapse solutions to be less problematic than the spatially closed one, which will be discussed in Sec. III.2.

We can also show that the nontrivial components of the extrinsic curvature tensors at the interfaces are given by

\displaystyle\tilde{K}^{(g)(\pm)}_{\tau\tau}=0,\quad\tilde{K}^{(g)(\pm)}_{ab}=R_{(g)}\theta_{ab},\quad\tilde{K}^{(f)(\pm)}_{\tau\tau}=0,\quad\tilde{K}^{(f)(\pm)}_{ab}=C_{0}R_{(f)}\theta_{ab},

(45)

respectively. Thus, by imposing the continuity of the extrinsic curvature tensors $[\tilde{K}^{(g)}_{\tau\tau}]=[\tilde{K}^{(g)}_{ab}]=[\tilde{K}^{(f)}_{\tau\tau}]=[\tilde{K}^{(f)}_{ab}]=0$ , we can ensure that the surface energy-momentum tensors at the interfaces vanish,

\displaystyle S^{(g)}_{\tau\tau}=S^{(g)}_{ab}=S^{(f)}_{\tau\tau}=S^{(f)}_{ab}=0,

(46)

and the exterior Schwarzschild regions (30) and interior dust-dominated FLRW universes (33) are smoothly joined across the interfaces [34, 36]. In the case where the interior regions are filled with matter fields with nonvanishing pressures, the propagation speed of sound waves is nonvanishing and thus nearby fluid elements do not evolve independently. In this case, in order to study the interior solution, one needs to solve a set of partial differential equations, instead of ordinary differential equations.

Before closing this subsection, we should mention that it is straightforward to extend the collapsing solution discussed in this subsection (and the one in next subsection, i.e., the Oppenheimer-Snyder model [34, 36]) in the presence of the nonzero values of the effective cosmological constants [see Eq. (99)] as in the case of GR [see e.g., Refs. [38]], where the exterior spacetimes are asymptotically de Sitter with the effective cosmological constants (99). However, for gravitational collapse of astrophysical interest, the effect of the effective cosmological constants would be negligible. Thus, for our purpose, it is sufficient to focus on the case of the asymptotically flat exterior solutions.

III.2 Collapse in the spherically symmetric spacetimes in the spatially closed coordinates

We consider the spherical collapse model of a pressureless dust where the spacetimes inside the matter surfaces are described by the spatially closed FLRW universes. In contrast to the spatially flat case, in the spatially closed case gravitational collapse starts from a finite distance. We note that the spatially closed case is an extension of the Oppenheimer-Snyder model in GR [34] to the case of MTBG.

III.2.1 The exterior solution

In order to discuss gravitational collapse of the spatially closed matter surfaces in both the physical and fiducial sectors, we assume that the spacetime metrics in the exterior regions are described by the Schwarzschild solutions written in the spatially closed coordinates

	$\displaystyle g_{\mu\nu}^{(+)}dx^{\mu}dx^{\nu}=-dt^{2}+\frac{1}{1-q_{0,(g)}}\left(dr+\sqrt{\frac{r_{(g)}}{r}-q_{0,(g)}}dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b},$		(47)
	$\displaystyle f_{\mu\nu}^{(+)}dx^{\mu}dx^{\nu}=C_{0}^{2}\left[-b^{2}dt^{2}+\frac{1}{1-q_{0,(f)}}\left(dr+\sqrt{\frac{r_{(f)}}{r}-Qq_{0,(f)}}dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b}\right].$		(48)

where $q_{0,(g)}$ and $q_{0,(f)}$ are constants satisfying

	$\displaystyle-1<-\frac{r_{(g)}}{r}<-q_{0,(g)}<0,$		(49)
	$\displaystyle-1<-\frac{r_{(f)}}{Qr}<-q_{0,(f)}<0,$		(50)

and $Q>0$ is also a constant. We also assume that the Lagrange multipliers are given by Eq. (25).

In the self-accelerating branch, the equations for $\lambda$ and ${\bar{\lambda}}$ are satisfied if

\displaystyle q_{0,(f)}=q_{0,(g)},

(51)

as well as the condition (29). The constraint equation for $\lambda^{r}$ is automatically satisfied. The equations of motion for the metric variables finally reduce to Eqs. (31) and (32), and

\displaystyle q_{0,(g)}{\tilde{\alpha}}^{2}\left(Q-b^{2}\right)=0.

(52)

Thus, for the spatially closed coordinates $q_{0,(g)}\neq 0$ , we obtain

\displaystyle Q=b^{2}.

(53)

Summarizing the above, the exterior metrics (47) and (48) reduce to

	$\displaystyle g_{\mu\nu}^{(+)}dx^{\mu}dx^{\nu}=-dt^{2}+\frac{1}{1-q_{0,(g)}}\left(dr+\sqrt{\frac{r_{(g)}}{r}-q_{0,(g)}}dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b},$		(54)
	$\displaystyle f_{\mu\nu}^{(+)}dx^{\mu}dx^{\nu}=C_{0}^{2}\left[-b^{2}dt^{2}+\frac{1}{1-q_{0,(g)}}\left(dr+\sqrt{\frac{r_{(f)}}{r}-b^{2}q_{0,(g)}}dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b}\right].$		(55)

III.2.2 The interior solution

We then consider spatially closed FLRW metrics as the interior solution [36, 35], written in comoving coordinates $(t_{c},r_{c})$ , which can be written as

	$\displaystyle g^{(-)}_{\mu\nu}dx^{\mu}dx^{\nu}$	$\displaystyle=-dt_{c}^{2}+a(t)^{2}\left[\frac{dr_{c}^{2}}{1-r_{c}^{2}}+r_{c}^{2}\theta_{ab}d\theta^{a}d\theta^{b}\right],$		(56)
	$\displaystyle f^{(-)}_{\mu\nu}dx^{\mu}dx^{\nu}$	$\displaystyle=-C_{0}^{2}{\tilde{b}}(t)^{2}dt_{c}^{2}+a_{f}(t)^{2}\left[\frac{dr_{c}^{2}}{1-r_{c}^{2}}+r_{c}^{2}\theta_{ab}d\theta^{a}d\theta^{b}\right],$		(57)

and the time-reparametrization invariance allows us to change variables into $t_{c}=t$ and $r_{c}=\frac{r}{a(t)}$ , as to bring the metric in the form similar to the one of GP coordinates as

	$\displaystyle g_{\mu\nu}^{(-)}dx^{\mu}dx^{\nu}=-dt^{2}+\frac{1}{1-\frac{r^{2}}{a(t)^{2}}}\left(dr-r\frac{\dot{a}(t)}{a(t)}dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b},$		(58)
	$\displaystyle f_{\mu\nu}^{(-)}dx^{\mu}dx^{\nu}=C_{0}^{2}\left[-\tilde{b}(t)^{2}dt^{2}+\frac{1}{1-\frac{r^{2}}{a(t)^{2}}}\left(dr-r\frac{\dot{a}(t)}{a(t)}dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b}\right],$		(59)

where as for the spatially flat case, in MTBG, for the self-accelerating branch, the constraints impose that

a_{f}=C_{0}\,a.

(60)

In the self-accelerating branch, the constant $C_{0}$ is uniquely determined by the parameters of the Lagrangian, as Eq (29). In any case, either in comoving coordinates $(t_{c},r_{c})$ or in GP-like coordinates $(t,r)$ , the Friedmann equations remain the same. In particular, the (constant) proportionality between the two scale factors leads to the following result

H_{f}:=\frac{a_{f,t}}{C_{0}\tilde{b}(t)a_{f}}=\frac{a_{,t}}{C_{0}\tilde{b}(t)a}=\frac{H}{C_{0}\tilde{b}(t)}\,\qquad{\rm or}\qquad\tilde{b}(t)=\frac{H}{C_{0}H_{f}}\,,

(61)

where a $a_{,t}\equiv da/dt$ , etc. If we use the Friedmann equations of motion for both the metrics, we find

	$\displaystyle\frac{3(1+{a}_{,t}^{2})}{a^{2}}=\kappa^{2}\frac{\rho_{(g,0)}}{a^{3}},\qquad\frac{3(\tilde{b}^{2}+{a}_{,t}^{2})}{\tilde{b}^{2}a^{2}}=\frac{C_{0}^{2}\kappa^{2}}{{\tilde{\alpha}}^{2}}\frac{\rho_{(f,0)}}{a^{3}},$
	$\displaystyle 1+a_{,t}^{2}+2aa_{,tt}=0,\qquad\tilde{b}^{3}-2aa_{,t}{\tilde{b}}_{,t}+\tilde{b}(a_{,t}^{2}+2aa_{,tt})=0,$		(62)

so that $\tilde{b}(t)$ is, in general, a function of time, which can be written as

\tilde{b}(t)=\frac{1}{C_{0}}\sqrt{\frac{\frac{\kappa^{2}\rho_{(g,0)}}{a^{3}}-\frac{3}{a^{2}}}{\frac{\kappa^{2}\rho_{(f,0)}}{\tilde{\alpha}^{2}a^{3}}-\frac{3}{C_{0}^{2}a^{2}}}}.

(63)

Since ${\tilde{b}}=H/(H_{f}C_{0})$ , if we want to have a finite (and nonzero) value of ${\tilde{b}}$ when $H$ (or $H_{f}$ ) vanishes, one needs to require that the initial time of the collapse, usually set at $\dot{a}=0$ , is the same for both the metrics. This initial starting-from-rest collapse time is then defined as $a_{,t}=0$ (or $a_{f,t}=0$ ) and $a=a_{{\rm max}}\neq 0$ . In this case, we have

\kappa^{2}\rho_{(g,0)}-3a_{{\rm max}}=0=\kappa^{2}\frac{C_{0}^{2}\rho_{(f,0)}}{\tilde{\alpha}^{2}}-3a_{{\rm max}}\,,

(64)

indicating that

\displaystyle\rho_{(f,0)}=\frac{\tilde{\alpha}^{2}}{C_{0}^{2}}\rho_{(g,0)}\,,

(65)

which, in general, corresponds to a fine-tuned value for the matter densities in the fiducial sector with respect to the one in the physical sector, as $C_{0}$ and $\tilde{\alpha}$ are specified by the parameters of the theory, as already stated above. In this case, we obtain

\tilde{b}(t)=\frac{1}{C_{0}}\sqrt{\frac{\frac{\kappa^{2}\rho_{(g,0)}}{a^{3}}-\frac{3}{a^{2}}}{\frac{\kappa^{2}\rho_{(f,0)}}{\tilde{\alpha}^{2}a^{3}}-\frac{3}{C_{0}^{2}a^{2}}}}=1\,.

(66)

Thus, the value of ${\tilde{b}}$ is no longer a function of time but actually equal to unity. On the other hand, we paid the price of a problematic fine-tuning condition among matter contents in the different sectors, namely, Eq. (65), where, at least in the self-accelerating branch, the value $\tilde{\alpha}/C_{0}$ corresponds to a constant in the theory. In particular, as was happening in the spatially flat case, matter should be present in both the sectors.

The exact solution of Eq. (III.2.2) is given by the parametrization

\displaystyle t(\eta)

\displaystyle=

\displaystyle\frac{\kappa^{2}}{6}\rho_{0,(g)}\left(\eta+\sin\eta\right),\qquad a(\eta)=\frac{\kappa^{2}}{6}\rho_{0,(g)}\left(1+\cos\eta\right),\qquad{\tilde{b}}(\eta)=\frac{\tilde{\alpha}\sqrt{\rho_{0,(g)}(1-\cos\eta)}}{\sqrt{2C_{0}^{2}\rho_{0,(f)}-\rho_{0,(g)}{\tilde{\alpha}}^{2}(1+\cos\eta)}},

(67)

where $\eta$ corresponds to the conformal time $dt=a(\eta)d\eta$ . Since the exterior $b$ is constant, the interior $\tilde{b}$ has to be constant. Imposing Eq. (65), we also obtain $\tilde{b}=1$ .

This situation is much worse than that in the spatially flat case, in which we only needed to require that $\rho_{(g,0)}>0$ and $\rho_{(f,0)}>0$ . This result of the spatially closed collapse might merely indicate that a general configuration of matter in both the sectors will not follow the kind of unique collapse described in this section. How to avoid this situation? One possibility, still keeping the spatially closed picture, would be that, similar to what we have encountered in the spatially flat collapse, the spatially closed collapse does not start from a configuration where $H$ and/or $H_{f}$ vanish, but collapse has already started some time before the particular time $t_{i}$ we start describing its dynamics by means of this model. Then, both $H$ and $H_{f}$ do not vanish and are both negative for $t_{i}\leq t$ . This possibility is viable, however, it leads, in general, to a time-dependence of the function $\tilde{b}(t)$ , which cannot in general be matched with a stationary configuration for the known exterior solutions of MTBG. In fact, in MTBG, the Birkhoff theorem does not hold in general and, at least in principle, we are bound to look for all the possible spherically symmetric solutions allowed by the theory. Therefore, there could be still unknown, in general time-dependent solutions that could accommodate a matching with the $\tilde{b}(t)$ imposed by Eq. (66).

In any case, in this section, we will only discuss $\tilde{b}(t)=1$ and describe what this fine-tuned configuration leads to. Hence, on using the GP-like coordinates described already for the spatially flat case, we can consider the form of the two interior metrics given as follows:

	$\displaystyle g_{\mu\nu}^{(-)}dx^{\mu}dx^{\nu}=-dt^{2}+\frac{1}{1-\frac{r^{2}}{a(t)^{2}}}\left(dr-r\frac{\dot{a}(t)}{a(t)}dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b},$		(68)
	$\displaystyle f_{\mu\nu}^{(-)}dx^{\mu}dx^{\nu}=C_{0}^{2}\left[-dt^{2}+\frac{1}{1-\frac{r^{2}}{a_{f}(t)^{2}}}\left(dr-r\frac{\dot{a}_{f}(t)}{a_{f}(t)}dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b}\right].$		(69)

With the background equation (III.2.2), the interior metrics (68) and (69) reduce to

	$\displaystyle g_{\mu\nu}^{(-)}dx^{\mu}dx^{\nu}$	$\displaystyle=$	$\displaystyle-dt^{2}+\frac{1}{1-\frac{r^{2}}{a(t)^{2}}}\left(dr+\sqrt{\frac{\kappa^{2}}{3}\frac{\rho_{0,(g)}}{a^{3}}r^{2}-\frac{r^{2}}{a(t)^{2}}}dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b},$		(70)
	$\displaystyle f_{\mu\nu}^{(-)}dx^{\mu}dx^{\nu}$	$\displaystyle=$	$\displaystyle C_{0}^{2}\left[-dt^{2}+\frac{1}{1-\frac{r^{2}}{a(t)^{2}}}\left(dr+\sqrt{\frac{C_{0}^{2}\kappa^{2}}{3{\tilde{\alpha}}^{2}}\frac{\rho_{0,(f)}}{a^{3}}r^{2}-\frac{r^{2}}{a(t)^{2}}}dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b}\right].$		(71)

III.2.3 Matching conditions at the interfaces

Matching of the temporal component of the exterior metric (55) and interior metric (71) at the matter interfaces yields

\displaystyle b={\tilde{b}}=1.

(72)

We consider the exterior metrics as in

	$\displaystyle g^{(+)}_{\mu\nu}dx^{\mu}dx^{\nu}$	$\displaystyle=-dt^{2}+\frac{1}{1-q_{0,(g)}}\left(dr+\xi\,dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b}\,,$		(73)
	$\displaystyle f^{(+)}_{\mu\nu}dx^{\mu}dx^{\nu}$	$\displaystyle=C_{0}^{2}\left[-dt^{2}+\frac{1}{1-q_{0,(g)}}\left(dr+\xi_{2}\,dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b}\right]\,,$		(74)

where we have also defined $N^{r}=\xi:=\sqrt{\frac{r_{(g)}}{r}-q_{0,(g)}}$ and $N^{r}_{(f)}=\xi_{2}:=\sqrt{\frac{r_{(f)}}{r}-q_{0,(g)}}$ . This ansatz satisfies the MTBG equations of motion for the vanishing effective cosmological constants in both the sectors that we have neglected in these physical collapse scenarios. Let us now consider the four velocities for both the sectors, at the interface, as given by

\displaystyle u_{(g)}^{\mu}=\left(\dot{T}_{(g)},\dot{R}_{(g)},0,0\right),\qquad u_{(f)}^{\mu}=\left(\dot{T}_{(f)},\dot{R}_{(f)},0,0\right)\,,

(75)

where a dot, in this subsection, represents a derivative with respect to the proper times $\tau_{(g)}$ (or $\tau_{(f)}$ ), e.g. $\dot{R}_{(g)}(\tau_{(g)})=R_{(g),\tau}=dR_{(g)}/d\tau_{(g)}$ . Their normalization conditions, $g_{\mu\nu}u_{(g)}^{\mu}u_{(g)}^{\nu}=-1$ and $f_{\mu\nu}u_{(f)}^{\mu}u_{(f)}^{\nu}=-1$ , give

	$\displaystyle\dot{T}_{(g)}$	$\displaystyle=\frac{\xi\dot{R}_{(g)}+\sqrt{(1-q_{0,(g)})(\dot{R}_{(g)}^{2}-\xi^{2}+1-q_{0,(g)})}}{1-\xi^{2}-q_{0,(g)}}>0\,,$		(76)
	$\displaystyle\dot{T}_{(f)}$	$\displaystyle=\frac{\xi_{2}\dot{R}_{(f)}C_{0}+\sqrt{(1-q_{0,(g)})(C_{0}^{2}\dot{R}_{(f)}^{2}-\xi_{2}^{2}+1-q_{0,(g)})}}{\left(1-\xi_{2}^{2}-q_{0,(g)}\right)C_{0}}>0\,,$		(77)

whereas we consider $\dot{R}_{(g)}<0$ , $\dot{R}_{(f)}<0$ . We will make use of these equations later on again.

Next we define the normal to the interface as $n_{\mu}^{(g)}dx^{\mu}\propto(dr-R_{(g),t}\,dt)$ , which is normalized to $g_{\mu\nu}n_{(g)}^{\mu}n_{(g)}^{\nu}=1$ , and the interface equation is described by $r-R_{(g)}(t)=0$ . Because of $R_{(g),t}=dR_{(g)}/dt=\dot{R}_{(g)}/\dot{T}_{(g)}$ , one can see that $n_{\mu}^{(g)}u_{(g)}^{\mu}=0$ on the interface. Along the same lines, we can introduce $n_{\mu}^{(f)}dx^{\mu}\propto(dr-R_{(f),t}\,dt)$ , which is now normalized to $f_{\mu\nu}n_{(f)}^{\mu}n_{(f)}^{\nu}=1$ , and the relative interface equation is $r-R_{(f)}(t)=0$ . In this case we have that $R_{(f),t}=dR_{(f)}/dt=\dot{R}_{(f)}/\dot{T}_{(f)}$ .

For the physical metric, we define the projection tensors to the interface as $e_{(g)}^{\mu}{}_{\tau}\partial_{\mu}=\dot{T}_{(g)}\partial_{t}+\dot{R}_{(g)}\partial_{r}$ and $e_{(g)}^{\mu}{}_{a}\partial_{\mu}=\partial_{a}$ , satisfying $n_{\mu}^{(g)}e_{(g)}^{\mu}{}_{i}=0$ , ( $i=\tau,a$ ), and calculate $h_{\mu\nu}^{(g)}=g_{\mu\nu}-n_{\mu}^{(g)}n_{\nu}^{(g)}$ , together with $K_{\mu\nu}^{(g)}=h_{\mu}^{(g)}{}^{\alpha}h_{\nu}^{(g)}{}^{\beta}n_{\alpha;\beta}^{(g)}$ . Finally, we can build $h_{ab}^{(g)}=\left.e_{(g)}^{\mu}{}_{a}e_{(g)}^{\nu}{}_{b}h_{\mu\nu}^{(g)}\right|_{t=T_{(g)},r=R_{(g)}}$ and $K_{ab}^{(g)}=\left.e_{(g)}^{\mu}{}_{a}e_{(g)}^{\nu}{}_{b}K_{\mu\nu}^{(g)}\right|_{t=T_{(g)},r=R_{(g)}}$ . At the interface, the nontrivial components to consider for the matching conditions are the following ones:

	$\displaystyle h_{\tau\tau}^{(g)}$	$\displaystyle=-1,\qquad h_{ab}^{(g)}=R_{(g)}^{2}\theta_{ab}\,,$		(78)
	$\displaystyle K_{\tau\tau}^{(g)}$	$\displaystyle=\mathcal{K}(R_{(g)}\,,\dot{R}_{(g)},\ddot{R}_{(g)}),\qquad K_{ab}^{(g)}=R_{(g)}\theta_{ab}\sqrt{\dot{R}_{(g)}^{2}-\xi^{2}+1-q_{0,(g)}}\,,$		(79)

where in the expression of $\mathcal{K}$ , the values of $\dot{T}_{(g)}$ and $\ddot{T}_{(g)}$ , present because $\mathcal{K}$ depends also on $R_{(g),tt}$ and $R_{(g),t}$ , have been replaced by using Eq. (76).

For the fiducial metric, we can follow an analogous path, namely defining the projectors $e_{(f)}^{\mu}{}_{\tau}\partial_{\mu}=\dot{T}_{(f)}\partial_{t}+\dot{R}_{(f)}\partial_{r}$ and $e_{(f)}^{\mu}{}_{a}\partial_{\mu}=\partial_{a}$ , satisfying $n_{\mu}^{(f)}e_{(f)}^{\mu}{}_{i}=0$ , ( $i=\tau,a$ ). After calculating the expressions for $h_{\mu\nu}^{(f)}=f_{\mu\nu}-n_{\mu}^{(f)}n_{\nu}^{(f)}$ , and $K_{\mu\nu}^{(f)}=h_{\mu}^{(f)}{}^{\alpha}h_{\nu}^{(f)}{}^{\beta}n_{\alpha;\beta}^{(f)}$ , one can build $h_{ab}^{(f)}=\left.e_{(f)}^{\mu}{}_{a}e_{(f)}^{\nu}{}_{b}h_{\mu\nu}^{(f)}\right|_{t=T_{(f)},r=R_{(f)}}$ and $K_{ab}^{(f)}=\left.e_{(f)}^{\mu}{}_{a}e_{(f)}^{\nu}{}_{b}K_{\mu\nu}^{(f)}\right|_{t=T_{(f)},r=R_{(f)}}$ evaluating them on the interface. Then we find

	$\displaystyle h_{\tau\tau}^{(f)}$	$\displaystyle=-1,\qquad h_{ab}^{(f)}=C_{0}^{2}R_{(f)}^{2}\theta_{ab}\,,$		(80)
	$\displaystyle K_{\tau\tau}^{(f)}$	$\displaystyle=\mathcal{K}_{f}(R_{(f)}\,,\dot{R}_{(f)},\ddot{R}_{(f)})\,,\qquad K_{ab}^{(f)}=C_{0}R_{(f)}\theta_{ab}\sqrt{C_{0}^{2}\dot{R}_{(f)}^{2}-\xi_{2}^{2}+1-q_{0,(g)}}\,.$		(81)

Notice that, at this point, we have still to set the junction conditions to match the interior/exterior solutions.

In fact, as for the internal metrics, we have

	$\displaystyle g^{(-)}_{\mu\nu}dx^{\mu}dx^{\nu}$	$\displaystyle=-dt^{2}+\frac{1}{1-\frac{r^{2}}{a^{2}}}\left(dr-r\frac{a_{,t}}{a}\,dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b}\,,$		(82)
	$\displaystyle f^{(-)}_{\mu\nu}dx^{\mu}dx^{\nu}$	$\displaystyle=C_{0}^{2}\left[-dt^{2}+\frac{1}{1-\frac{r^{2}}{a^{2}}}\left(dr-r\frac{a_{,t}}{a}\,dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b}\right],$		(83)

where $a_{,t}=da/dt$ . These expressions describe the case of a pair of spatially closed FLRW universes. For the case of $b=1$ , as we have seen above, $\rho_{(g,0)}$ and $\rho_{(f,0)}$ are proportional to each other, so that the two Friedmann equations are equivalent to each other, giving effectively only one single constraint, which can be written as Eq. (III.2.2). For these metrics, the four velocities of matter are given by³³3To see this, we can perform a coordinate change as $t_{c}=t$ and $r_{c}=r/a$ , that is, $dt_{c}=dt$ and $dr_{c}=dr/a-r\,(a_{,t}/a^{2})dt$ . Then in comoving coordinates $v_{(g)}^{\mu}\partial_{\mu}=\partial_{t_{c}}$ . However, we can see that $\partial_{t_{c}}=\partial_{t}+ra_{,t}/a\,\partial_{r}$ (and $\partial_{r_{c}}=a\partial_{r}$ ), giving the result in Eq. (84).

	$\displaystyle v_{(g)}^{\mu}\partial_{\mu}$	$\displaystyle=\partial_{t}+r\frac{a_{,t}}{a}\,\partial_{r}\,,$		(84)
	$\displaystyle v_{(f)}^{\mu}\partial_{\mu}$	$\displaystyle=\frac{1}{C_{0}}\left(\partial_{t}+r\frac{a_{,t}}{a}\,\partial_{r}\right).$		(85)

We also define the normalized normal vector as $n_{\mu}^{(g)}\propto(dr-r_{c0}a_{,t}dt)$ , whereas $n_{\mu}^{(f)}\propto(dr-r_{c0f}a_{,t}dt)$ . Then $v_{(g)}^{\mu}n_{\mu}^{(g)}=0$ on the interface defined by $r=r_{c0}a$ , where matter is moving on constant comoving coordinate $r_{c}$ , and similarly by $r=r_{c0f}a$ in the fiducial sector. Choosing $e_{(g)}^{\mu}{}_{\tau}\partial_{\mu}=\partial_{t}+r\frac{a_{,t}}{a}\,\partial_{r}$ and $e_{(g)}^{\mu}{}_{a}\partial_{\mu}=\partial_{a}$ , we can define $h_{ab}^{(g)}$ and $K_{ab}^{(g)}$ as was done for the exterior metrics. Then we find

	$\displaystyle h_{\tau\tau}^{(g)}$	$\displaystyle=-1,\qquad h_{ab}^{(g)}=r_{c0}^{2}a^{2}\theta_{ab}\,,$		(86)
	$\displaystyle K_{\tau\tau}^{(g)}$	$\displaystyle=0,\qquad K_{ab}^{(g)}=r_{c0}a\theta_{ab}\sqrt{1-r_{c0}^{2}}\,.$		(87)

We then discuss the junction conditions for the physical metric first. In this case, we first find, from the $h_{ab}^{(g)}$ expressions in Eqs. (78) and (86), that $R_{(g)}=r_{c0}a\,$ . Then the components of $K_{ab}^{(g)}$ given in Eqs. (79) and (87) lead to $\dot{R}_{(g)}^{2}-\xi(R_{(g)})^{2}-q_{0,(g)}+r_{c0}^{2}=0$ or

\dot{R}_{(g)}=-\sqrt{\xi(R_{(g)})^{2}+q_{0,(g)}-r_{c0}^{2}}.

(88)

Using this relation, we find that $K_{\tau\tau}^{(g)}=\mathcal{K}=0$ also for the exterior solution.

Since we want $u_{(g)}^{\mu}$ to coincide at the interface with $v_{(g)}^{\mu}$ , then, from their zeroth component, we need to set $\dot{T}_{(g)}=1$ . Using this last condition, together with the expression for $\dot{R}_{(g)}$ given in Eq. (88) into Eq. (76), we find that the condition $q_{0,(g)}=r_{c0}^{2}$ needs to be imposed. Furthermore, we have

\frac{da}{dt}=a_{,t}=\frac{1}{r_{c0}}\,R_{(g),t}=\frac{1}{r_{c0}}\,\frac{\dot{R}_{(g)}}{\dot{T}_{(g)}}=\frac{1}{r_{c0}}\,\dot{R}_{(g)}\,,

(89)

which, once inserted into the Friedmann equation, together with the expression for $\dot{R}_{(g)}$ given in Eq. (88), makes the Friedmann equation impose that $(\rho_{(g,0)}r_{c0}^{3}-3M_{P}^{2}r_{(g)})/R_{(g)}^{3}=0$ and sets the value of $r_{c0}$ in terms of the other parameters. Furthermore, we can also see that, on the interface, $v_{(g)}^{r}=R_{(g)}\frac{a_{,t}}{a}=R_{(g),t}=\dot{R}_{(g)}=u_{(g)}^{r}$ , which holds true since $\dot{T}_{(g)}=1$ .

Let us now focus on the fiducial sector. For the fiducial metric, we choose $e_{(f)}^{\mu}{}_{\tau}\partial_{\mu}=C_{0}^{-1}[\partial_{t}+r\frac{a_{,t}}{a}\,\partial_{r}]$ and $e_{(f)}^{\mu}{}_{a}\partial_{\mu}=\partial_{a}$ as the projection tensors. On using the expression for $n_{(f)}^{\mu}$ , we can define $h_{ab}^{(f)}$ and $K_{ab}^{(f)}$ as done for the exterior metrics. The results of these calculations can be written as

	$\displaystyle h_{\tau\tau}^{(f)}$	$\displaystyle=-1,\qquad h_{ab}^{(f)}=C_{0}^{2}r_{c0f}^{2}a^{2}\theta_{ab}\,,$		(90)
	$\displaystyle K_{\tau\tau}^{(f)}$	$\displaystyle=0\,\qquad K_{ab}^{(f)}=C_{0}r_{c0f}a\theta_{ab}\sqrt{1-r_{c0f}^{2}}\,,$		(91)

which, on using the junction conditions for $h_{ab}^{(f)}$ and $K_{ab}^{(f)}$ , lead to

	$\displaystyle R_{(f)}$	$\displaystyle=r_{c0f}a\,,$		(92)
	$\displaystyle\dot{R}_{(f)}$	$\displaystyle=-C_{0}^{-1}\sqrt{\xi_{2}(R_{(f)})^{2}+q_{0,(g)}-r_{c0f}^{2}}\,.$		(93)

This relation, once inserted into Eq. (81), gives $K_{\tau\tau}^{(f)}=0$ for the exterior metric.

Following exactly the same procedure taken with the physical metric, namely by setting $u_{(f)}^{\mu}$ to coincide with $v_{(f)}^{\mu}$ at the interface, ⁴⁴4As for the $r$ -component, we have $u_{(f)}^{r}=\dot{R}_{(f)}$ , whereas $v_{(f)}^{r}=\left.\frac{r}{C_{0}}\frac{a_{,t}}{a}\right|_{t=T_{(f)},r=R_{(f)}}=\frac{R_{(f),t}}{C_{0}}=\frac{\dot{R}_{(f)}}{C_{0}\dot{T}_{f}}=\dot{R}_{(f)}$ , so that all the components of $u_{(f)}^{\mu}=v_{(f)}^{\mu}$ coincide on the interface. we need to set $\dot{T}_{(f)}=C_{0}^{-1}$ , which, in turn, leads to $q_{0,(g)}=r_{c0f}^{2}$ from Eqs. (77) and (93). Then on using the result that $r_{c0}^{2}=q_{0,(g)}=r_{c0f}^{2}$ , we can see that $3r_{(g)}=\kappa^{2}\rho_{(g,0)}r_{c0}^{3}=\kappa^{2}\rho_{(g,0)}r_{c0f}^{3}=3r_{(f)}$ from the Friedmann equation given in Eq. (III.2.2). Then we are bound to conclude that $r_{(f)}=r_{(g)}$ . Therefore, the fine-tuned interior matter content leads to fine-tuned values for the Schwarzschild radii for the exterior metrics.

III.3 Absence of cusps of the constant-time hypersurfaces

Before closing this section, we check the continuity of the normal vectors to the $t={\rm constant}$ hypersurfaces across the matter interfaces in both the spatially flat and closed cases and in both the physical and fiducial sectors, i.e., the continuity of the derivatives of the Stückelberg fields along the normals across the matter interfaces. Because of the absence of the four-dimensional diffeomorphism invariance, in general, different time slicings would provide different solutions, and hence this condition is not obvious and has to be checked for each solution. If this condition fails to be satisfied, cusps would be formed on the interfaces, which would give rise to extra boundary contributions to the evolution of the metric and matter fields. In our cases, the nontrivial components of the normals to a $t={\rm constant}$ hypersurface

\displaystyle m_{\mu}^{(g)}dx^{\mu}=-dt,\qquad m_{\mu}^{(f)}dx^{\mu}=-C_{0}bdt,

(94)

where $b=1$ for the spatially closed case. We then find that at the matter interfaces

	$\displaystyle n_{\alpha}^{(g)}m^{\alpha}_{(g)}\Big{\|}_{(T_{(g)},R_{(g)})}=0,\qquad n_{\alpha}^{(f)}m^{\alpha}_{(f)}\Big{\|}_{(T_{(f)},R_{(f)})}=0,$
	$\displaystyle e^{\alpha}_{(g)\tau}m_{\alpha}^{(g)}\Big{\|}_{(T_{(g)},R_{(g)})}=-1,\qquad e^{\alpha}_{(f)\tau}m_{\alpha}^{(f)}\Big{\|}_{(T_{(f)},R_{(f)})}=-1,$
	$\displaystyle e^{\alpha}_{(g)a}m_{\alpha}^{(g)}\Big{\|}_{(T_{(g)},R_{(g)})}=0,\qquad e^{\alpha}_{(f)a}m_{\alpha}^{(f)}\Big{\|}_{(T_{(f)},R_{(f)})}=0.$		(95)

Hence, the normals to the $t={\rm constant}$ hypersurfaces are continuous across the matter interfaces in both the sectors. This implies that there is no cusp formation in the $t={\rm constant}$ hypersurfaces at the matter interfaces.

In summary, the results in this section suggest that, in the self-accelerating branch of MTBG with matter composed of two independent species (13), gravitational collapse proceeds as that in the two copies of GR. Thus, in order to observe the deviation from the case of two copies of GR, we should investigate the linear perturbations about the spherically symmetric solutions. In the next section, we will study the odd-parity perturbations of the Schwarzschild $-$ de Sitter solutions written in the spatially flat coordinates in the self-accelerating branch of MTBG.

IV Odd-parity perturbations of Schwarzschild $-$ de Sitter solutions in MTBG

In this section, we investigate the linear perturbations of the Schwarzschild $-$ de Sitter solutions written in the spatially flat coordinates in the self-accelerating branch of MTBG. The linear perturbations in the static and spherically symmetric black holes are decomposed into the odd- and even-parity sectors. In this paper, we focus on the odd-parity sector of the perturbations.

IV.1 Odd-parity black hole perturbations in the spatially flat coordinates

The perturbed physical and fiducial metrics in the odd-parity sectors about the Schwarzschild $-$ de Sitter solutions written in the spatially flat coordinates are, respectively, given by

$\displaystyle g_{\mu\nu}dx^{\mu}dx^{\nu}$	$\displaystyle=$	$\displaystyle-dt^{2}+\left(dr+N^{r}(r)dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b}$	(96)
	$\displaystyle+$	$\displaystyle 2\sum_{\ell,m}\left(r^{2}H_{t}(t,r)dt+rH_{r}(t,r)(dr+N^{r}(r)dt)\right)E_{a}{}^{b}\partial_{b}Y_{\ell m}(\theta^{c})d\theta^{a}$
	$\displaystyle+$	$\displaystyle r^{2}\sum_{\ell,m}H_{3}(t,r)E_{(a}{}^{c}{\tilde{\nabla}}_{b)}{\tilde{\nabla}}_{c}Y_{\ell m}(\theta^{c})d\theta^{a}d\theta^{b},$
$\displaystyle f_{\mu\nu}dx^{\mu}dx^{\nu}$	$\displaystyle=$	$\displaystyle C_{0}^{2}\Big{[}-b^{2}dt^{2}+\left(dr+N_{(f)}^{r}(r)dt\right)^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b}$	(97)
	$\displaystyle+$	$\displaystyle 2\sum_{\ell,m}\left(r^{2}K_{t}(t,r)dt+rK_{r}(t,r)(dr+N_{(f)}^{r}(r)dt)\right)E_{a}{}^{b}\partial_{b}Y_{\ell m}(\theta^{c})d\theta^{a}$
	$\displaystyle+$	$\displaystyle r^{2}\sum_{\ell,m}K_{3}(t,r)E_{(a}{}^{c}{\tilde{\nabla}}_{b)}{\tilde{\nabla}}_{c}Y_{\ell m}(\theta^{c})d\theta^{a}d\theta^{b}\Big{]},$

where $N^{r}(r)$ and $N^{r}_{(f)}(r)$ are given by

\displaystyle N^{r}(r)=\sqrt{\frac{\Lambda_{(g)}r^{2}}{3}+\frac{r_{(g)}}{r}},\qquad N_{(f)}^{r}(r)=\sqrt{\frac{r^{2}C_{0}^{2}b^{2}\Lambda_{(f)}}{3}+\frac{r_{(f)}}{r}},

(98)

$b>0$ is a constant, $Y_{\ell m}(\theta^{c})$ are spherical harmonics with the angular multipole and magnetic moments with $-\ell\leq m\leq\ell$ , respectively, ${\tilde{\nabla}}_{a}$ represents the covariant derivative with respect to the unit two-sphere with the metric $\theta_{ab}$ , and $E_{ab}$ represents the totally antisymmetric tensor on the two-sphere satisfying ${\tilde{\nabla}}^{a}E_{ab}=0$ . $\Lambda_{(g)}$ and $\Lambda_{(f)}$ are the effective cosmological constants defined as $\Lambda_{(g)}=3H_{(g)}^{2}$ and $\Lambda_{(f)}=3H_{(f)}^{2}$ on a cosmological de Sitter expansion with Hubble parameters $H_{(g)}$ for the physical metric and $H_{(f)}$ for the fiducial one, whose values can be rewritten as

\displaystyle\Lambda_{(g)}=\frac{m^{2}\left(c_{4}-2C_{0}^{3}c_{1}-3C_{0}^{2}c_{2}\right)}{2},\qquad\Lambda_{(f)}=\frac{\left(C_{0}^{2}c_{0}+2C_{0}c_{1}+c_{2}\right)m^{2}}{2C_{0}^{2}{\tilde{\alpha}}^{2}}.

(99)

For the existence of the background Schwarzschild $-$ de Sitter solutions written in the spatially flat coordinates in the self-accelerating branch, the constants $c_{0}$ , $c_{1}$ , $c_{2}$ , $c_{3}$ , and $c_{4}$ satisfy Eq. (29), provided that we add to the right-hand side of Eq. (31) the quantity $2\Lambda_{(g)}/m^{2}$ and the quantity $\frac{2C_{0}^{3}\Lambda_{(f)}{\tilde{\alpha}}^{2}}{m^{2}}$ to the right-hand side of Eq. (32). This procedure gives the values written in Eq. (99). In the rest, with use of Eqs. (29) and (99), the constants $c_{0}$ , $c_{3}$ , and $c_{4}$ are related to the other coupling constants $(c_{1},c_{2})$ and the constant $C_{0}$ . We also assume that $C_{0}c_{1}+c_{2}\neq 0$ and $\beta C_{0}\neq 1$ . As noted in Sec. III.1, the constant $C_{0}b$ measures the ratio of the lapse functions between the fiducial and physical sectors. As we will see, in the case of $b=1$ the odd-parity perturbations in the two sectors effectively become massless and decoupled.

The perturbed three-dimensional metrics in the odd-parity sector are, respectively, given by

$\displaystyle\gamma_{ij}dy^{i}dy^{j}$	$\displaystyle=$	$\displaystyle dr^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b}$	(100)
	$\displaystyle+$	$\displaystyle 2r\sum_{\ell,m}H_{r}(t,r)E_{a}{}^{b}\partial_{b}Y_{\ell m}(\theta^{c})drd\theta^{a}+r^{2}\sum_{\ell,m}H_{3}(t,r)E_{(a}{}^{c}{\tilde{\nabla}}_{b)}{\tilde{\nabla}}_{c}Y_{\ell m}(\theta^{c})d\theta^{a}d\theta^{b},$	(100)
$\displaystyle\phi_{ij}dy^{i}dy^{j}$	$\displaystyle=$	$\displaystyle C_{0}^{2}\Big{[}dr^{2}+r^{2}\theta_{ab}d\theta^{a}d\theta^{b}$	(101)
	$\displaystyle+$	$\displaystyle 2r\sum_{\ell,m}K_{r}(t,r)E_{a}{}^{b}\partial_{b}Y_{\ell m}(\theta^{c})drd\theta^{a}+r^{2}\sum_{\ell,m}K_{3}(t,r)E_{(a}{}^{c}{\tilde{\nabla}}_{b)}{\tilde{\nabla}}_{c}Y_{\ell m}(\theta^{c})d\theta^{a}d\theta^{b}\Big{]}.$	(101)

The odd-parity perturbation of the Lagrange multipliers is given by

\displaystyle\lambda=0,\quad{\bar{\lambda}}=0,\quad\lambda^{r}=0,\quad\lambda^{a}=\sum_{\ell,m}\Lambda(t,r)E^{ab}\partial_{b}Y_{\ell m}(\theta^{c}).

(102)

Note that in the odd-parity sector the $\ell=0$ mode is absent. We first focus on the $\ell\geq 2$ modes, and then the dipolar modes $\ell=1$ separately. For the dipolar mode $\ell=1$ , the $H_{3}$ and $K_{3}$ perturbations are not present automatically.

IV.2 The higher-multipole modes $\ell\geq 2$

First, we focus on the modes $\ell\geq 2$ . Under the joint foliation-preserving spatial gauge transformation

\displaystyle t\to t,\qquad x^{i}\to x^{i}+\xi^{i}(t,x^{i}),

(103)

which for the odd-parity perturbations are explicitly given by

\displaystyle\xi^{r}=0,\qquad\xi^{a}=\sum_{\ell,m}\Xi(t,r)E^{ab}\partial_{b}Y_{\ell m}(\theta^{c}),

(104)

the metric perturbations are transformed as

	$\displaystyle{\bar{\delta}}H_{t}=-\dot{\Xi}+N^{r}\Xi^{\prime},\qquad{\bar{\delta}}H_{r}=-r\Xi^{\prime},\qquad{\bar{\delta}}H_{3}=-2\Xi,$
	$\displaystyle{\bar{\delta}}K_{t}=-\dot{\Xi}+N_{(f)}^{r}\Xi^{\prime},\qquad{\bar{\delta}}K_{r}=-r\Xi^{\prime},\qquad{\bar{\delta}}K_{3}=-2\Xi,$		(105)

where in this section dots and primes represent the derivatives with respect to $t$ and $r$ , respectively. Thus, the combinations of

\displaystyle H_{r}-K_{r},\qquad H_{3}-K_{3},

(106)

are gauge invariant. For the convenience, we introduce the gauge-invariant variables ${\bar{H}}_{t}$ , ${\bar{H}}_{r}$ ${\bar{K}}_{t}$ , ${\bar{K}}_{r}$ , and ${\bar{K}}_{3}$ by

$\displaystyle H_{t}$	$\displaystyle=$	$\displaystyle{\bar{H}}_{t}+\frac{1}{2}\left(\dot{H}_{3}-N^{r}H_{3}^{\prime}\right),$
$\displaystyle H_{r}$	$\displaystyle=$	$\displaystyle{\bar{H}}_{r}+\frac{1}{2}rH_{3}^{\prime},$
$\displaystyle K_{t}$	$\displaystyle=$	$\displaystyle{\bar{K}}_{t}+\frac{1}{2}\left(\dot{K}_{3}-N_{(f)}^{r}K_{3}^{\prime}\right),$
$\displaystyle K_{r}$	$\displaystyle=$	$\displaystyle{\bar{K}}_{r}+\frac{1}{2}rK_{3}^{\prime},$
$\displaystyle K_{3}$	$\displaystyle=$	$\displaystyle{\bar{K}}_{3}+H_{3}.$	(107)

The perturbation of the Lagrange multiplier is gauge invariant by itself.

Because of the degeneracy among the different $m$ modes in the spherically symmetric backgrounds, we may use the Legendre polynomials $P_{\ell}(\cos\theta)$ instead of the spherical harmonics $Y_{\ell m}(\theta^{c})$ . Expanding the MTBG action (4) on the Schwarzschild $-$ de Sitter backgrounds written in the spatially flat coordinates with the algebraic conditions (29) and (99) up to the second order of the perturbations and integrating over the angular directions $\theta^{a}$ , we obtain the second-order action of the odd-parity perturbations besides the $\ell=1$ mode,

\displaystyle\delta_{(2)}S

\displaystyle=

\displaystyle\sum_{\ell\neq 1}\frac{\ell(\ell+1)}{2\ell+1}\int dtdrL_{\ell}\left[H_{t},H_{r},H_{3},K_{t},K_{r},K_{3},\Lambda\right].

(108)

Note that at the level of the linearized perturbations there is no mixing between different $\ell$ modes, and so we minimize the second-order action for each $\ell$ mode. With use of Eq. (IV.2), the second-order action (108) can be rewritten in terms of the gauge-invariant variables as

\displaystyle\delta_{(2)}S

\displaystyle=

\displaystyle\sum_{\ell\neq 1}\frac{\ell(\ell+1)}{2\ell+1}\int dtdrL_{\ell}\left[\bar{H}_{t},\bar{H}_{r},\bar{K}_{t},\bar{K}_{r},\bar{K}_{3},H_{3},\Lambda\right].

(109)

Varying Eq. (109) under the algebraic conditions (29) and (99), we obtain the equations of motion for the gauge-invariant quantities $\bar{H}_{t}$ , $\bar{H}_{r}$ , $\bar{K}_{t}$ , $\bar{K}_{r}$ , $\bar{K}_{3}$ , and $\Lambda$ , as Eqs. (166), (167), (168), (169), (170), and (171), respectively. We note that Eqs. (166)-(171) are not singular at the event horizons $r=r_{(g)}$ and $r=r_{(f)}$ , and the limit to the de Sitter metrics $r_{(g)}\to 0$ and $r_{(f)}\to 0$ can be taken smoothly. The remaining mode $H_{3}$ corresponds to the gauge mode and hence the equation of motion for $H_{3}$ becomes trivial. Thus, without any loss of generality, we may set $H_{3}=0$ , and in Eq. (109) $L_{\ell}\left[\bar{H}_{t},\bar{H}_{r},\bar{K}_{t},\bar{K}_{r},\bar{K}_{3},H_{3},\Lambda\right]\to L_{\ell}\left[\bar{H}_{t},\bar{H}_{r},\bar{K}_{t},\bar{K}_{r},\bar{K}_{3},\Lambda\right]$ .

In order to derive the master equations for the odd-parity perturbations, we introduce the new variables $\chi_{h}$ and $\chi_{k}$ and define the new second-order action [39]⁵⁵5We follow here the same method used in the context of the $f(R,{\cal G})$ theories ( ${\cal G}$ being the Gauss-Bonnet scalar). Even though the method in this paper is equivalent to the one in [39], the results will differ. In particular, in the $f(R,{\cal G})$ theories, the odd-mode master variable was sufficient in order to entirely describe the dynamics of all the odd-mode fields. In MTBG, we will see that this is not the case. In fact, we will see that the shadowy modes cannot be, in general, integrated out in terms of the introduced master variables and that the same shadowy modes will act as sources for the dynamics of the same master variables. This feature might be shared also by other theories that introduce shadowy modes, a phenomenon that could be worthy of further investigation. by

$\displaystyle\delta_{(2)}S^{\prime}$	$\displaystyle=$	$\displaystyle\delta_{(2)}S-\sum_{\ell\neq 1}\frac{\ell(\ell+1)}{2\ell+1}\int dtdr\Big{\{}-A_{1}(r)\left(\dot{\bar{H}}_{r}+A_{2}(r)\bar{H}_{r}^{\prime}+A_{3}(r)\bar{H}_{t}^{\prime}+A_{4}(r)\bar{H}_{r}-\frac{\chi_{h}}{r}\right)^{2}$	(110)
		$\displaystyle+B_{1}(r)\left(\dot{\bar{K}}_{r}+B_{2}(r)\bar{K}_{r}^{\prime}+B_{3}(r)\bar{K}_{t}^{\prime}+B_{4}(r)\bar{K}_{r}-\frac{\chi_{k}}{r}\right)^{2}\Big{\}}$
	$\displaystyle=$	$\displaystyle\sum_{\ell\neq 1}\frac{\ell(\ell+1)}{2\ell+1}\int dtdr\Big{\{}L_{\ell}\left[\bar{H}_{t},\bar{H}_{r},\bar{K}_{t},\bar{K}_{r},\bar{K}_{3},\Lambda\right]-A_{1}(r)\left(\dot{\bar{H}}_{r}+A_{2}(r)\bar{H}_{r}^{\prime}+A_{3}(r)\bar{H}_{t}^{\prime}+A_{4}(r)\bar{H}_{r}-\frac{\chi_{h}}{r}\right)^{2}$
		$\displaystyle-B_{1}(r)\left(\dot{\bar{K}}_{r}+B_{2}(r)\bar{K}_{r}^{\prime}+B_{3}(r)\bar{K}_{t}^{\prime}+B_{4}(r)\bar{K}_{r}-\frac{\chi_{k}}{r}\right)^{2}\Big{\}}.$

Varying the new second-order action $\delta_{(2)}S^{\prime}$ , Eq. (110), with respect to $\chi_{h}$ and $\chi_{k}$ , we obtain

	$\displaystyle\dot{\bar{H}}_{r}+A_{2}(r)\bar{H}_{r}^{\prime}+A_{3}(r)\bar{H}_{t}^{\prime}+A_{4}(r)\bar{H}_{r}-\frac{\chi_{h}}{r}=0,$		(111)
	$\displaystyle\dot{\bar{K}}_{r}+B_{2}(r)\bar{K}_{r}^{\prime}+B_{3}(r)\bar{K}_{t}^{\prime}+B_{4}(r)\bar{K}_{r}-\frac{\chi_{k}}{r}=0.$		(112)

Solving Eqs. (111) and (112) in terms of $\chi_{h}$ and $\chi_{k}$ and then substituting them back into Eq. (110), the new second-order action (110) reduces to the original one (109).

Varying the new second-order action $\delta_{(2)}S^{\prime}$ , Eq. (110), with respect to $\bar{H}_{t}$ , $\bar{H}_{r}$ , $\bar{K}_{t}$ , $\bar{K}_{r}$ , we obtain the equations of motion for them. Requiring that the derivatives of $\bar{H}_{t}$ , $\bar{H}_{r}$ , $\bar{K}_{t}$ , $\bar{K}_{r}$ vanish in their equations, so that $\bar{H}_{t}$ , $\bar{H}_{r}$ , $\bar{K}_{t}$ , $\bar{K}_{r}$ become auxiliary variables of the second-order action (110), we determine the coefficients in Eq. (110),

	$\displaystyle A_{1}(r)=r^{2},\quad A_{2}(r)=-N^{r},\quad A_{3}(r)=-r,\quad A_{4}(r)=\frac{3}{2r}N^{r}-\frac{\Lambda_{(g)}r}{2N^{r}},$
	$\displaystyle B_{1}(r)=\frac{{\tilde{\alpha}}^{2}C_{0}^{2}}{b}r^{2},\quad B_{2}(r)=-N_{(f)}^{r},\quad B_{3}(r)=-r,\quad B_{4}(r)=\frac{3}{2r}N_{(f)}^{r}-\frac{\Lambda_{(f)}C_{0}^{2}b^{2}r}{2N_{(f)}^{r}}.$		(113)

The equations of motion for $\bar{H}_{t}$ , $\bar{H}_{r}$ , $\bar{K}_{t}$ , and $\bar{K}_{r}$ then relate them to the variables $\chi_{h}$ and $\chi_{k}$ as

	$\displaystyle\bar{H}_{t}$	$\displaystyle=-\frac{r\chi_{h}^{\prime}+2\chi_{h}}{(\ell+2)(\ell-1)r},$
	$\displaystyle\bar{H}_{r}$	$\displaystyle=\frac{1}{2\left[(1-b)C_{0}^{2}(C_{0}c_{1}+c_{2})m^{2}r^{2}-2(\ell+2)(\ell-1)\right]}$
		$\displaystyle\times\left[C_{0}\left(C_{0}c_{1}+c_{2}\right)m^{2}r^{2}\left\{(1-b)C_{0}\left(r\bar{K}_{3}^{\prime}+2\bar{K}_{r}\right)+2r(1-C_{0}\beta)\Lambda^{\prime}\right\}-4\left(N^{r}\left(2\chi_{h}+r\chi_{h}^{\prime}\right)-r\dot{\chi}_{h}\right)\right],$
	$\displaystyle\bar{K}_{t}$	$\displaystyle=-\frac{r\chi_{k}^{\prime}+2\chi_{k}}{(\ell+2)(\ell-1)r},$
	$\displaystyle\bar{K}_{r}$	$\displaystyle=\frac{1}{2bC_{0}\left[(1-b)(C_{0}c_{1}+c_{2})m^{2}r^{2}-2(\ell+2)(\ell-1)b{\tilde{\alpha}}^{2}\right]}$
		$\displaystyle\times\left[b\left(C_{0}c_{1}+c_{2}\right)m^{2}r^{2}\left\{\left((1-b)C_{0}\left(-r\bar{K}_{3}^{\prime}+2\bar{H}_{r}\right)-2r(1-C_{0}\beta)\Lambda^{\prime}\right)\right\}-4C_{0}{\tilde{\alpha}}^{2}\left(N_{(f)}^{r}\left(2\chi_{k}+r\chi_{k}^{\prime}\right)-r\dot{\chi}_{k}\right)\right].$

Since in the above expression (IV.2) $\bar{H}_{r}$ and $\bar{K}_{r}$ are coupled, we still need to solve for each of them, although we omit to show them explicitly. $\chi_{h}$ and $\chi_{k}$ can be interpreted as the master variables for the odd-parity perturbations for $\ell\geq 2$ . Substituting Eq. (IV.2) into Eqs. (111) and (112), we obtain

	$\displaystyle\chi_{h}-\left(\frac{3}{2}N^{r}-\frac{r^{2}\Lambda_{(g)}}{2N^{r}}\right)\bar{H}_{r}+r\left(N^{r}\bar{H}_{r}^{\prime}+r\bar{H}_{t}^{\prime}-\dot{\bar{H}}_{r}\right)=0,$		(115)
	$\displaystyle\chi_{k}-\left(\frac{3}{2}N_{(f)}^{r}-\frac{C_{0}^{2}b^{2}r^{2}\Lambda_{(f)}}{2N_{(f)}^{r}}\right)\bar{K}_{r}+r\left(N_{(f)}^{r}\bar{K}_{r}^{\prime}+r\bar{K}_{t}^{\prime}-\dot{\bar{K}}_{r}\right)=0,$		(116)

and moreover substituting Eq. (IV.2) into Eqs. (115) and (116), we obtain the master equations for $\chi_{h}$ and $\chi_{k}$ .

Varying the action (110) with Eq. (IV.2) with respect to $\Lambda$ and $K_{3}$ leads to the equations of motion for $\Lambda$ and $\bar{K}_{3}$ , which are given by

	$\displaystyle r^{2}\Lambda^{\prime\prime}+4r\Lambda^{\prime}-\left(\ell+2\right)\left(\ell-1\right)\Lambda=0,$		(117)
	$\displaystyle r^{2}\bar{K}_{3}^{\prime\prime}+4r\bar{K}_{3}^{\prime}-\left(\ell+2\right)\left(\ell-1\right)\bar{K}_{3}=-2r\left[\left(\bar{K}_{r}^{\prime}-\bar{H}_{r}^{\prime}\right)+\frac{3}{r}\left(\bar{K}_{r}-\bar{H}_{r}\right)\right],$		(118)

which are also obtained by combining Eqs. (170) and (171). Eliminating $\bar{H}_{r}$ and $\bar{K}_{r}$ in Eq. (118) with Eq. (IV.2), the equation (118) reduces to the form

$\displaystyle r^{2}\bar{K}_{3}^{\prime\prime}$	$\displaystyle+\left\{4r-\frac{2r^{3}\left({\tilde{\alpha}}^{2}C_{0}^{2}b+1\right)m^{2}\left(b-1\right)\left(C_{0}c_{1}+c_{2}\right)}{(C_{0}^{2}{\tilde{\alpha}}^{2}b+1)\,r^{2}(C_{0}c_{1}+c_{2})(b-1)m^{2}+2b(\ell^{2}+\ell-2){\tilde{\alpha}}^{2}}\right\}\bar{K}_{3}^{\prime}$
	$\displaystyle-\left\{\left(\ell+2\right)\left(\ell-1\right)+\frac{m^{2}r^{2}\left(b-1\right)\left({\tilde{\alpha}}^{2}bC_{0}^{2}+1\right)\left(C_{0}c_{1}+c_{2}\right)}{2{\tilde{\alpha}}^{2}b}\right\}\bar{K}_{3}$
	$\displaystyle={\cal F}\left[\dot{\chi}_{h}^{\prime},\chi_{h}^{\prime\prime},\dot{\chi}_{h},\chi_{h}^{\prime},\chi_{h};\dot{\chi}_{k}^{\prime},\chi_{k}^{\prime\prime},\dot{\chi}_{k},\chi_{k}^{\prime},\chi_{k};\Lambda^{\prime\prime},\Lambda^{\prime},\Lambda\right],$	(119)

where ${\cal F}$ represents a linear combination of the variables shown in the arguments, which is not explicitly written. Since Eqs. (117) and (IV.2) are elliptic equations, $\Lambda$ and $\bar{K}_{3}$ are fixed under the given spatial boundary conditions. Thus, $\Lambda$ and $\bar{K}_{3}$ correspond to the shadowy modes [27, 29]. In Eq. (117), $\Lambda$ is not coupled to the other variables and can be solved as

\displaystyle\Lambda=g_{1}(t)r^{-2-\ell}+g_{2}(t)r^{-1+\ell},

(120)

where $g_{1}(t)$ and $g_{2}(t)$ are integration constants. For instance, requiring the regularity at the spatial infinity $r\to\infty$ for the $\ell\geq 2$ modes, we may impose $g_{2}(t)=0$ . In the rest, however, we leave $g_{1}(t)$ and $g_{2}(t)$ as the general functions of the time $t$ .

IV.2.1 The case of $b=1$

In the case of $b=1$ , the substitution of Eq. (IV.2) into Eqs. (115) and (116) yields the master equations

	$\displaystyle\ddot{\chi}_{h}-2N^{r}\dot{\chi}_{h}^{\prime}+[\left(N^{r}\right)^{2}-1]\chi_{h}^{\prime\prime}-\frac{\Lambda_{(g)}r^{2}+3\left(N^{r}\right)^{2}}{2r\,N^{r}}\dot{\chi}_{h}$
	$\displaystyle+\frac{\Lambda_{(g)}r^{2}+\left(N^{r}\right)^{2}-2}{r}\chi_{h}^{\prime}+\frac{2\Lambda_{(g)}r^{2}-4\left(N^{r}\right)^{2}+\ell^{2}+\ell}{r^{2}}\chi_{h}$
	$\displaystyle=\frac{\left(C_{0}c_{1}+c_{2}\right)m^{2}r\left(1-\beta C_{0}\right)C_{0}\{[\Lambda_{(g)}r^{2}+3\left(N^{r}\right)^{2}]\Lambda^{\prime}+2r\,N^{r}(\Lambda^{\prime\prime}N^{r}-\dot{\Lambda}^{\prime})\}}{4N^{r}},$		(121)
	$\displaystyle\ddot{\chi}_{k}-2N_{(f)}^{r}\dot{\chi}_{k}^{\prime}+[(N_{(f)}^{r})^{2}-1]\chi_{k}^{\prime\prime}-\frac{C_{0}^{2}\Lambda_{(f)}r^{2}+3(N_{(f)}^{r})^{2}}{2rN_{(f)}^{r}}\dot{\chi}_{k}$
	$\displaystyle+\frac{C_{0}^{2}\Lambda_{(f)}r^{2}+(N_{(f)}^{r})^{2}-2}{r}\chi_{k}^{\prime}+\frac{2C_{0}^{2}\Lambda_{(f)}r^{2}-4(N_{(f)}^{r})^{2}+\ell^{2}+\ell}{r^{2}}\chi_{k}$
	$\displaystyle=-\frac{m^{2}\left(C_{0}c_{1}+c_{2}\right)(1-\beta C_{0})r\left(\Lambda^{\prime}C_{0}^{2}\Lambda_{(f)}r^{2}+2(N_{(f)}^{r})^{2}r\Lambda^{\prime\prime}+3\Lambda^{\prime}(N_{(f)}^{r})^{2}-2N_{(f)}^{r}r\dot{\Lambda}^{\prime}\right)}{4N_{(f)}^{r}C_{0}{\tilde{\alpha}}^{2}}.$		(122)

The left-hand sides of Eqs. (IV.2.1) and (IV.2.1) correspond to the operators for the master equations in the two copies of GR. In our case, Eqs. (IV.2.1) and (IV.2.1), together with Eqs. (117) and (IV.2), form a closed system. In Eqs. (IV.2.1) and (IV.2.1), the two dynamical modes $\chi_{h}$ and $\chi_{k}$ are decoupled and sourced by the shadowy mode $\Lambda$ . From Eqs. (IV.2.1) and (IV.2.1), we find that in both the sectors the odd-parity perturbations propagate with the speed of light.

With the explicit form of the solution for $\Lambda$ , Eq. (120), the equations for the odd-parity perturbations reduce to

	$\displaystyle\ddot{\chi}_{h}-2N^{r}\dot{\chi}_{h}^{\prime}+[\left(N^{r}\right)^{2}-1]\chi_{h}^{\prime\prime}-\frac{\Lambda_{(g)}r^{2}+3\left(N^{r}\right)^{2}}{2r\,N^{r}}\dot{\chi}_{h}+\frac{\Lambda_{(g)}r^{2}+\left(N^{r}\right)^{2}-2}{r}\chi_{h}^{\prime}+\frac{2\Lambda_{(g)}r^{2}-4\left(N^{r}\right)^{2}+\ell^{2}+\ell}{r^{2}}\chi_{h}$
	$\displaystyle=\frac{m^{2}\left(C_{0}c_{1}+c_{2}\right)(1-\beta C_{0})C_{0}r}{2N^{r}}\left\{g_{1}\left(\ell+\frac{3}{2}\right)\left(N^{r}\right)^{2}\left(\ell+2\right)r^{-\ell-3}+\dot{g}_{1}N^{r}\left(\ell+2\right)r^{-\ell-2}-\frac{\Lambda_{(g)}g_{1}\left(\ell+2\right)r^{-\ell-1}}{2}\right.$
	$\displaystyle+\left.\left(\ell-1\right)\left(g_{2}\left(N^{r}\right)^{2}\left(\ell-\frac{1}{2}\right)r^{\ell-2}+\frac{\Lambda_{(g)}g_{2}r^{\ell}}{2}-r^{\ell-1}\dot{g}_{2}N^{r}\right)\right\},$		(123)
	$\displaystyle\ddot{\chi}_{k}-2N_{(f)}^{r}\dot{\chi}_{k}^{\prime}+[(N_{(f)}^{r})^{2}-1]\chi_{k}^{\prime\prime}-\frac{C_{0}^{2}\Lambda_{(f)}r^{2}+3(N_{(f)}^{r})^{2}}{2rN_{(f)}^{r}}\dot{\chi}_{k}$
	$\displaystyle+\frac{C_{0}^{2}\Lambda_{(f)}r^{2}+(N_{(f)}^{r})^{2}-2}{r}\chi_{k}^{\prime}+\frac{2C_{0}^{2}\Lambda_{(f)}r^{2}-4(N_{(f)}^{r})^{2}+\ell^{2}+\ell}{r^{2}}\chi_{k}$
	$\displaystyle=-\frac{m^{2}\left(C_{0}c_{1}+c_{2}\right)\left(1-\beta C_{0}\right)r}{2N_{(f)}^{r}C_{0}{\tilde{\alpha}}^{2}}\left\{(N_{(f)}^{r})^{2}\left(\ell+\frac{3}{2}\right)\left(\ell+2\right)g_{1}r^{-\ell-3}+\dot{g}_{1}N_{(f)}^{r}\left(\ell+2\right)r^{-\ell-2}-\frac{C_{0}^{2}\Lambda_{(f)}g_{1}\left(\ell+2\right)r^{-\ell-1}}{2}\right.$
	$\displaystyle+\left.\left(\ell-1\right)\left(g_{2}(N_{(f)}^{r})^{2}\left(\ell-\frac{1}{2}\right)r^{\ell-2}+\frac{C_{0}^{2}\Lambda_{(f)}g_{2}r^{\ell}}{2}-r^{\ell-1}\dot{g}_{2}N_{(f)}^{r}\right)\right\}.$		(124)

IV.2.2 The case of $b\neq 1$

In the case of $b\neq 1$ , after eliminating $\bar{H}_{t}$ , $\bar{H}_{r}$ , $\bar{K}_{t}$ , and $\bar{K}_{r}$ in Eqs. (115) and (116) with Eq. (IV.2) and combining them appropriately, we obtain the master equations where the kinetic terms are diagonalized. Although the master equations are too long to be shown explicitly, their general structure is given by

	$\displaystyle\ddot{\chi}_{h}+M_{1}(r)\chi_{h}+M_{2}(r)\chi_{k}+L_{1}(r)\dot{\chi}_{h}^{\prime}+L_{2}(r)\chi_{h}^{\prime\prime}+L_{3}(r)\dot{\chi}_{h}+L_{4}(r)\chi_{h}^{\prime}+P_{1}(r)\dot{\chi}_{k}^{\prime}+P_{2}(r)\chi_{k}^{\prime\prime}+P_{3}(r)\dot{\chi}_{k}+P_{4}(r)\chi_{k}^{\prime}$
	$\displaystyle=Q_{1}(r)\dot{\bar{K}}_{3}^{\prime}+Q_{2}(r)\bar{K}_{3}^{\prime\prime}+Q_{3}(r)\bar{K}_{3}^{\prime}+S_{1}(r)\dot{\Lambda}^{\prime}+S_{2}(r)\Lambda^{\prime\prime}+S_{3}(r)\Lambda^{\prime},$		(125)
	$\displaystyle\ddot{\chi}_{k}+{\tilde{M}}_{1}(r)\chi_{k}+{\tilde{M}}_{2}(r)\chi_{h}+\tilde{L}_{1}(r)\dot{\chi}_{k}^{\prime}+\tilde{L}_{2}(r)\chi_{k}^{\prime\prime}+\tilde{L}_{3}(r)\dot{\chi}_{k}+\tilde{L}_{4}(r)\chi_{k}^{\prime}+{\tilde{P}}_{1}(r)\dot{\chi}_{h}^{\prime}+{\tilde{P}}_{2}(r)\chi_{h}^{\prime\prime}+{\tilde{P}}_{3}(r)\dot{\chi}_{h}+{\tilde{P}}_{4}(r)\chi_{h}^{\prime}$
	$\displaystyle={\tilde{Q}}_{1}(r)\dot{\bar{K}}_{3}^{\prime}+{\tilde{Q}}_{2}(r)\bar{K}_{3}^{\prime\prime}+{\tilde{Q}}_{3}(r)\bar{K}_{3}^{\prime}+{\tilde{S}}_{1}(r)\dot{\Lambda}^{\prime}+{\tilde{S}}_{2}(r)\Lambda^{\prime\prime}+{\tilde{S}}_{3}(r)\Lambda^{\prime},$		(126)

where the background-dependent coefficients $M_{1,2}(r)$ , $L_{1,2,3,4}(r)$ , $P_{1,2,3,4}(r)$ , $Q_{1,2,3}(r)$ , $S_{1,2,3}(r)$ , ${\tilde{M}}_{1,2}(r)$ , ${\tilde{L}}_{1,2,3,4}(r)$ , ${\tilde{P}}_{1,2,3,4}(r)$ , ${\tilde{Q}}_{1,2,3}(r)$ , $\tilde{S}_{1,2,3}(r)$ are the pure functions of $r$ . In the limit of $b\to 1$ , $M_{2}(r)$ , $P_{1,2,3,4}(r)$ , $Q_{1,2,3}(r)$ , ${\tilde{M}}_{2}(r)$ , ${\tilde{P}}_{1,2,3,4}(r)$ , and ${\tilde{Q}}_{1,2,3}(r)$ vanish, and Equations. (IV.2.2) and (IV.2.2) reduce to Eqs. (IV.2.1) and (IV.2.1), respectively.

On the other hand, in the general case of $b\neq 1$ the above coefficients do not vanish, and the two dynamical modes $\chi_{h}$ and $\chi_{k}$ are coupled to each other, respectively. Eqs. (IV.2.2) and (IV.2.2), together with Eqs. (117) and (IV.2), form the closed system. While the two modes $\chi_{h}$ and $\chi_{k}$ are dynamical, the remaining ones $\bar{K}_{3}$ and $\Lambda$ are shadowy. From Eqs.(IV.2.1) and (IV.2.1), the dynamical $\chi_{h}$ and $\chi_{k}$ are coupled and sourced by the shadowy modes $\Lambda$ and $\bar{K}_{3}$ . From Eqs. (IV.2.2) and (IV.2.2), we define the mass matrix by

\displaystyle{\cal M}=\begin{pmatrix}M_{1}(r)&M_{2}(r)\\ {\tilde{M}}_{2}(r)&{\tilde{M}}_{1}(r)\end{pmatrix}.

(127)

In the large distance limit $r\gg{\rm max}(r_{(g)},r_{(f)})$ and assuming astrophysical scales, i.e., neglecting the cosmological constants terms, we find

	$\displaystyle M_{1}(r)$	$\displaystyle=$	$\displaystyle M_{1,0}+{\cal O}\left(\frac{1}{r}\right),\qquad M_{2}(r)=-M_{1,0}+{\cal O}\left(\frac{1}{r}\right),$
	$\displaystyle{\tilde{M}}_{1}(r)$	$\displaystyle=$	$\displaystyle{\tilde{M}}_{1,0}+{\cal O}\left(\frac{1}{r}\right),\qquad{\tilde{M}}_{2}(r)=-{\tilde{M}}_{1,0}+{\cal O}\left(\frac{1}{r}\right),$		(128)

where

\displaystyle M_{1,0}:=\frac{m^{2}C_{0}^{2}(C_{0}c_{1}+c_{2})(b-1)\ell(\ell+1)}{2(\ell+2)(\ell-1)},\qquad{\tilde{M}}_{1,0}:=\frac{m^{2}(C_{0}c_{1}+c_{2})b(b-1)\ell(\ell+1)}{2(\ell+2)(\ell-1){\tilde{\alpha}}^{2}}.

(129)

This squared mass is of order of $m^{2}\sim H_{0}^{2}$ for cosmological implementation of this MTBG. Only for astrophysical purposes we can neglect its contribution to the dynamics of the fields.

So far, we have considered the parameter $b$ as a free parameter, as it is not fixed by any background equation of motion for the spatially flat coordinates. However, as we have seen in Sec. III, the collapse fixes the value of $b$ at least for the cases we have found a solution. Indeed, in the spatially closed case, $b$ was set to unity, whereas for the spatially flat case, $b$ was determined by the amount of matter energy densities in the physical and fiducial sectors. However, we do not know all the collapse solutions, and therefore there could be other possibilities for both the interior and exterior metrics. In particular, as we have already discussed in the spatially closed collapsing case, a time dependence on $b$ might appear that could determine a deviation from a stationary configuration for the exterior metric too.

Furthermore, we might wonder what boundary conditions we should impose on the shadowy modes. We believe that these topics have a nontrivial answer, mostly because it may depend on the system we consider. In fact, we know that the theory possesses shadowy modes, and as such, we should expect that the behavior of the solution to depend on the boundary conditions which we fix on them. These boundary conditions might differ depending on the environment of the solution itself. For instance, if we were to consider a typical astrophysical system, before reaching cosmological scales, we would encounter the inhomogeneities induced by the matter distribution, such as the baryonic matter contribution coming from the galaxy and then the dark matter halo distribution. Therefore, in this case we would need to match the astrophysical scale solution with an environment filled with matter. On top of that, we would need to fix boundary conditions for the shadowy modes at these same scales.

Instead if we were to describe an approximately spherically symmetric large cluster of galaxies, then it would make sense to choose cosmological boundary conditions to impose not only for the shadowy modes but also for the value of $b$ . Because of these possibilities, we will leave the parameter $b$ as a free parameter.

IV.3 Conditions for the absence of ghost and gradient instabilities

In this subsection, we are going to study the conditions for the absence of the ghost and gradient instabilities for the odd modes. After eliminating $\bar{H}_{t}$ , $\bar{H}_{r}$ , $\bar{K}_{t}$ , and $\bar{K}_{r}$ with Eq. (IV.2), the second-order action (110) can be rewritten as the functional of the two dynamical modes $\chi_{h}$ , $\chi_{k}$ and two shadowy modes $\bar{K}_{3}$ , $\Lambda$ . The shadowy modes cannot be integrated out, in general, and they source the dynamical fields $\chi_{h}$ and $\chi_{k}$ . So, in the following, we will also assume that the shadowy modes $\Lambda$ and $\bar{K}_{3}$ will not be given boundary conditions as to drastically change the dynamics of the equations of motion for $\chi_{h}$ and $\chi_{k}$ . On neglecting a possible strong backreaction of the shadowy modes on $\chi_{h,k}$ assuming appropriate boundary conditions are imposed, we can find the kinetic matrix for $\chi_{h,k}$ .

The no-ghost conditions can be obtained by imposing the positivity of the eigenvalues of the kinetic matrix, ⁶⁶6To study the sign of the eigenvalues, it is sufficient to diagonalize the symmetric kinetic matrix via a simple congruence diagonalization method using a unitriangular matrix (see e.g. Ref. [40]), which corresponds to make a field redefinition, with determinant equal to unity, for the perturbation variables. Then, the positivity conditions on the sign of eigenvalues are equivalent to imposing the diagonal elements to be positive. which is constituted by the coefficients of the $\dot{\chi}_{h}^{2}$ , $\dot{\chi}_{k}^{2}$ , and $\dot{\chi}_{h}\dot{\chi}_{k}$ terms in the reduced second-order action. In the following, we will also consider the case of the large $\omega$ - $k_{r}$ - $\ell$ case, so as to assume that time and $r$ derivatives give large contributions as well as $\ell\gg 1$ , e.g., $r|\chi^{\prime}_{h}|\gg|\chi_{h}|$ , etc. When evaluated in the high- $\ell$ regime, the off-diagonal components of the kinetic matrix are suppressed by $1/\ell^{2}$ , and the no-ghost conditions in this regime are given by

\frac{C_{0}^{2}\tilde{\alpha}^{2}r^{2}}{4b^{3}}>0\,,\qquad\frac{r^{2}}{4}>0\,,

(130)

which are both trivially satisfied.

As for the speed of propagation, we can consider the leading contribution in the equations of motion for $\chi_{h}$ and $\chi_{k}$ , once more, on neglecting backreactions from the shadowy modes and assuming $m\simeq H_{0}$ . In this case the equations of motion simplify to

	$\displaystyle\ddot{\chi}_{h}$	$\displaystyle-2N^{r}\dot{\chi}^{\prime}_{h}+[(N^{r})^{2}-1]\,\chi^{\prime\prime}_{h}+\frac{\ell(\ell+1)}{r^{2}}\,\chi_{h}\simeq 0\,,$		(131)
	$\displaystyle\ddot{\chi}_{k}$	$\displaystyle-2N_{(f)}^{r}\dot{\chi}^{\prime}_{k}+[(N_{(f)}^{r})^{2}-b^{2}]\,\chi^{\prime\prime}_{k}+\frac{b^{2}\ell(\ell+1)}{r^{2}}\,\chi_{k}\simeq 0\,.$		(132)

The first case is simple, as it reduces to the leading contributions for a massless scalar propagating in the physical metric and the speed of propagation is evidently $c_{s,g}^{2}=1$ . As for the second equation, in order to make the analysis simpler, one can perform a coordinate transformation, namely,

\displaystyle dt=d\tau,\qquad dr=N_{(f)}^{r}\,(d\rho-d\tau)\,,

(133)

for which $r=r[\rho-\tau]$ . This leads to a metric in the so-called Lemaître coordinates as in

ds_{f}^{2}=C_{0}^{2}\,[-b^{2}d\tau^{2}+(N_{(f)}^{r})^{2}\,d\rho^{2}+r[\rho-\tau]^{2}\,\theta_{ab}d\theta^{a}d\theta^{b}]\,.

(134)

After little algebra, one finds

\displaystyle\frac{\partial}{\partial t}=\frac{\partial}{\partial\tau}+\frac{\partial}{\partial\rho}\,,\qquad\frac{\partial}{\partial r}=\frac{1}{N_{(f)}^{r}}\,\frac{\partial}{\partial\rho}\,,

(135)

where we have used the conditions holding for a coordinate basis, $dt(\partial_{t})=1$ , $dt(\partial_{r})=0$ , etc. In this case Eq. (132) reduces to

\frac{\partial^{2}\chi_{k}}{\partial\tau^{2}}-\frac{b^{2}}{(N_{(f)}^{r})^{2}}\,\frac{\partial^{2}\chi_{k}}{\partial\rho^{2}}+\frac{b^{2}\ell(\ell+1)}{r^{2}}\,\chi_{k}\simeq 0\,.

(136)

If we were assuming the mode traveling on the spacetime of the metric $ds_{f}^{2}=-d\tau^{2}+(N_{(f)}^{r})^{2}\,d\rho^{2}+r[\rho-\tau]^{2}\,\theta_{ab}d\theta^{a}d\theta^{b}$ , then the speed of both the radial and angular propagations is $c_{f}^{2}=b^{2}$ , the proper radial distance being $(N_{(f)}^{r})\,d\rho$ [41], which does not lead to any new condition to impose.

IV.4 The dipolar modes

IV.4.1 The general solution

For the dipolar mode $\ell=1$ , $K_{3}$ and $H_{3}$ automatically vanish. Then, the second-order action is a functional of $H_{t},H_{r},K_{t}$ and $K_{r}$ together with $\Lambda$ . However, the method employed in Sec. IV.2 can be implemented.

The second-order action for the $\ell=1$ mode is given by

\displaystyle\delta_{(2)}S_{\ell=1}

\displaystyle=

\displaystyle\frac{2}{3}\int dtdrL_{1}\left[H_{t},H_{r},K_{t},K_{r},\Lambda\right].

(137)

As in Sec. IV.2, we introduce the new variables $\chi_{h}$ and $\chi_{k}$ ,

	$\displaystyle\delta_{(2)}S^{\prime}_{\ell=1}$	$\displaystyle=$	$\displaystyle\frac{2}{3}\int dtdr\Big{\{}L_{1}\left[H_{t},H_{r},K_{t},K_{r},\Lambda\right]-A_{1}(r)\left(\dot{H}_{r}+A_{2}(r)H_{r}^{\prime}+A_{3}(r)H_{t}^{\prime}+A_{4}(r)H_{r}-\frac{\chi_{h}}{r}\right)^{2}$		(138)
			$\displaystyle-B_{1}(r)\left(\dot{K}_{r}+B_{2}(r)K_{r}^{\prime}+B_{3}(r)K_{t}^{\prime}+B_{4}(r)K_{r}-\frac{\chi_{k}}{r}\right)^{2}\Big{\}}.$		(138)

Varying the new second-order action $\delta_{(2)}S_{\ell=1}^{\prime}$ with respect to $\chi_{h}$ and $\chi_{k}$ , we obtain

	$\displaystyle\dot{H}_{r}+A_{2}(r)H_{r}^{\prime}+A_{3}(r)H_{t}^{\prime}+A_{4}(r)H_{r}-\frac{\chi_{h}}{r}=0,$		(139)
	$\displaystyle\dot{K}_{r}+B_{2}(r)K_{r}^{\prime}+B_{3}(r)K_{t}^{\prime}+B_{4}(r)K_{r}-\frac{\chi_{k}}{r}=0.$		(140)

Varying $\delta_{(2)}S_{\ell=1}^{\prime}$ with respect to $H_{t}$ , $H_{r}$ , $K_{t}$ , $K_{r}$ , we obtain the equations of motion for them. Requiring that the derivatives of $H_{t}$ , $H_{r}$ , $K_{t}$ , $K_{r}$ vanish in their equations, the coefficients in Eq. (138) are given by Eq. (IV.2).

In the case of $b\neq 1$ as well as $C_{0}c_{1}+c_{2}\neq 0$ , we can integrate out $H_{r}$ using its own algebraic equation of motion, as in

\displaystyle H_{r}=\frac{r\left(\beta C_{0}-1\right)\Lambda^{\prime}}{C_{0}\left(b-1\right)}+K_{r}+\frac{\left(2N^{r}\chi^{\prime}_{h}-2\dot{\chi}_{h}\right)r+4N^{r}\chi_{h}}{r^{2}C_{0}^{2}\left(b-1\right)\left(C_{0}c_{1}+c_{2}\right)m^{2}}\,.

(141)

After inserting this expression into Eq. (138), we have that $H_{t}$ , $K_{t}$ and $K_{r}$ are all Lagrange multipliers that enter only linearly as to set constraints on the other fields. In any case, we can now find all the relevant equations of motion and try to solve them. In particular, the constraint for $H_{t}$ reads

r\chi_{h}^{\prime}+2\chi_{h}=0\,,\qquad\chi_{h}=\frac{F_{h}(t)}{r^{2}}\,,

(142)

whereas the constraint for $K_{t}$ can be written and solved as

r\chi_{k}^{\prime}+2\chi_{k}=0\,,\qquad\chi_{k}=\frac{F_{k}(t)}{r^{2}}\,.

(143)

Using these solutions, the constraint for the field $K_{r}$ reads

C_{0}^{2}\tilde{\alpha}^{2}\dot{F}_{k}+b\dot{F}_{h}=0\,,\qquad{\rm where}\qquad F_{k}=F_{\mathit{k}0}-\frac{bF_{h}}{C_{0}^{2}\tilde{\alpha}^{2}}\,.

(144)

. On using these solutions in the equation of motion for $\Lambda$ , we find

r^{2}\Lambda^{\prime\prime}+4r\Lambda^{\prime}=0\,,\qquad\Lambda=g_{2}(t)+\frac{g_{1}(t)}{r^{3}}\,.

(145)

We can evaluate the following gauge-invariant combinations for the $\ell=1$ mode:

$\displaystyle G_{r}$	$\displaystyle:=H_{r}-K_{r},$	(146)
$\displaystyle G_{t}$	$\displaystyle:=H_{t}-K_{t}+\frac{N^{r}}{r}\,H_{r}-\frac{N_{(f)}^{r}}{r}\,K_{r}\,,$	(147)
$\displaystyle G_{h}$	$\displaystyle:=H^{\prime}_{t}-\frac{\dot{H}_{r}}{r}+\left(\frac{N^{r}H_{r}}{r}\right)^{\prime}\,,$	(148)
$\displaystyle G_{k}$	$\displaystyle:=K^{\prime}_{t}-\frac{\dot{K}_{r}}{r}+\left(\frac{N_{(f)}^{r}K_{r}}{r}\right)^{\prime}\,.$	(149)

From Eqs. (141), (142) and (143), we obtain

\displaystyle G_{r}=\frac{3\left(1-\beta C_{0}\right)g_{1}}{C_{0}\left(b-1\right)r^{3}}-\frac{2\dot{F}_{h}}{r^{3}C_{0}^{2}\left(b-1\right)\left(C_{0}c_{1}+c_{2}\right)m^{2}}.

(150)

Considering a linear combination of the equations of motion for $\chi_{h}$ and $\chi_{k}$ , we find

G^{\prime}_{t}=\frac{F_{\mathit{k0}}}{r^{4}}-\frac{\left(\tilde{\alpha}^{2}C_{0}^{2}+b\right)F_{h}}{\tilde{\alpha}^{2}r^{4}C_{0}^{2}}+\frac{3\left(1-\beta C_{0}\right)\dot{g}_{1}}{C_{0}r^{4}\left(b-1\right)}-\frac{2\ddot{F}_{h}}{C_{0}^{2}r^{4}m^{2}\left(C_{0}c_{1}+c_{2}\right)\left(b-1\right)}\,,

(151)

which can be easily integrated. Finally, the equation of motion for $\chi_{k}$ gives

G_{k}=\frac{bF_{h}}{r^{4}C_{0}^{2}\tilde{\alpha}^{2}}-\frac{F_{\mathit{k}0}}{r^{4}}\,.

(152)

We can also find

G_{h}=-\frac{F_{h}}{r^{4}}\,.

(153)

IV.4.2 The slow-rotation limit of the Kerr $-$ de Sitter solution

By introducing a function of time $k_{0}(t)$ by

\displaystyle G_{r}=-\frac{k_{0}(t)}{r^{3}},

(154)

with Eq. (150), $g_{1}(t)$ can be eliminated as

\displaystyle g_{1}(t)=\frac{(b-1)C_{0}^{2}(C_{0}c_{1}+c_{2})m^{2}k_{0}(t)-2\dot{F}_{h}(t)}{3C_{0}(C_{0}c_{1}+c_{2})(C_{0}\beta-1)m^{2}}.

(155)

As the solution for the physical sector, we consider the slow-rotation limit of the Kerr $-$ de Sitter solution

\displaystyle H_{t}=-\frac{1}{1-(N^{r})^{2}}\left(\omega_{0}+\frac{2J_{(g)0}}{r^{3}}\right),\qquad H_{r}=\frac{rN^{r}}{1-(N^{r})^{2}}\left(\omega_{0}+\frac{2J_{(g)0}}{r^{3}}\right),

(156)

where $\omega_{0}$ and $J_{(g)0}$ are constants, which satisfies $H_{r}+rN^{r}H_{t}=0$ . This solution describing the slow-rotation limit of the Kerr $-$ de Sitter solution has been found as follows. First, we consider the limit for small rotations in Boyer-Lindquist coordinates for the Kerr $-$ de Sitter metric, which reduces to the Schwarzschild $-$ de Sitter contribution written in the Schwarzschild coordinates plus a term of the form $ds^{2}\ni-2a(N^{r})^{2}(1-z^{2})dt_{s}d\varphi_{s}$ , where $a$ is the spin angular momentum per unit mass, $z:=\cos\theta$ , and $t_{s}$ is the time of the Schwarzschild coordinates. Since we want the metric in the GP coordinates, we make the coordinate transformation $dt_{s}=dt+F^{\prime}\,dr$ with $F^{\prime}=-N^{r}/(1-(N^{r})^{2})$ together with a time-dependent shift in the angular variable, as in $\varphi_{s}=\varphi+\omega_{0}\,t$ , describing a reference system rotating uniformly, not due to gravitational effects. These transformations lead to a metric that can be expanded in terms of both $\omega_{0}$ and $a$ . At lowest order, we find the Schwarzschild $-$ de Sitter metric written in GP coordinates, as expected. At linear level, we find two contributions that can be written as a contribution to the shift, $g_{t\varphi}=(\omega_{0}r^{2}-a\,(N^{r})^{2})(1-z^{2})$ , and an element $g_{r\varphi}=a\xi^{3}(1-z^{2})/(1-(N^{r})^{2})$ . In turn, matching these elements to the perturbation variables for the dipolar modes in the physical sector gives the modes $\tilde{H}_{t}=\frac{\left(-\omega_{0}(N^{r})^{2}+\omega_{0}\right)r^{3}+2J_{(g)0}}{\left(({N^{r}})^{2}-1\right)r^{3}}$ and $\tilde{H}_{r}=-\frac{2J_{(g)0}N^{r}}{r^{2}\left(({N^{r}})^{2}-1\right)}$ , where we have replaced the spin angular momentum per unit mass $a$ with $-2J_{(g)0}/r_{(g)}$ . On performing an $r$ -dependent gauge transformation defined by $\Xi^{\prime}(r)=\frac{N^{r}\omega_{0}}{({N^{r}})^{2}-1}$ , we find the fields⁷⁷7It should be noted that the combination $H_{t}+N^{r}H_{r}/r$ for an $r$ -dependent gauge transformation becomes a gauge-invariant quantity; as such, it is also equal to ${\tilde{H}}_{t}+N^{r}{\tilde{H}}_{r}/r=-\omega_{0}-\frac{2J}{r^{3}}$ . in the form given in Eq. (156).

The constant $J_{(g)0}$ is then linked to the angular momentum in the slow-rotation limit of a Kerr $-$ de Sitter black hole. With Eq. (154), the solution for $K_{r}$ is then

\displaystyle K_{r}=\left(\omega_{0}+\frac{2J_{(g)0}}{r^{3}}\right)\frac{rN^{r}}{1-(N^{r})^{2}}+\frac{k_{0}(t)}{r^{3}}.

(157)

The equation of motion for $\chi_{h}$ provides

\displaystyle F_{h}=-6J_{(g)0}.

(158)

The equation of motion for $\chi_{k}$ can then be solved as

\displaystyle K_{t}=-\frac{2J_{(f)0}}{r^{3}}-\frac{N^{r}N^{r}_{(f)}}{1-(N^{r})^{2}}\left(\omega_{0}+\frac{2J_{(g)0}}{r^{3}}\right)-\frac{1}{r^{3}}\left(\frac{1}{r}N_{(f)}^{r}k_{0}(t)+\frac{\dot{k}_{0}(t)}{3}\right)-k_{3}(t),

(159)

where the constant $J_{(f)0}$ is introduced by fixing

\displaystyle F_{k0}=-\frac{6}{C_{0}^{2}{\tilde{\alpha}}^{2}}\left(bJ_{(g)0}+C_{0}^{2}{\tilde{\alpha}}^{2}J_{(f)0}\right),

(160)

and $k_{3}(t)$ is an integration constant. In the limit of $r\to\infty$ , the perturbed metrics behave as

	$\displaystyle g_{ta}\to-\omega_{0}r^{2}-\frac{2J_{(g)0}}{r},\qquad g_{ra}\to-\sqrt{\frac{3}{\Lambda_{(g)}}}\omega_{0}r,$		(161)
	$\displaystyle f_{ta}\to-k_{3}(t)r^{2}-\frac{1}{r}\left(2J_{(f)0}+\frac{k_{0}^{\prime}(t)}{3}\right),\qquad f_{ra}\to-\sqrt{\frac{3}{\Lambda_{(g)}}}\omega_{0}r.$		(162)

For the above solution, the gauge-invariant quantities (146)-(149) reduce to

\displaystyle G_{r}=-\frac{k_{0}}{r^{3}},\qquad G_{t}=\frac{6(J_{(f)0}-J_{(g)0})+k_{0}^{\prime}(t)}{3r^{3}}-\left(\omega_{0}-k_{3}(t)\right),\qquad G_{h}=\frac{6J_{(g)0}}{r^{4}},\qquad G_{k}=\frac{6J_{(f)0}}{r^{4}}.

(163)

The leading terms $\omega_{0}$ and $k_{3}(t)$ are fixed by the boundary conditions at the spatial infinity. Then, the constants $J_{(g)0}$ and $J_{(f)0}$ correspond to the angular momenta in the physical and fiducial sectors. Thus, the general solution for the $\ell=1$ mode in the fiducial sector deviates from the slow-rotation approximation of the Kerr $-$ de Sitter metric. Instead, in the case where the fiducial metric is given by the slow-rotation limit of the Kerr $-$ de Sitter metric, the general solution in the physical sector deviates from the slow-rotation limit of the Kerr $-$ de Sitter metric.

By setting $k_{0}=0$ , then $G_{r}=0$ , we obtain

\displaystyle K_{t}=H_{t}+\frac{2}{r^{3}}\left(\frac{1-N^{r}N^{r}_{(f)}}{1-(N^{r})^{2}}J_{(g)0}-J_{(f)0}\right)+\left(\frac{1-N^{r}N^{r}_{(f)}}{1-(N^{r})^{2}}\omega_{0}-k_{3}(t)\right).

(164)

In the special case where $N^{r}_{(f)}=N^{r}$ and $J_{(f)0}=J_{(g)0}$ ,

\displaystyle K_{t}=H_{t}+\left(\omega_{0}-k_{3}(t)\right).

(165)

which coincide up to the difference between $\omega_{0}$ and $k_{3}(t)$ fixed by the boundary conditions at the spatial infinity. Namely, in this case, we obtain the identical Kerr $-$ de Sitter metrics in the slow-rotation limit in both the sectors.

V Conclusions

As a continuation of our previous work [30], we have studied the dynamical processes of spherically symmetric systems in MTBG. In particular, we have focused on gravitational collapse of pressureless dust and odd-parity perturbations of static and spherically symmetric Schwarzschild $-$ de Sitter black holes. Throughout the paper, we have focused on the self-accelerating branch satisfying the condition (29).

In Sec. III, we have found the exact solutions describing gravitational collapse of pressureless dust, where the interior spacetime geometries are written with the spatially flat and closed FLRW universes, respectively, and the exterior vacuum solutions are described by the Schwarzschild solutions with specific time slicings. For simplicity, we foliated the interior geometry by homogeneous and isotropic spacetimes. For a spatially flat interior universe, we foliated the exterior geometry by a time-independent spatially flat space, while for a spatially curved interior universe, we foliated the exterior geometry by a time-independent spacetime with deficit solid angle. Despite the rather restrictive choice of foliations, we have successfully found interesting classes of exact solutions that represent gravitational collapse in MTBG.

While in the spatially flat case the matter surfaces start to collapse from the spatial infinities with zero initial velocities or from finite distances with nonzero initial velocities, in the spatially closed case they started from finite radii with vanishing initial velocities. The spatially closed case corresponds to an extension of the Oppenheimer-Snyder model to the case of MTBG. Since the Lagrange multipliers are trivial and continuous across the interfaces of the collapsing matter in both the physical and fiducial sectors, the junction conditions across them remain the same as those in the two copies of GR.

In the spatially flat case, under the above-mentioned tuning of the initial conditions of the collapse, we have obtained solutions in which gravitational collapse in the physical and fiducial sectors proceed independently and the gravitational radii in the two sectors may be different. On the other hand, in the spatially closed case, under certain tuning of the matter energy densities and Schwarzschild radii between the two sectors, we have found exact solutions of gravitational collapse. Needless to say, these tunings reflect the fact that we restrict our considerations to a rather specific choice of time slicings. It is expected that relaxing the restriction on the choice of time slicings, e.g. to those with time-dependent spatial metrics, should allow for more general solutions. Furthermore, since in general the Birkhoff theorem does not hold in MTBG, yet-unknown non-Schwarzschild exterior solutions may be found and then even wider classes of global solutions representing gravitational collapse may be constructed. It is also worthwhile studying spherically-symmetric inhomogeneous dust collapse, i.e., analog of the Lemaître-Tolman-Bondi model [42, 43, 44], in the context of MTBG and investigating the homogeneous limit. Moreover, it would be also important to investigate deviations from the spherical GR metric solutions in the weak-field limit and relate them to the parametrized post-Newtonian (PPN) parameters (see e.g. Ref. [3]). Since the absence of the Birkhoff theorem is due to the presence of the effective energy-momentum tensor arising from the variation of the interaction and constraint terms of the Lagrangian (5) and (II.2), in the weak-field limit the deviation from the GR solution would be proportional to $m^{2}$ , which may give rise to a new bound on the graviton mass parameter $m$ with the Solar System experiments. We leave the explicit PPN formulation for future work.

In Sec. IV, we studied odd-parity perturbations of the Schwarzschild $-$ de Sitter black holes in the self-accelerating branch of MTBG. As the background solutions, we considered the Schwarzschild $-$ de Sitter solutions written in the spatially flat coordinates. In the equations of motion for the odd-parity perturbations, there were the contributions of the interaction terms and the Lagrange multipliers in addition to those from the two copies of the Einstein-Hilbert terms. We distinguished the dipolar mode with the angular multipole moment $\ell=1$ and the higher-multipolar modes $\ell\geq 2$ . In order to discuss the odd-parity perturbations with $\ell\geq 2$ , we introduced the master variables. For the modes $\ell\geq 2$ , the system of the odd-parity perturbations in the self-accelerating branch is composed of the two dynamical modes and the two shadowy modes. The dynamical modes correspond to the two master variables of the metric perturbations, while the shadowy modes correspond to the gauge-invariant perturbation of the Lagrange multiplier and one of the two metric perturbations in the angular directions. We note that the other odd-parity metric perturbation in the angular directions is a gauge mode and can be set to be zero. The equation of motion for one of the two shadowy modes can be easily integrated as it is not coupled to the other modes, and the other shadowy mode is sourced by the two dynamical modes and the first shadowy mode. The two dynamical modes are also coupled to each other and sourced by the shadowy modes.

In the case where the ratio of the lapse functions between the physical and fiducial sectors is equal to $C_{0}$ , which is a constant determined by the parameters of the theory, the two dynamical modes are decoupled. Since they are still sourced by one of the shadowy modes, their behavior is different from that in the two copies of GR and depends on the boundary conditions of the shadowy mode. In the case where the ratio of the lapse functions between the physical and fiducial sectors differs from $C_{0}$ , the two dynamical modes are coupled to each other and sourced by both the shadowy modes. In the high frequency and short wavelength limits, we showed that the perturbations do not suffer from ghost or gradient instabilities. On the other hand, in the dipolar sector of $\ell=1$ , we found that the two copies of the slow-rotation limit of the Kerr $-$ de Sitter metrics, in general, are not a solution in the self-accelerating branch of MTBG, unless the masses and angular momenta of black holes and the effective cosmological constants are tuned to be the same in both the sectors. Therefore, deviation from GR is expected for rotating black holes in the self-accelerating branch of MTBG. Readers should also refer to the last paragraph in Sec. IV of [37] for a similar statement in the context of MTMG.

The analysis of the even-parity black hole perturbations in the self-accelerating branch will be left for future work. The analysis of the same issues in the normal branch, i.e., gravitational collapse and black hole perturbations, will also be left for future work. We expect that the spacetime dynamics in the normal branch would differ from that in the case of the two copies of GR at both the background and perturbation levels. We hope to come back to these issues in a future publication.

ACKNOWLEDGMENTS

The work of A.D.F. was supported by Japan Society for the Promotion of Science Grants-in-Aid for Scientific Research No. 20K03969. M.M. was supported by the Portuguese national fund through the Fundação para a Ciência e a Tecnologia in the scope of the framework of the Decree-Law 57/2016 of August 29, changed by Law 57/2017 of July 19, and the Centro de Astrofísica e Gravitação through the Project No. UIDB/00099/2020. M.M. also would like to thank Yukawa Institute for Theoretical Physics for their hospitality under the Visitors Program of FY2022. The work of S.M. was supported in part by Japan Society for the Promotion of Science Grants-in-Aid for Scientific Research No. 17H02890, No. 17H06359, and by World Premier International Research Center Initiative, MEXT, Japan. M.O. would like to thank Yukawa Institute for Theoretical Physics for its institutional support.

Appendix A Equations of motion for the odd-parity perturbations in the spatially flat coordinates in MTBG

We show the equations of motion for the odd-parity perturbations about the Schwarzschild $-$ de Sitter solutions written in the spatially flat coordinates in the self-accelerating branch of MTBG.

A.1 The $\ell\geq 2$ modes

The equations of motion for the gauge invariants $\bar{H}_{t}$ , $\bar{H}_{r}$ , $\bar{K}_{t}$ , $\bar{K}_{r}$ , and $\bar{K}_{3}$ in the odd-parity perturbations for the $\ell\geq 2$ modes are given by

$\displaystyle 0$	$\displaystyle=$	$\displaystyle\dot{\bar{H}}_{r}^{\prime}-N^{r}(r)\bar{H}_{r}^{\prime\prime}-r\bar{H}_{t}^{\prime\prime}+\frac{3}{r}\dot{\bar{H}}_{r}$	(166)
		$\displaystyle-\frac{2}{r}\left(N^{r}(r)+rN^{r}{}^{\prime}(r)\right)\bar{H}_{r}^{\prime}-4\bar{H}_{t}^{\prime}+\frac{1}{r}\left(\ell-1\right)\left(\ell+2\right)\bar{H}_{t}+\left(\frac{2N^{r}(r)}{r^{2}}-\frac{2N^{r}{}^{\prime}(r)}{r}-N^{r}{}^{\prime\prime}(r)\right)\bar{H}_{r},$	(166)
$\displaystyle 0$	$\displaystyle=$	$\displaystyle\ddot{\bar{H}}_{r}-2N^{r}(r)\dot{\bar{H}}_{r}^{\prime}+\left(N^{r}(r)\right)^{2}\bar{H}_{r}^{\prime\prime}-r\dot{\bar{H}}_{t}^{\prime}+rN^{r}(r)\bar{H}_{t}^{\prime\prime}-\left(N^{r}{}^{\prime}(r)+\frac{2}{r}N^{r}(r)\right)\dot{\bar{H}}_{r}$	(167)
		$\displaystyle+\frac{2N^{r}(r)}{r}\left(N^{r}(r)+rN^{r}{}^{\prime}(r)\right)\bar{H}_{r}^{\prime}+4N^{r}(r)\bar{H}_{t}^{\prime}$
		$\displaystyle+\left[\frac{(\ell+2)(\ell-1)}{r^{2}}+(2C_{0}^{3}c_{1}+3C_{0}^{2}c_{2}-c_{4})m^{2}+N^{r}(r)\left(N^{r}{}^{\prime\prime}(r)+\frac{6}{r}N^{r}{}^{\prime}(r)\right)\right]\bar{H}_{r}$
		$\displaystyle+\frac{m^{2}C_{0}^{2}\left(b-1\right)\left(C_{0}c_{1}+c_{2}\right)}{2}\left(\bar{H}_{r}-\bar{K}_{r}-\frac{r}{2}\bar{K}_{3}^{\prime}\right)+\frac{m^{2}C_{0}\left(C_{0}c_{1}+c_{2}\right)\left(1-\beta C_{0}\right)}{2}r\Lambda^{\prime},$
$\displaystyle 0$	$\displaystyle=$	$\displaystyle\dot{\bar{K}}_{r}^{\prime}-N^{r}_{(f)}(r)\bar{K}_{r}^{\prime\prime}-r\bar{K}_{t}^{\prime\prime}+\frac{3}{r}\dot{\bar{K}}_{r}$	(168)
		$\displaystyle-\frac{2}{r}\left(N_{(f)}^{r}(r)+rN_{(f)}^{r}{}^{\prime}(r)\right)\bar{K}_{r}^{\prime}-4\bar{K}_{t}^{\prime}+\frac{1}{r}\left(\ell-1\right)\left(\ell+2\right)\bar{K}_{t}+\left(\frac{2N_{(f)}^{r}(r)}{r^{2}}-\frac{2N_{(f)}^{r}{}^{\prime}(r)}{r}-N_{(f)}^{r}{}^{\prime\prime}(r)\right)\bar{K}_{r},$	(168)
$\displaystyle 0$	$\displaystyle=$	$\displaystyle\ddot{\bar{K}}_{r}-2N^{r}_{(f)}(r)\dot{\bar{K}}_{r}^{\prime}+\left(N^{r}_{(f)}(r)\right)^{2}\bar{K}_{r}^{\prime\prime}-r\dot{\bar{K}}_{t}^{\prime}+rN^{r}_{(f)}(r)\bar{K}_{t}^{\prime\prime}-\left(\frac{2N^{r}_{(f)}(r)}{r}+N^{r}_{(f)}{}^{\prime}(r)\right)\dot{\bar{K}}_{r}$	(169)
		$\displaystyle+\frac{2N^{r}_{(f)}}{r}\left(N^{r}_{(f)}(r)+rN^{r}_{(f)}{}^{\prime}(r)\right)\bar{K}_{r}^{\prime}+4N^{r}_{(f)}(r)\bar{K}_{t}^{\prime}$
		$\displaystyle+\left[\frac{b^{2}(\ell+2)(\ell-1)}{r^{2}}-b^{2}\frac{(c_{0}C_{0}^{2}+2C_{0}c_{1}+c_{2})m^{2}}{{\tilde{\alpha}}^{2}}+N^{r}_{(f)}(r)\left(N^{r}_{(f)}{}^{\prime\prime}(r)+\frac{6}{r}N^{r}_{(f)}{}^{\prime}(r)\right)\right]\bar{K}_{r}$
		$\displaystyle-\frac{m^{2}b\left(b-1\right)(C_{0}c_{1}+c_{2})}{2{\tilde{\alpha}}^{2}}\left(\bar{H}_{r}-\bar{K}_{r}-\frac{1}{2}r\bar{K}_{3}^{\prime}\right)-\frac{m^{2}b\left(c_{1}C_{0}+c_{2}\right)\left(1-\beta C_{0}\right)}{2C_{0}{\tilde{\alpha}}^{2}}r\Lambda^{\prime},$
$\displaystyle 0$	$\displaystyle=$	$\displaystyle m^{2}\left[(1-b)C_{0}\left(r^{2}{\bar{K}_{3}}^{\prime\prime}+4r{\bar{K}_{3}}^{\prime}-(\ell+2)(\ell-1)\bar{K}_{3}+2r\left(\bar{K}_{r}^{\prime}-\bar{H}_{r}^{\prime}\right)+6\left(\bar{K}_{r}-\bar{H}_{r}\right)\right)\right.$	(170)
		$\displaystyle\left.+2(1-C_{0}\beta)\left(r^{2}\Lambda^{\prime\prime}+4r\Lambda^{\prime}-\left(\ell+2\right)\left(\ell-1\right)\Lambda\right)\right],$	(170)
$\displaystyle 0$	$\displaystyle=$	$\displaystyle m^{2}\left[r^{2}{\bar{K}_{3}}^{\prime\prime}+4r{\bar{K}_{3}}^{\prime}-(\ell+2)(\ell-1)\bar{K}_{3}+2r\left(\bar{K}_{r}^{\prime}-\bar{H}_{r}^{\prime}\right)+6\left(\bar{K}_{r}-\bar{H}_{r}\right)\right].$	(171)

A.2 The dipolar mode

The equations of motion for the odd-parity perturbations for the dipolar $\ell=1$ mode are given by

$\displaystyle 0$	$\displaystyle=$	$\displaystyle\dot{H}_{r}^{\prime}-N^{r}(r)H_{r}^{\prime\prime}-rH_{t}^{\prime\prime}+\frac{3}{r}\dot{H}_{r}$	(172)
		$\displaystyle-\frac{2}{r}\left(N^{r}(r)+rN^{r}{}^{\prime}(r)\right)H_{r}^{\prime}-4H_{t}^{\prime}+\left(\frac{2N^{r}(r)}{r^{2}}-\frac{2N^{r}{}^{\prime}(r)}{r}-N^{r}{}^{\prime\prime}(r)\right)H_{r},$	(172)
$\displaystyle 0$	$\displaystyle=$	$\displaystyle\ddot{H}_{r}-2N^{r}(r)\dot{H}_{r}^{\prime}+\left(N^{r}(r)\right)^{2}H_{r}^{\prime\prime}-r\dot{H}_{t}^{\prime}+rN^{r}(r)H_{t}^{\prime\prime}-\left(N^{r}{}^{\prime}(r)+\frac{2}{r}N^{r}(r)\right)\dot{H}_{r}$	(173)
		$\displaystyle+\frac{2N^{r}(r)}{r}\left(N^{r}(r)+rN^{r}{}^{\prime}(r)\right)H_{r}^{\prime}+4N^{r}(r)H_{t}^{\prime}$
		$\displaystyle+\left[(2C_{0}^{3}c_{1}+3C_{0}^{2}c_{2}-c_{4})m^{2}+N^{r}(r)\left(N^{r}{}^{\prime\prime}(r)+\frac{6}{r}N^{r}{}^{\prime}(r)\right)\right]H_{r}$
		$\displaystyle+\frac{m^{2}C_{0}^{2}\left(b-1\right)\left(C_{0}c_{1}+c_{2}\right)}{2}\left(H_{r}-K_{r}\right)+\frac{m^{2}C_{0}\left(C_{0}c_{1}+c_{2}\right)\left(1-\beta C_{0}\right)}{2}r\Lambda^{\prime},$
$\displaystyle 0$	$\displaystyle=$	$\displaystyle\dot{\bar{K}}_{r}^{\prime}-N^{r}_{(f)}(r)K_{r}^{\prime\prime}-rK_{t}^{\prime\prime}+\frac{3}{r}\dot{K}_{r}$	(174)
		$\displaystyle-\frac{2}{r}\left(N_{(f)}^{r}(r)+rN_{(f)}^{r}{}^{\prime}(r)\right)K_{r}^{\prime}-4K_{t}^{\prime}+\left(\frac{2N_{(f)}^{r}(r)}{r^{2}}-\frac{2N_{(f)}^{r}{}^{\prime}(r)}{r}-N_{(f)}^{r}{}^{\prime\prime}(r)\right)K_{r},$	(174)
$\displaystyle 0$	$\displaystyle=$	$\displaystyle\ddot{K}_{r}-2N^{r}_{(f)}(r)\dot{K}_{r}^{\prime}+\left(N^{r}_{(f)}(r)\right)^{2}K_{r}^{\prime\prime}-r\dot{K}_{t}^{\prime}+rN^{r}_{(f)}(r)K_{t}^{\prime\prime}-\left(\frac{2N^{r}_{(f)}(r)}{r}+N^{r}_{(f)}{}^{\prime}(r)\right)\dot{K}_{r}$	(175)
		$\displaystyle+\frac{2N^{r}_{(f)}}{r}\left(N^{r}_{(f)}(r)+rN^{r}_{(f)}{}^{\prime}(r)\right)K_{r}^{\prime}+4N^{r}_{(f)}(r)K_{t}^{\prime}$
		$\displaystyle+\left[-b^{2}\frac{(c_{0}C_{0}^{2}+2C_{0}c_{1}+c_{2})m^{2}}{{\tilde{\alpha}}^{2}}+N^{r}_{(f)}(r)\left(N^{r}_{(f)}{}^{\prime\prime}(r)+\frac{6}{r}N^{r}_{(f)}{}^{\prime}(r)\right)\right]K_{r}$
		$\displaystyle-\frac{m^{2}b\left(b-1\right)(C_{0}c_{1}+c_{2})}{2{\tilde{\alpha}}^{2}}\left(H_{r}-K_{r}\right)-\frac{m^{2}b\left(c_{1}C_{0}+c_{2}\right)\left(1-\beta C_{0}\right)}{2C_{0}{\tilde{\alpha}}^{2}}r\Lambda^{\prime},$
$\displaystyle 0$	$\displaystyle=$	$\displaystyle m^{2}\left(C_{0}c_{1}+c_{2}\right)\left(1-\beta C_{0}\right)\left[r\left(H_{r}^{\prime}-K_{r}^{\prime}\right)+3\left(H_{r}-K_{r}\right)\right].$	(176)

References

Clifton et al. [2012] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, Modified Gravity and Cosmology, Phys. Rept. 513, 1 (2012), arXiv:1106.2476 [astro-ph.CO] .
Berti et al. [2015] E. Berti et al., Testing General Relativity with Present and Future Astrophysical Observations, Class. Quant. Grav. 32, 243001 (2015), arXiv:1501.07274 [gr-qc] .
Will [2014] C. M. Will, The Confrontation between General Relativity and Experiment, Living Rev. Rel. 17, 4 (2014), arXiv:1403.7377 [gr-qc] .
Berti et al. [2018a] E. Berti, K. Yagi, H. Yang, and N. Yunes, Extreme Gravity Tests with Gravitational Waves from Compact Binary Coalescences: (II) Ringdown, Gen. Rel. Grav. 50, 49 (2018a), arXiv:1801.03587 [gr-qc] .
Berti et al. [2018b] E. Berti, K. Yagi, and N. Yunes, Extreme Gravity Tests with Gravitational Waves from Compact Binary Coalescences: (I) Inspiral-Merger, Gen. Rel. Grav. 50, 46 (2018b), arXiv:1801.03208 [gr-qc] .
Fierz and Pauli [1939] M. Fierz and W. Pauli, On relativistic wave equations for particles of arbitrary spin in an electromagnetic field, Proc. Roy. Soc. Lond. A 173, 211 (1939).
Boulware and Deser [1972] D. G. Boulware and S. Deser, Can gravitation have a finite range?, Phys. Rev. D 6, 3368 (1972).
de Rham et al. [2011] C. de Rham, G. Gabadadze, and A. J. Tolley, Resummation of Massive Gravity, Phys.Rev.Lett. 106, 231101 (2011), arXiv:1011.1232 [hep-th] .
Hassan and Rosen [2012] S. Hassan and R. Rosen, Bimetric gravity from ghost-free massive gravity, J. High Energy Phys. 2012 (02), 126, arXiv:1109.3515 [hep-th] .
Yamashita et al. [2014] Y. Yamashita, A. De Felice, and T. Tanaka, Appearance of Boulware–Deser ghost in bigravity with doubly coupled matter, Int. J. Mod. Phys. D 23, 1443003 (2014), arXiv:1408.0487 [hep-th] .
de Rham et al. [2014] C. de Rham, L. Heisenberg, and R. H. Ribeiro, Ghosts and matter couplings in massive gravity, bigravity and multigravity, Phys. Rev. D 90, 124042 (2014), arXiv:1409.3834 [hep-th] .
Gumrukcuoglu et al. [2015] A. E. Gumrukcuoglu, L. Heisenberg, S. Mukohyama, and N. Tanahashi, Cosmology in bimetric theory with an effective composite coupling to matter, JCAP 04, 008, arXiv:1501.02790 [hep-th] .
De Felice and Mukohyama [2016a] A. De Felice and S. Mukohyama, Minimal theory of massive gravity, Phys. Lett. B 752, 302 (2016a), arXiv:1506.01594 [hep-th] .
De Felice and Mukohyama [2016b] A. De Felice and S. Mukohyama, Phenomenology in minimal theory of massive gravity, JCAP 04, 028, arXiv:1512.04008 [hep-th] .
Gumrukcuoglu et al. [2012] A. E. Gumrukcuoglu, C. Lin, and S. Mukohyama, Cosmological perturbations of self-accelerating universe in nonlinear massive gravity, JCAP 03, 006, arXiv:1111.4107 [hep-th] .
De Felice et al. [2012] A. De Felice, A. E. Gümrükçüoğlu, and S. Mukohyama, Massive gravity: nonlinear instability of the homogeneous and isotropic universe, Phys.Rev.Lett. 109, 171101 (2012), arXiv:1206.2080 [hep-th] .
Fasiello and Tolley [2012] M. Fasiello and A. J. Tolley, Cosmological perturbations in Massive Gravity and the Higuchi bound, JCAP 11, 035, arXiv:1206.3852 [hep-th] .
De Felice and Mukohyama [2017] A. De Felice and S. Mukohyama, Graviton mass might reduce tension between early and late time cosmological data, Phys. Rev. Lett. 118, 091104 (2017), arXiv:1607.03368 [astro-ph.CO] .
Hagala et al. [2021] R. Hagala, A. D. Felice, D. F. Mota, and S. Mukohyama, Non-linear dynamics of the minimal theory of massive gravity, Astron. Astrophys. 653, A148 (2021), arXiv:2011.14697 [astro-ph.CO] .
De Felice et al. [2021a] A. De Felice, S. Mukohyama, and M. C. Pookkillath, Minimal theory of massive gravity and constraints on the graviton mass, JCAP 12 (12), 011, arXiv:2110.01237 [astro-ph.CO] .
De Felice et al. [2018a] A. De Felice, F. Larrouturou, S. Mukohyama, and M. Oliosi, Black holes and stars in the minimal theory of massive gravity, Phys. Rev. D 98, 104031 (2018a), arXiv:1808.01403 [gr-qc] .
De Felice et al. [2021b] A. De Felice, F. Larrouturou, S. Mukohyama, and M. Oliosi, Minimal Theory of Bigravity: construction and cosmology, JCAP 04, 015, arXiv:2012.01073 [gr-qc] .
Garcia-Saenz et al. [2021] S. Garcia-Saenz, A. Held, and J. Zhang, Destabilization of Black Holes and Stars by Generalized Proca Fields, Phys. Rev. Lett. 127, 131104 (2021), arXiv:2104.08049 [gr-qc] .
Silva et al. [2022] H. O. Silva, A. Coates, F. M. Ramazanoğlu, and T. P. Sotiriou, Ghost of vector fields in compact stars, Phys. Rev. D 105, 024046 (2022), arXiv:2110.04594 [gr-qc] .
Demirboğa et al. [2022] E. S. Demirboğa, A. Coates, and F. M. Ramazanoğlu, Instability of vectorized stars, Phys. Rev. D 105, 024057 (2022), arXiv:2112.04269 [gr-qc] .
Manita et al. [2022] Y. Manita, K. Aoki, T. Fujita, and S. Mukohyama, Spin-2 dark matter from anisotropic Universe in bigravity, (2022), arXiv:2211.15873 [gr-qc] .
De Felice et al. [2018b] A. De Felice, D. Langlois, S. Mukohyama, K. Noui, and A. Wang, Generalized instantaneous modes in higher-order scalar-tensor theories, Phys. Rev. D 98, 084024 (2018b), arXiv:1803.06241 [hep-th] .
De Felice et al. [2020] A. De Felice, A. Doll, and S. Mukohyama, A theory of type-II minimally modified gravity, JCAP 09, 034, arXiv:2004.12549 [gr-qc] .
De Felice et al. [2021c] A. De Felice, S. Mukohyama, and K. Takahashi, Nonlinear definition of the shadowy mode in higher-order scalar-tensor theories, JCAP 12 (12), 020, arXiv:2110.03194 [gr-qc] .
Minamitsuji et al. [2022] M. Minamitsuji, A. De Felice, S. Mukohyama, and M. Oliosi, Static and spherically symmetric general relativity solutions in minimal theory of bigravity, Phys. Rev. D 105, 123026 (2022), arXiv:2204.08217 [gr-qc] .
De Felice et al. [2021d] A. De Felice, S. Mukohyama, and M. C. Pookkillath, Addressing $H_{0}$ tension by means of VCDM, Phys. Lett. B 816, 136201 (2021d), arXiv:2009.08718 [astro-ph.CO] .
De Felice et al. [2021e] A. De Felice, A. Doll, F. Larrouturou, and S. Mukohyama, Black holes in a type-II minimally modified gravity, JCAP 03, 004, arXiv:2010.13067 [gr-qc] .
De Felice et al. [2022] A. De Felice, K.-i. Maeda, S. Mukohyama, and M. C. Pookkillath, Comparison of two theories of Type-IIa minimally modified gravity, Phys. Rev. D 106, 024028 (2022), arXiv:2204.08294 [gr-qc] .
Oppenheimer and Snyder [1939] J. R. Oppenheimer and H. Snyder, On Continued gravitational contraction, Phys. Rev. 56, 455 (1939).
Kanai et al. [2011] Y. Kanai, M. Siino, and A. Hosoya, Gravitational collapse in Painleve-Gullstrand coordinates, Prog. Theor. Phys. 125, 1053 (2011), arXiv:1008.0470 [gr-qc] .
Blau [2017] M. Blau, Lecture notes on general relativity (2017).
De Felice et al. [2019] A. De Felice, F. Larrouturou, S. Mukohyama, and M. Oliosi, On the absence of conformally flat slicings of the Kerr spacetime, Phys. Rev. D 100, 124044 (2019), arXiv:1908.03456 [gr-qc] .
Markovic and Shapiro [2000] D. Markovic and S. Shapiro, Gravitational collapse with a cosmological constant, Phys. Rev. D 61, 084029 (2000), arXiv:gr-qc/9912066 .
De Felice et al. [2011] A. De Felice, T. Suyama, and T. Tanaka, Stability of Schwarzschild-like solutions in f(R,G) gravity models, Phys. Rev. D 83, 104035 (2011), arXiv:1102.1521 [gr-qc] .
Lipschutz [2013] S. Lipschutz, Schaum’s Outline of Theory and Problems of Linear Algebra, 5th ed. (McGraw-Hill, New York, 2013).
Takahashi and Motohashi [2021] K. Takahashi and H. Motohashi, Black hole perturbations in DHOST theories: master variables, gradient instability, and strong coupling, JCAP 08, 013, arXiv:2106.07128 [gr-qc] .
Lemaitre [1933] G. Lemaitre, The expanding universe, Annales Soc. Sci. Bruxelles A 53, 51 (1933).
Tolman [1934] R. C. Tolman, Effect of imhomogeneity on cosmological models, Proc. Nat. Acad. Sci. 20, 169 (1934).
Bondi [1947] H. Bondi, Spherically symmetrical models in general relativity, Mon. Not. Roy. Astron. Soc. 107, 410 (1947).

	$\displaystyle n_{\alpha}^{(g)}m^{\alpha}_{(g)}\Big{\|}_{(T_{(g)},R_{(g)})}=0,\qquad n_{\alpha}^{(f)}m^{\alpha}_{(f)}\Big{\|}_{(T_{(f)},R_{(f)})}=0,$
	$\displaystyle e^{\alpha}_{(g)\tau}m_{\alpha}^{(g)}\Big{\|}_{(T_{(g)},R_{(g)})}=-1,\qquad e^{\alpha}_{(f)\tau}m_{\alpha}^{(f)}\Big{\|}_{(T_{(f)},R_{(f)})}=-1,$
	$\displaystyle e^{\alpha}_{(g)a}m_{\alpha}^{(g)}\Big{\|}_{(T_{(g)},R_{(g)})}=0,\qquad e^{\alpha}_{(f)a}m_{\alpha}^{(f)}\Big{\|}_{(T_{(f)},R_{(f)})}=0.$		(95)