Spherical scalar collapse in a type-II minimally modified gravity

We adopt the spherically symmetric ansatz

where $\alpha$ and $\beta$ now relate to the lapse and the shift, respectively, and $\Phi$ is the areal radius. In parallel, we also write

where $d\Omega_{2}^{2}$ is the metric of the unit $2$-sphere. In order to simplify the equations of motion, we will use the trace of the extrinsic curvature that now reads

where $\partial_{\perp}=(1/\alpha)(\partial_{t}-\beta\partial_{r})$. Writing $K^{i}_{j}\equiv\gamma^{ik}K_{kj}$, we also introduce the following variables

Given this ansatz, the position of an AH $r=r_{\rm AH}$, can be found by solving²²2A derivation of this condition can be found in appendix A. For more details, please consult subsection III B in [31] and/or subchapters 2.4 and 5.1.7 in [32].

where

In terms of the introduced variables, the equations of motion can then be divided into two sets: one set with the dependence on the Lagrange multipliers $\lambda$ and $\tilde{\lambda}$ removed, and the other set determining $\lambda$ and $\tilde{\lambda}$.

The first set of equations of motion, i.e. those independent of the Lagrange multipliers $\lambda$ and $\tilde{\lambda}$, consists of constraints, dynamical equations, and non-dynamical equations. The constraints are

where equation (14) comes from the Hamiltonian constraint, equation (15) from the momentum constraint, and equation (16) follows from the variation with respect to $\tilde{\lambda}$. The dynamical equations are

Two non-dynamical equations are derived from the definitions of $(a,Q,K)$ and give the first spatial derivative of the lapse and the shift as

The remaining non-dynamical equations are obtained by requiring that the time derivatives of the constraints be consistent with the spatial derivatives of the dynamical equations as

Finally, the extra equations that determine the Lagrange multipliers $\lambda$ and $\tilde{\lambda}$ are

We are interested in evolving the nine variables $(\psi,P,\Phi,Q,\phi,\alpha,\beta,K,a)$, and note that $\lambda$ and $\tilde{\lambda}$ do not appear in any of the equations (14)-(25). Hence, determining them will not be necessary for the later numerical integration; we shall not consider them anymore.

In actual numerical studies, we restrict our considerations to the situation where the solution far from the center approaches a flat spacetime in the standard Minkowski slicing without excitations of the scalar field $\psi$. In vacuum (or for $\psi=\text{const.}$), the theory admits a flat spacetime in the standard Minkowski slicing without spontaneous breaking of the spatial diffeomorphism invariance (and thus with $\tilde{\lambda}=0$) if and only if $\phi=\phi_{0}$, where $\phi_{0}$ is a constant that simultaneously satisfies

One can consider equation (27) as the definition of $\phi_{0}$ and equation (28) as the condition setting the effective cosmological constant to zero. The latter can be achieved by tuning a constant in $V(\phi)$. Since equation (16) says that $\phi$ is independent of $r$, the assumed asymptotic condition that the solution should approach the flat spacetime far from the center implies that $\phi=\phi_{0}$ everywhere, where $\phi_{0}$ is defined by equation (27), and that the effective cosmological constant should be set to zero as in equation (28).

Upon introducing this change into equations (14)-(25), the equations of motion simplify. We note that the previously dynamical equation (21) for $\phi$ downgrades to a non-dynamical one for $a$, explicitly

Because of this, we no longer include equation (25) in the set of equations to be solved; it automatically follows from other equations.

Furthermore, the equation (24) for $K$ simplifies to

which can be integrated once to give $\partial_{r}K(t,r)=f_{1}(t)/[\alpha(t,r)\Phi(t,r)]$, where $f_{1}(t)$ is an arbitrary function of $t$. The regularity of the solution at the center of spherical symmetry, where $\Phi(t,r)=0$, requires that $\alpha(t,r=0)\neq 0$ and that $\partial_{r}K(t,r=0)=0$. Therefore, the regularity of the center sets $f_{1}(t)=0$ and thus $\partial_{r}K(t,r)=0$ and $K(t,r)=f_{0}(t)$, where $f_{0}(t)$ is an arbitrary function of $t$. As stated previously, we assume that the solution far from the center approaches a flat spacetime in the standard Minkowski slicing without excitations of the scalar field $\psi$. This in particular implies that $K(t,r)$ should vanish at the spatial infinity, meaning that $f_{0}(t)=0$ and that $K(t,r)=0$. In this way we have successfully integrated out the so-called shadowy mode, avoiding the corresponding boundary value problem.

We also point out that, as $K$ has now reduced to a constant (actually zero), equation (13) can be re-written

In addition, $K$ vanishing implies from equation (26) that $\lambda$ also becomes $0$. Accordingly, the solutions to the equations (14)-(25) after this simplification will, in fact, be exact GR solutions; we are interested in whether or not an AH appears, and if it does, before or after a breakdown of the time foliation and formation of a singularity. We reiterate that in GR, a change of the spacetime foliation could reverse the order of these events. However, different foliations in VCDM correspond to physically different solutions. Any singularity or breakdown of the time foliation appearing in this theory cannot be circumvented by a change of the time slicing.

For actual numerical simulation, one can downgrade $Q$ from a dynamical variable to a non-dynamical one, abandoning equation (20) and instead using equation (15) to determine $Q$. The simplifications in the previous subsection along with such a downgrade give the equations of motion in their final form. Only the following constraint equation remains

while the dynamical equations reduce to

and the non-dynamical equations become

where equation (39) comes from equation (15). For convenience, equation (32) can be reformulated to define the quantity $\mathcal{C}$ as

In this case, if equation (32) is indeed well respected, the quantity $\mathcal{C}$ remains identically null. Any deviations from $0$, which is bound to happen once numerically implemented, will allow us to watch and quantify the convergence of the numerical integration (see section IV).

We need to be careful about the boundary condition at the center of spherical symmetry, where $\Phi=0$. By redefinition of the radial coordinate $r$, we set

so that the center of spherical symmetry is at $r=0$ and that we consider the region $r\geq 0$. For the regularity of the center, we require ($\psi,P,Q,\alpha$) to be even functions of $r$ and $(\Phi,\beta,a)$ to be odd functions of $r$. Moreover, equations (32) and (39) respectively imply that

By Taylor expanding the variables $(\psi,P,\Phi,Q,\alpha,\beta,a)$ with respect to $r$ around $r=0$ and inserting these, one can confirm that the right-hand sides of all relevant equations (32)-(39) are well-defined in the limit $r\to+0$, implying that they can also be Taylor expanded with respect to $r$. As a consistency check, one can show that the right-hand sides of equations (32), (33), (34), (37), and (38) are even functions of $r$, while the right-hand sides of equations (35), (36), and (39) are odd functions of $r$. Furthermore, by using the $r\to+0$ limits of equations (37) and (32) respectively, one can also show that the first $r$-derivative of the right-hand side of equation (35) vanishes at $r=0$.

We now introduce an outer boundary at $r=r_{\rm b}$ and hereafter consider the region $0\leq r\leq r_{\rm b}$. In order to evolve $(\psi,P,\Phi)$ with the dynamical equations (33)-(35), we need to impose appropriate boundary conditions on the outer boundary as well. For example, if $r_{\rm b}$ is large enough then we can set

Up to here, the overall normalization of $\alpha$ is not fixed, corresponding to the fact that the theory enjoys time reparametrization symmetry. We fix it by demanding that

so that the time $t$ agrees with the proper time measured by an observer at the outer boundary. The equations (37) and (38) for $\beta$ and the acceleration $a$ do not require additional boundary conditions; the fact that $\beta$ and $a$ are odd functions of $r$ already imposes $\beta(t,r=0)=0=a(t,r=0)$.

We consider the following initial conditions at the initial time $t=t_{0}$:

and

where $\psi_{0}(r)$ is a free function of $r$, which we, for concreteness, take to be an even Gaussian wave packet of the form

The constants $A$, $r_{0}$ and $s$ determine the amplitude, starting position, and width of the initial wave packet respectively. Having defined $\psi_{0}$, the constraint (32) implies that $\Phi_{0}(r)$ is the solution to

with the boundary conditions given in subsection III.4, i.e.

The dynamics of the collapse are fully described by three dynamical equations for $(\psi,P,\Phi)$ and four non-dynamical ones for $(\alpha,\beta,a,Q)$³³3These form a strongly hyperbolic system of partial differential equations, since after having integrated out the shadowy mode, they were shown to be the same as in GR.. We solve the system described in subsection III.3 with the boundary conditions and initial conditions formulated in subsections III.4 and III.5. In order to remain consistent, the non-dynamical equations must be solved at every time evaluation. Starting from the initial time $t=t_{0}$, the initial condition must satisfy the constraint (32) and the non-dynamical equations (36)-(39). We then evolve $(\psi,P,\Phi)$ to the next time step $t=t_{1}$ with the dynamical equations (33)-(35). In order to obtain $(\alpha,\beta,a,Q)$ at $t=t_{1}$, we integrate the non-dynamical equations (36)-(39) once more, now on the time slice $t=t_{1}$. At this point, the numerical accuracy of the solution may be gauged by checking how well the data on the time slice $t=t_{1}$ satisfies the constraint (32). The numerical scheme can then be advanced to the next time step by repeating this procedure.

Unless stated otherwise, the initial field $\psi$ is given by an even Gaussian wave packet (defined in equation (47)) of height $A=0.15$ and a width characterized by $s=4$ positioned at $r_{0}=10$. An outer boundary $r_{\mathrm{b}}$ of the system is placed far from the origin $r=0$ and the initial wave packet, to prevent it from disturbing the collapsing process and to best simulate an asymptotically flat spacetime. In our case, $r_{\mathrm{b}}=80$ was used. Larger or slightly smaller values yield equivalent results. Space and time are uniformly discretized into intervals of size $\Delta r$ and $\Delta t$, respectively. We decide here to keep the ratio $\Delta t/\Delta r=0.2$ constant. A smaller ratio would also have been acceptable, but if $\Delta t$ is too large relative to $\Delta r$, the system may become unstable.

To prevent any numerical instabilities from jeopardizing the computations, we also adopt two additional safety measures. Firstly — and most importantly — we use a series expansion near the boundaries to approximate the first and last integration points of the nondynamical fields (see appendix B). Secondly, we include Kreiss-Oliger dissipation terms [33]. These two measures serve to mitigate any possible emergent high-frequency modes near the boundaries and around the scalar field profile. As is standard, the artificial dissipation is chosen to be two orders higher than the convergence order of the spatial numerical derivatives.

The simulation was executed with several independently written codes, using different integration schemes, and all results matched. As an extra sanity check, the quantity $\mathcal{C}$, defined in equation (40), was continuously and carefully surveilled to ensure it remained sufficiently small throughout the simulation. A more quantitative examination of $\mathcal{C}$ was also carried out, as figure 1 demonstrates. While the shape of $\mathcal{C}$ remains unchanged (left plot), its overall amplitude is inversely proportional to the square of the time step $\Delta t$ or, equivalently, the spatial gap $\Delta r$ (right plot). The quadratic convergence is expected. All codes — two of which are, in fact, the same codes as in [34, 35] and one of the two was also used in [36, 37] — exploit second-order methods to integrate the non-dynamical equations in space, and at least second-order methods for the integration in time. Explicitly, one used the trapezoidal method in space and the iterated Crank-Nicolson method with two iterations in time, while the other two used Heun’s method in space and higher order methods, mainly RK4, in time.

We now turn to the physical quantities thusly obtained⁴⁴4We remind the reader that as a consequence of having integrated out the shadowy mode, the integrated equations (32)-(39) are the same as in GR. Naturally, the results are then the same as in GR as far as the foliation-independent parts are concerned. Unlike in GR, however, the foliation-dependent parts are also physical in VCDM and thus will be investigated with care.. Figure 2 shows the areal radius $\Phi$ against the proper distance $r$ for some selected times $t$. Notice that one of the displayed times, in this figure and the subsequent ones, is $t\approx 14.1$. As we shall explain in subsection IV.3, this corresponds to the AH formation time. At $t=0$, we observe that $\Phi\approx r$. As time proceeds, it only retains this behavior near $r=0$ (in accordance with subsection III.4), falling and settling down into a plateau elsewhere. The height of this plateau, $\Phi_{\text{pl.}}\approx 0.65$, echoes in the other quantities, as the subsequent figures show. Henceforth, we shall use $\Phi$ as a radial coordinate to illustrate the evolution of all other quantities.

The other two dynamical quantities $\psi$ and $P$ are given for different times in figure 3. For the prescribed initial conditions, the scalar field $\psi$ starts by splitting into two parts. One moves out to $r=\infty$ and vanishes, while the other falls to the origin $r=0$ under its self-gravity and eventually settles near $r=0$. Naturally, this foreshadows an AH formation (see subsection IV.3). A similar behavior is observed for the related variable $P$ ($\equiv\partial_{\perp}\psi$, as defined in equation (11)).

In parallel, the results for the nondynamical fields are given in figure 4. One can observe how the lapse $\alpha$ and the shift $\beta$ collapse to zero below the threshold value of $\Phi_{\text{pl.}}$. Far from the origin, an asymptotically flat value is manifestly recovered for both quantities. The variables $a$ and $Q$ simultaneously spike around the same threshold value of $\Phi_{\text{pl.}}$.

The quantity $g^{\Phi\Phi}$ (defined in equation (31)) is depicted for different times in the left plot of figure 5. As the behavior of $\psi$ already suggested (see figure 3), the collapse does in fact yield an AH, which occurs when $g^{\Phi\Phi}$ crosses zero. This condition is satisfied at $t\approx 14.1$, as the right plot of figure 5 shows. Furthermore, in figure 4, one can deduce that $\alpha$ and $\beta$ have not yet collapsed down to $0$ at the time of the crossing, implying that the AH indeed forms before any singularity or breakdown of the time slicing does. Eventually, as the lapse and shift evolve towards zero inside the AH, $g^{\Phi\Phi}$ settles with a minimum around $\Phi_{\text{pl.}}\approx 0.65$ and an outermost crossing of zero at $\Phi\approx 0.88$. This constitutes the main achievement of this work, namely that the gravitational collapse of a massless scalar field in VCDM results in the manifestation of an AH before the breakdown of the time foliation. The final stage of the collapse is a static black hole.

Until now, we have kept fixed the initial scalar field parameters of $\psi_{0}$ as defined in equation (47) to $A=0.15$ and $s=4$. However, depending on the values assigned to these two parameters, an AH may or may not form. The parameter space where an AH does appear is illustrated in figure 6. A decrease in $A$ and/or an increase of $s$ may both “dilute” the scalar field enough to allow it to bounce back around the origin before any AH appears or the time slicing breaks down. An example of this is shown in the left plot of figure 7, where we have chosen $s=8$ and plotted $\psi$ at different times throughout the evolution. The right plot of the same figure makes the non-formation of an AH evident: $g^{\Phi\Phi}$ never crosses zero.