On the nature of Hawking's incompleteness for the Einstein-vacuum equations: The regime of moderately spatially anisotropic initial data

In the mathematical physics literature, there are heuristic arguments, going back three decades, suggesting that for an open set of initially smooth solutions to the Einstein-vacuum equations in high dimensions, stable, approximately monotonic curvature singularities can dynamically form along a spacelike hypersurface. In this article, we study the Cauchy problem and give a rigorous proof of this phenomenon in sufficiently high dimensions, thereby providing the first constructive proof of stable curvature blowup (without symmetry assumptions) along a spacelike hypersurface as an effect of pure gravity. Our proof applies to an open subset of regular initial data satisfying the assumptions of Hawking's celebrated"singularity"theorem, which shows that the solution is geodesically incomplete but does not reveal the nature of the incompleteness. Specifically, our main result is a proof of the dynamic stability of the Kasner curvature singularity for a subset of Kasner solutions whose metrics exhibit only moderately (as opposed to severely) spatially anisotropic behavior. Of independent interest is our method of proof, which is more robust than earlier approaches in that i) it does not rely on approximate monotonicity identities and ii) it accommodates the possibility that the solution develops very singular high-order spatial derivatives, whose blowup rates are allowed to be, within the scope of our bootstrap argument, much worse than those of the base-level quantities driving the fundamental blowup. For these reasons, our approach could be used to obtain similar blowup results for various Einstein-matter systems in any number of spatial dimensions for solutions corresponding to an open set of moderately spatially anisotropic initial data, thus going beyond the nearly spatially isotropic regime treated in earlier works.


Introduction
Hawking's celebrated "singularity" theorem (see, for example, [59,Theorem 9.5.1]) shows that an interestingly large set of initial data for the Einstein-vacuum 1 equations leads to geodesically incomplete solutions. The chief drawback of this result is that the proof is by contradiction and therefore does not reveal the nature of the incompleteness; see Subsect. 1.4 for further discussion. In this article, for an open 2 subset of regular initial data in high spatial dimensions that satisfy the assumptions of Hawking's theorem, we show that the solution's incompleteness is due to the formation of a Big Bang, that is, the formation of a curvature singularity along a spacelike hypersurface. For more detailed statements of our main results, readers can jump ahead to Theorem 1.1 for a rough summary or to Theorem 11.1 for the precise versions.
Before proceeding, we note that the Einstein-vacuum equations in D spatial dimensions are where Ric is the Ricci curvature of the Lorentzian spacetime metric g (which has signature (−, +, ⋯, +)), and R = Ric α α is 3 its scalar curvature. Throughout, M denotes the 1 + D-dimensional spacetime manifold, on which the system (1.1) is posed. In this article, the spacetimes under study are cosmological, i.e., they have compact spacelike Cauchy hypersurfaces. In particular, the solutions that we study are such that M = I × T D , where I is an interval of time and D is the standard D-dimensional torus (i.e., [0, 1] D with the endpoints identified).
Our work provides the first constructive proof of stable curvature blowup along a spacelike hypersurface as an effect of pure gravity (i.e., without the presence of matter) for Einstein's equations in more than one spatial dimension without symmetry assumptions. The present work can be viewed as complementary to our previous works [50,51,53], in which we proved similar results in the presence of scalar field or stiff fluid 4 matter. As we explain below, out of necessity, we had to devise a new and more robust analytic framework to treat the Einstein-vacuum equations, since many of the special structures that we exploited in [50,51,53] are not available in the vacuum case.
For the singularities that we study here, the blowup is rather "controlled" in the sense that the solutions exhibit approximately monotonic behavior as the singularity is approached. For example, relative to an appropriately constructed time coordinate t, the Kretschmann scalar Riem αβγδ Riem αβγδ blows up like Ct −4 as t ↓ 0 (see Lemma 10.1 for the precise statement), where Riem is the Riemann curvature of g. A main theme of this paper is that the monotonic behavior is not just a curiosity, but rather it lies at the heart of our analysis. Heuristic evidence for the existence of a large family of (non-spatially homogeneous) approximately monotonic spacelike-singularity-forming Einsteinvacuum solutions for large D goes back more than 30 years to the work [26], which was preceded by other works of a similar vein that we review in Subsubsect. 1.4.1. In [26] and in the present work as well, the largeness of D is used only to ensure the existence of background solutions with certain quantitative properties (in our case Kasner solutions whose exponents verify the condition (1.7d)); for the Einstein-vacuum equations in low space dimensions, the only obstacle to the existence of such solutions is the Hamiltonian constraint equation (1.2a). Given such a background solution, the rest of our analysis is essentially dimensionally independent. 5 Although our main theorem applies only when D ≥ 38, our approach here is of interest in itself since it is more robust compared to prior works on stable blowup for various Einstein-matter systems, and since it has further applications, for example in three spatial dimensions; see the end of Subsect. 1.2 for further discussion on this point.
1.1. The evolution problem and the initial data. Before further describing our results, we first provide some background material on the evolution problem for Einstein's equations. The foundational works [15,16] collectively imply that the Einstein-vacuum equations (1.1) have an evolution problem formulation in which all sufficiently regular initial data verifying the constraint equations (1.2a)-(1.2b) launch a unique 6 maximal classical solution (M (M ax) , g (M ax) ), known as the maximal globally hyperbolic development of the data. Below we will discuss the evolution equations. For now, we recall that the initial data are (Σ 1 ,g,k), where Σ 1 is as D-dimensional manifold,g is a Riemannian metric on Σ 1 , andk is a symmetric type 0 2 Σ 1 -tangent tensorfield. The subscript "1" on Σ 1 emphasizes that in our main theorem, Σ 1 will be identified with a hypersurface of constant time 1.g andk represent, respectively, the first and second fundamental form of Σ 1 , viewed as a Riemannian submanifold of the spacetime to be constructed. The constraints are the well-known In analyzing solutions to (1.1), we will use the well-known constant-mean-curvature-transported spatial coordinates gauge. The corresponding PDE system involves hyperbolic evolution equations for the first and second fundamental forms of the constant-time hypersurface Σ t coupled to an elliptic PDE for the lapse n > 0, which is defined by g(∂ t , ∂ t ) = −n 2 . See Subsubsect. 2.1.1 for a review of this gauge. Here we only note that the elliptic PDE for n is essential for synchronizing the singularity across space. To employ this gauge in the context of our main theorem, we assume that the initial data verify the constant-mean-curvature (CMC from now on) condition k a a = −1. 3) should not be viewed as a further restriction on the initial data since for the near-Kasner solutions that we study, we can always construct a CMC hypersurface lying near the initial data hypersurface. This can be achieved by making minor modifications to the arguments given in [50], where we constructed a CMC hypersurface in a closely related context.

Some additional context and preliminary comments on the new ideas in the proof.
The aforementioned results [15,16], while philosophically of great importance, reveal very little information about the nature of the maximal globally hyperbolic development. In our main theorem, we derive sharp information about the maximal globally hyperbolic development of an open set of nearly spatially homogeneous initial data (i.e., initial data that vary only slightly in space at each fixed time) given along the D-dimensional torus Σ 1 = T D , where D ≥ 38. The largeness of D allows us to consider initial data that are only moderately (as opposed to severely) spatially anisotropic, in the sense that the data are close to a Kasner solution (1.5) whose Kasner exponents verify (1.7d) (see below for further discussion). Broadly speaking, the main analytic themes of our work can be summarized as follows (see Subsect. 1.5 for a more detailed outline of our proof): The amount of spatial anisotropy exhibited by the solutions under study is tied to the strength of various nonlinear error terms that depend on spatial derivatives. Below a certain threshold, the error terms become subcritical (in the sense of the strength of their singularity) with respect to the main terms. This allows us to give a perturbative proof of stable blowup through a bootstrap argument, where the blowup is driven by "ODE-type" singular terms that are linear in time derivatives, and, at least at the low derivative levels, the nonlinear spatial derivative terms are strictly weaker than the singular linear time derivative terms. Put differently, at the low derivative levels, the spatial derivative terms remain negligible throughout the evolution, a phenomenon that, in the general relativity literature, is sometimes referred to as asymptotically velocity term dominated (AVTD) behavior. In contrast, at the high derivative levels, the spatial derivatives can be very singular; bounding the maximum strength of this singular behavior and showing that it does not propagate down into the low derivative levels together constitute the main technical challenge of the proof.
We stress that in the absence of matter, spatially homogeneous and isotropic Big Bang solutions do not exist; they are precluded by the Hamiltonian constraint equation (1.2a). Hence, it is out of necessity that our main theorem concerns spatially anisotropic solutions. In our previous works [50,51,53], we proved related stable blowup results for the Einstein-scalar field and Einstein-stiff fluid systems in three spatial dimensions. In those works, the presence of matter allowed for the existence of spatially homogeneous, spatially isotropic Big-Bang-containing solutions, specifically, the famous Friedmann-Lemaître-Robertson-Walker (FLRW) solutions, whose perturbations we studied. The approximate spatial isotropy of the perturbed solutions lied at the heart of the analysis of [50,51,53], and we therefore had to develop a new approach to handle the moderately spatially anisotropic solutions under study here. We now quickly highlight the main new contributions of the present work; in Subsect. 1.5, we provide a more in depth overview of how they fit into our analytical framework.
• In [50,51,53], the presence and the precise structure of the matter model was used in several key places, leaving open the possibility that the blowup was essentially "caused" by the matter. In contrast, our present work shows that for suitable initial data (which exist for D ≥ 38), stable curvature blowup can occur in general relativity as an effect of pure gravity. • Our previous works [50,51,53] relied on approximate monotonicity identities, which were L 2type integral identities in which the special structure of the matter models and the spatially isotropic nature of the FLRW background solutions were exploited to exhibit miraculous cancellations. The cancellations were such that dangerous singular error integrals with large coefficients were replaced, up to small error terms, with coercive ones. That is, the identities led to the availability of friction-type spacetime integrals that were used to control various error terms up to the singularity. The net effect was that we were able to prove that the very high spatial derivatives of the solution are not much more singular than the low-order derivatives; this was a fundamental ingredient in the analytic framework that we used to control error terms. In contrast, in the present article, we completely avoid relying on approximate monotonicity identities, which might not even exist for the kinds of solutions that we study here. This leads, at the high spatial derivative levels, to very singular energy estimates near the Big Bang, which our proof is able to accommodate; see Subsect. 1.5 for an outline of the main ideas behind the estimates. For these reasons, our approach here is significantly more robust compared to [50,51,53] and could be used, for example, to substantially enlarge the set of initial data for the Einstein-scalar field and Einstein-stiff fluid systems in three spatial dimensions that are guaranteed to lead to curvature blowup. More precisely, it could be used to prove stable blowup for perturbations of generalized Kasner solutions (i.e., non-vacuum Kasner solutions) to these systems that exhibit moderate spatial anisotropy. We stress that it might not be possible to extend these stable blowup results to Einstein-matter systems that feature g-timelike characteristics, such as the Euler-Einstein system under a general equation of state; in [50,51,53], the analysis relied on the fact that for the matter models studied, the characteristics of the system are exactly the null hypersurfaces of g, a feature that is tied to the following crucial structural property: the absence of time-derivative-involving error-terms in various equations. Readers can consult [51] for further discussion on this point.

1.3.
Kasner solutions and a rough summary of the main theorem. As we have mentioned, in our main theorem, we consider initial data given on where D is the standard D-dimensional torus (i.e., [0, 1] D with the endpoints identified) and, relative to the coordinates that we use in proving our main results, Σ 1 will be identified with a spacelike hypersurface of constant time 1, i.e., Σ 1 = {1} × T D . Our assumption on the topology of Σ 1 is for convenience and is not of fundamental importance; we expect that similar results hold for other spatial topologies, much in the same way that the blowup results of [50,51] in the case of T 3 spatial topology were extended to the case of S 3 spatial topology in [53] (see Subsubsect. 1.4.4 for further discussion of these works). Our main theorem addresses the evolution of perturbations of the initial data (at time 1) of the Kasner solutions [32] g ∶= −dt ⊗ dt +g, (t, x) ∈ (0, ∞) × T D , (1.5) whereg 6) and the constants q i ∈ [−1, 1], known as the Kasner exponents, are constrained by In the rest of the paper, we will also use the notatioñ to denote the Kasner mixed second fundamental form, i.e.,k i j = − 1 2 (g −1 ) ia ∂ tgaj relative to the coordinates featured in (1.5)-(1.6). Aside from exceptional cases in which one of the q i are equal to 1 and the rest equal to 0, Kasner solutions have Big Bang singularities at {t = 0} where their Kretschmann scalar Riem αβγδ Riem αβγδ blows up like t −4 (see the proof of Lemma 10.1 for a proof of this fact). We now briefly summarize our main results. See Theorem 11.1 for the precise statements. More precisely, perturbations (belonging to a suitable high-order Sobolev space) of the Kasner initial data along Σ 1 = {t = 1} launch a perturbed solution that also has a Big Bang singularity in the past of Σ 1 . In particular, relative to a set of CMC-transported spatial coordinates normalized by k a a = −t −1 , where k i j is the (mixed) second fundamental form of Σ t , the perturbed solution's Kretschmann scalar blows up like Ct −4 . Hence, the past of Σ 1 in the maximal (classical) globally hyperbolic development of the data is foliated by a family of spacelike CMC hypersurfaces Σ t . Furthermore, every past-directed causal geodesic that emanates from Σ 1 crashes into the singular hypersurface Σ 0 in finite affine parameter time. That is, the perturbed solutions are geodesically incomplete to the past, and the incompleteness coincides with curvature blowup. Remark 1.2 (Additional information about the solution). Theorems 1.1 (and 11.1) reveal only the most fundamental aspects of the singularity formation. It is possible to derive substantial additional information about the solution using the estimates that we prove in this paper. For example, one could show that various time-rescaled variables, such as the time-rescaled second fundamental form components tk i j (t, x), converge to regular functions of x as t ↓ 0, which is a manifestation of the AVTD behavior mentioned above. For brevity, we have chosen not to state or prove such additional results in this article since, thanks to the a priori estimates yielded by Prop. 9.2, their statements and proofs are similar to the ones given in [50,Theorem 2]. This bound for q emerges from considering the strength of various error terms; see Subsect. 1.5 for further discussion on this point.

Previous breakdown results for Einstein's equations.
The study of the breakdown of solutions in general relativity was ignited by the classic Hawking-Penrose "singularity" theorems (see e.g. [28,40], and the discussion in [29]), which show that for appropriate matter models and spatial topologies, a compellingly large set of initial conditions (with non-empty interior relative to a suitable function space topology) leads to geodesically incomplete solutions. In particular, Hawking's theorem (see [59,Theorem 9.5.1]) guarantees that under assumptions satisfied by the initial data featured in our main theorem, all past-directed timelike geodesics are incomplete. Although these theorems are robust with respect to the kinds of initial data and matter models to which they apply, they are "soft" in that they do not reveal the nature of the incompleteness, leaving open the possibilities that i) it is tied the blowup of some invariant quantity, such as a spacetime curvature scalar, or ii) it is due to some other more sinister phenomenon, such as the formation of a Cauchy horizon (beyond which the solution cannot be classically extended in a unique sense due to lack of sufficient information). In the wake of the Hawking-Penrose theorems, many authors studied the nature of the breakdown, though, with only a few key exceptions, the works produced were heuristic, numerical, concerned initial data with symmetry, or yielded information only about "special" solutions (as opposed to an open set of solutions corresponding to regular initial data given along a spacelike Cauchy hypersurface). In the next several subsubsections, we review some of these works.

1.4.1.
Heuristic and numerical work. The famous-but-controversial work [11] of Belinskii-Khalatnikov-Lifshitz gave heuristic arguments suggesting that for the Einstein-vacuum equations in three spatial dimensions, the "generic" (in an unspecified sense) solution that breaks down should i) be local near the breakdown, i.e., be well-modeled by a family of Bianchi IX 7 ODE solutions 8 that are parameterized by space; ii) become highly oscillatory near the breakdown-points like the Bianchi IX solutions; and iii) the breakdown points should collectively form a spacelike singularity. We refer readers to the recent monograph [9] for a detailed, modern account of [11] and related works. Although the work [11] has stimulated a lot of research activity, as of the present, it is not clear to what extent the heuristic picture it painted holds true. The picture is almost certainly not true in a generic sense, for the recent breakthrough work [22] shows, assuming only a widely believed (but not yet proven) quantitative version of the stability of the Kerr black hole family of solutions to the Einstein-vacuum equations, that the Cauchy horizon inside the black hole is dynamically stable. Specifically, in [22], the authors showed that the metric is C 0 -extendible past the Cauchy horizon, which is a null hypersurface. This in particular contradicts the vision of spacelike singularities posited in [11]. In the opposite direction (i.e., in accordance with [11]), in three spatial dimensions, Ringström [46] confirmed the oscillatory picture 9 of solutions near singularities for solutions with Bianchi IX symmetry to the Einstein-vacuum equations and to the Euler-Einstein equations under the equations of state p = c 2 s ρ, where p is the pressure, ρ is the density, and the constant 0 < c s < 1 is the speed of sound. However, outside of the class of spatially homogeneous solutions, there are 7 Readers can consult [19] for an overview of the Bianchi IX class and other symmetry classes that we mention later. 8 More precisely, the authors argued that they should be well-modeled by the so-called "Mixmaster" [36] solutions, which are of the form g = −dt 2 is the fully antisymmetric symbol normalized by [123] = 1. 9 Ringström also studied the stiff fluid case c s = 1 in Bianchi IX symmetry and proved monotonic-type singularity formation results similar to the ones we obtained in [50,51,53] and in the present work.
currently no known examples of Einstein-vacuum solutions that exhibit the oscillatory behavior suggested by [11].
There are also heuristic results concerning the existence non-oscillatory solutions to various Einstein-matter systems that form a spacelike singularity. Specifically, in [10], Belinskii and Khalatnikov noted that if one considers the Einstein-scalar field system in three spatial dimensions, then the heuristic arguments given in [11] concerning oscillations no longer seem plausible. They argued that instead, the generic incomplete solution should exhibit monotonic behavior near the breakdown-points, which should still collectively form a spacelike singularity. In [8], Barrow gave a similar heuristic argument for the Einstein-stiff fluid system. Most relevant for our work here is the work [26] that we mentioned at the beginning, in which the authors noted that for D ≥ 10 (i.e., at least 11 spacetime dimensions), there is a nontrivial regime for the Einstein-vacuum equations in which heuristic non-oscillatory arguments along the lines of [8,10] go through. This suggests that indeed, in high spatial dimensions, stable blowup results of the type that we prove in our main theorem should hold. Note that there is a significant gap between the assumption D ≥ 10 needed for the heuristic arguments in [26] and the assumption D ≥ 38 that we make in our main theorem. In Subsect. 1.6, we will further discuss this gap.
We close this subsubsection by noting that the above works and others like them have stimulated a large number of numerical studies that have been designed to probe the validity of the heuristic predictions. We do not attempt to survey the vast literature here, but instead refer to the works [3,12,27,43], which serve as useful starting points for exploring the subject of numerical analysis in the context of singularities in general relativity.
1.4.2. Construction (but not the stability of ) of singular solutions. There are a variety of works showing the existence, but not the stability, of singularity-containing solutions to various Einsteinmatter systems that exhibit the same kind of AVTD behavior 10 near the singularity exhibited by the solutions in our main theorem. Typically, the constructions rely on deriving/solving a Fuchsian PDE system. Roughly speaking, a Fuchsian PDE is a system of the following form (where for convenience we restrict out attention to one spatial dimension): where u is the array of unknowns, A α and B are matrices (symmetric when the symmetric hyperbolic framework for energy estimates is invoked), and f is an array, all of which must satisfy a collection of technical assumptions. A typical Fuchsian analysis is based on splitting the solution as u = u 0 + w, where u 0 is the leading order part and w is an error term that one would like to show is small compared to u 0 as t ↓ 0. The leading order part u 0 must be guessed/solved for using an appropriate ansatz. Here we do not describe how this is typically accomplished; readers can consult [51] for a further description of Fuchsian analysis in the context of singular solutions in general relativity. Although the Fuchsian approach can sometimes be used to show the existence of singular solutions, it is inadequate for treating the true stability problem, i.e., the problem of posing initial data on a regular Cauchy hypersurface {t = const} with const > 0 and then solving the equations all the way down to the singular hypersurface {t = 0}. We now describe some of the Fuchsian-type existence results in more detail. Notable among these is the work of Andersson-Rendall [5], in which they constructed a family of spatially analytic Big-Bang-containing solutions to the Einstein-scalar field and Einstein-stiff fluid systems in three spatial dimensions that exhibit approximately monotonic behavior near the singularities. The family of solutions that the authors constructed was large in the sense that its number of degrees of freedom coincides with the number of free functions in Einstein initial data for the standard Cauchy problem. However, the results of [5] do not show the stability of the blowup under Sobolev-class perturbations of initial data given along a spacelike hypersurface near the expected singularity. Another notable work in the spirit of [5] is [25], in which the authors proved similar results for various Einstein-matter systems in various spatial dimensions, including for the Einstein-vacuum equations in 10 or more spatial dimension (i.e., D ≥ 10 in the notation of the present article). Note that the Einstein-vacuum result with D ≥ 10 supports the heuristic work [26] mentioned in the previous subsubsection.
There are many additional works that yield the construction of (but not the stability of) singularitycontaining solutions to select nonlinear Einstein-matter systems. We do not attempt to exhaustively survey the literature here, but we mention the following ones, which concern various symmetry reduced equations: [1,7,13,14,30,33,41,54].
There are also results in which the authors constructed singular solutions by essentially prescribing "singular data" on the singular hypersurface itself and then solving to the future; see, for example, [6,21,37,38,[55][56][57]. One drawback of this approach is that there are fewer degrees of freedom in the singular data compared to the initial data corresponding to a standard Cauchy problem on a regular spacelike hypersurface. Thus, the solutions constructed in this fashion are not generic; this is explained in more detail in [44, Section 6.1].
We close this subsubsection by mentioning that the Fuchsian techniques behind the above results have applications outside of general relativity. There are techniques for constructing solutions to large classes of hyperbolic Fuchsian PDEs; see, for example, [2]. Readers can consult Rendall's work [42] for a more detailed comparison of many of the results described above as well as application of the Fuchsian framework to prove the existence of singular solutions to the Einstein-vacuum equations with Gowdy symmetry.

1.4.3.
Constructive results under symmetry assumptions. There are a variety of works on symmetry reduced Einstein-matter systems in which the authors give a constructive proof of stable singularity formation. In the spatially homogeneous case, in which the equations reduce to a system of ODEs, there are many results, including the constructive work of Ringström [46] mentioned earlier. We do not attempt to survey the literature here; instead we direct readers to [44,58] for an overview of ODE blowup for solutions to Einstein's equations.
There are also constructive proofs of stable singularity formation for various symmetry reduced Einstein-matter systems in which the equations reduce to a 1 + 1-dimensional system of PDEs; see, for example, [20,31,48,49]. Chief among these are Christodoulou's remarkable works [17,18] on the Einstein-scalar field system in three spatial dimensions in spherical symmetry for 1-or 2-ended asymptotically flat data. In those works, he showed that the maximal globally hyperbolic future developments of generic data are future-inextendible as time-oriented Lorentzian manifolds with a C 0 metric; i.e., the breakdown is severe, at the level of the metric itself. See also the recent works of [34,35] on the spherically symmetric Einstein-Maxwell-(real) scalar field system with asymptotically flat 2-ended initial data, in which the authors proved that the maximal globally hyperbolic future developments of generic data are future-inextendible as time-oriented Lorentzian manifolds with a C 2 metric. This is an especially intriguing result in view of the fact that Dafermos-Rodnianski [23,24] showed that the statement is false for this system if one replaces C 2 with C 0 .

1.4.4.
Constructive results without symmetry assumptions. The only prior works exhibiting stable spacelike singularity formation for solutions to Einstein's equations without symmetry assumptions are [50,53], which are closely related to the present work. In [50,53], we showed that curvature blowup develops along a spacelike hypersurface for open sets of solutions to two Einstein-matter systems: the Einstein-scalar field system and the Einstein-stiff fluid system; see also our related work [51] concerning the linear analysis and Ringström's recent monograph [45], which concerns estimates for solutions to a large family of linear wave equations whose corresponding metrics model the behavior that can occur in solutions to Einstein's equations near cosmological singularities.
Specifically, in [50], we showed that in three spatial dimensions with spatial topology T 3 , the FLRW solution is nonlinearly stable in a neighborhood of its Big Bang singularity. In [53], we proved the same result in the case of S 3 spatial topology, a key new feature being that the solutions are not approximately spatially flat in the S 3 case (although they are nearly spatially homogeneous and isotropic). As we stressed earlier, the analytic framework of [50,51,53] was quite different from that of the present work, due to the availability of L 2 -based approximate monotonicity identities tied to the special structure of the matter models and the spatially isotropic nature of the FLRW background solutions.
1.5. An overview of the proof of the main theorem. In this subsection, we outline the ideas behind the proof of Theorem 11.1. As in our prior works [50,51,53], we analyze solutions relative CMC-transported spatial coordinates, in which the spacetime metric is decomposed into the lapse n and a Riemannian metric g on the constant time hypersurfaces Σ t as follows: (1.9) In such coordinates, the Einstein-vacuum evolution problem consists of the Hamiltonian and momentum constraint equations, hyperbolic evolution equations for the first fundamental form g of Σ t and the second fundamental form k of Σ t , and an elliptic PDE for n along Σ t , supplemented by initial data for g and k given along Σ 1 that satisfy the constraints; see Prop. 2.1 for the details.
Here we only note that the mixed second fundamental form k i j verifies ∂ t g ij = −2ng ia k a j and that we normalize t so that k a a = −t −1 , i.e., Σ t is a hypersurface of constant mean curvature −D −1 t −1 . The main part of the proof is showing that the solution (g, k, n) exists classically for (t, x) ∈ (0, 1] × T D , i.e., long enough to form a curvature singularity. The proofs that the Kretschmann scalar blows up as t ↓ 0 and that the spacetime is geodesically incomplete follow as straightforward consequences of estimates that we use in proving existence on (0, 1]×T D ; we refer readers to Sect. 10 for details on the nature of the breakdown, and we will not discuss them further here.
The main step in proving that the solution exists classically for (t, x) ∈ (0, 1] × T D is to obtain suitable a priori estimates showing that various norms along Σ t do not blow up until t = 0. For this reason, in our discussion here, we describe only the a priori estimates. At the heart of the proof lies the following task, whose importance we describe below: Showing that for perturbed solutions, the t-rescaled type 1 1 spatial Ricci components tRic i j are, at each fixed spatial point x ∈ T D , integrable in time over t ∈ (0, 1]. Here and throughout, Ric denotes the Ricci curvature of g.
The above task is essentially a quantified version of the following idea, which has its origins in the heuristic works discussed in Subsubsect. 1.4.1: For near-Kasner initial data, we can prove stable blowup in regimes where we can prove that time-derivative terms dominate spatial derivative terms in the equations, i.e., when we can prove that the AVTD behavior described in Subsubsect. 1.2 occurs. In practice, to prove the time-integrability of tRic i j and the many other estimates that we need to close the proof, we rely on a bootstrap argument involving various norms that capture the following behavior: the high-level derivatives are allowed to be substantially more singular than the low-level derivatives. Here we will not describe the logical flow of our bootstrap argument in detail, but rather only describe how the various estimates fit together consistently. Moreover, we will not focus on the "smallness assumptions" (i.e., near-Kasner assumptions) on the initial data that we need in our detailed proof, but rather only on the main feature of the analysis: the various powers of t that arise; in this subsection, we will simply denote all quantities that can be estimated in terms of the near-Kasner initial data by "Data". We stress already that our proof relies on commuting the evolution equations with up to N spatial derivatives, where N has to be chosen sufficiently large in a manner that we explain below.
We now recall that we are studying perturbations of a Kasner solution (1.5) whose exponents satisfy (1.7d) (see Subsect. 2.3 for a proof that such Kasner solutions exist when D ≥ 38). We also recall that for the Kasner solution (g,k,ñ), we haveñ ≡ 1, whileg andk are respectively given by (1.6) and (1.7c). Note also that the Kasner solution is spatially flat and thus the Ricci tensor of g vanishes. For the perturbed solution, to obtain the desired time integrability of the components tRic i j described in the previous paragraph, it clearly suffices to prove that Ric i j ≲ t −p for some constant p < 2. (1.10) We stress that (1.10) is essentially the same as the heuristic estimates featured in the works [8,10,10,26], and that the estimate was verified by the solutions that we studied in [50,51,53]. Before describing how we prove (1.10), we first outline the two main consequences that it affords: (1) n − 1 → 0 as t ↓ 0.
(2) The components tk i j remain bounded as t ↓ 0. The proof of (1) follows from a simple application of the maximum principle to the elliptic PDE t 2 g ab ∇ a ∇ b n = (n − 1) + t 2 nR, where R = Ric a a is the scalar curvature of g; see the proof of (5.1) for the details. The proof of (2) essentially follows from freezing the spatial point and integrating the following evolution equation (see equation (6.3)) from time t to time 1: ∂ t (tk i j ) = tRic i j + ⋯, and from using a few additional estimates that allow one to show that, like tRic i j , the error terms denoted by ⋯ are integrable over t ∈ (0, 1]; see the proof of Prop. 6.1 for the details. We now return to the crucial issue of proving that Ric i j ≲ t −p for some constant p < 2, a bound that is tied to all aspects of the proof. To achieve this, we first, in view of (1.7d), fix a constant q such that max i=1,⋯,D q i < q < 1 6.
To control Ric i j , we rely on the following estimates, whose proof we will describe below: max Note that the bounds (1.11) represent an absolute worst-case scenario, in which all components of g and g −1 are allowed to be slightly more singular than the most singular component of the background Kasner spatial metric and its inverse. As will become clear, the bounds (1.11) are the most fundamental ones in our analysis. One might think of (1.11) as allowing for the "complete mixing" of the Kasner exponents in the perturbed solutions; this is the most glaring spot in the proof that has potential for improvements in future studies. To control Ric i j , we first express it in terms of the Christoffel symbols of the transported spatial coordinates (see (2.12b) for the precise expression). For the solutions under study, the most singular term in the component Ric i j is not a top-order term, but rather lower-order terms (i.e., the last two products on RHS (2.12b)) that have the following schematic form: g −3 (∂g) 2 , where ∂ denotes the spatial gradient with respect to the transported spatial coordinates. An interpolation argument, which heavily relies on (1.11) and which we explain below, yields that for large N (where we again stress that N is the maximum number of times that we commute the equations with spatial derivatives), the low-level spatial derivatives of the components of g, including ∂g, are only slightly more singular than the components of g itself. Thus, in view of (1.11), we see that g −3 (∂g) 2 is only slightly more singular than (t −2q ) 5 = t −10q . Since q < 1 6, we conclude that the product g −3 (∂g) 2 is less singular than t −2 , consistent with the desired bound (1.10).
Remark 1.4. Since g −3 (∂g) 2 is fifth-order in (g, ∂g), the discussion above suggests that the proof should close assuming only q < 1 5. However, in the top-order energy estimates, we encounter some below-top-order error terms that seem to prevent the proof from closing unless we assume q < 1 6; we encounter such error terms, for example, in the proof of (7.13b).
Having outlined how to obtain the desired bound for Ric i j , we can, as we described above, show that as t ↓ 0, n − 1 → 0 and tk i j remains bounded. Then, given these bounds for n and k, we can integrate the evolution equations ∂ t g ij = −2ng ia k a j and ∂ t g ij = 2ng ia k j a and use a near-Kasner assumption on the initial data to obtain, through standard arguments, the desired estimates (1.11) (see Prop. 6.1 and its proof for the precise statements).
We now return to the issue of the interpolation argument mentioned above, which we used to show, for example, that ∂g is only slightly more singular than g. To appreciate the role of interpolation, it is essential to already know that the best energy estimates we are able to obtain allow for the following rather singular behavior: The top-order derivatives of g can be as large (in L 2 ) as Data × t −(A+1) , where A > 0 is a large universal constant, independent of the maximum number of times (i.e., N) that we commute the equations with spatial derivatives. That is, under appropriate bootstrap assumptions, we can prove only that 11 (1.12) Then standard Sobolev interpolation (see (4.4)) yields the following bound for the components of ∂g: . (1.13) From (1.13), we see that even if the top-order homogeneous norm ∂g ḢN (Σt) is as singular as Data × t −(A+1) , as long as we take N to be sufficiently large (relative to A and D), the singular nature of ∂g L ∞ (Σt) will not be much worse than that of g L ∞ (Σt) , i.e., no worse than t −(2q+α) , where α is small; see the beginning of Subsect. 4.4 for further discussion of this issue.
It remains for us to discuss the top-order energy estimates for g and k. In reality, these must be complemented with top-order elliptic estimates for n, but we will ignore this (relatively standard) issue in this subsection in order to condense our summary of the proof. Let ⃗ I be a top-order spatial derivative multi-index, i.e., ⃗ I = N. Using the evolution equations and integration by parts, and using appropriate bootstrap assumptions to control error terms, we are able to derive an estimate of the following form, valid for t ∈ (0, 1] (see Prop. 7.1 for the details and Subsect. 3.3 for the definition of the norm ⋅ L 2 g (Σt) ): (1.14) ds + ⋯, where ⋯ denotes error terms that, while they require some care to handle, are less singular. We stress the following points: • The constant C * on RHS (1.14) is universal, i.e., independent of N and A. Roughly, C * is generated by the coefficients of the most singular linear terms in the evolution equations and the elliptic PDE for n; the most singular terms are all tied, either directly or indirectly, to the top-order derivatives of g and k. In our prior works [50,51,53], we were able to show that C * is small or vanishing, thanks to the approximate monotonicity identities that we mentioned earlier in the paper. For the solutions under study here, C * can be large (and we do not bother to carefully track the precise value of C * ). • We inserted the time weights t A+1 "by hand" into the energies in equation (1.14). We note that when proving (1.14) (roughly, by taking the time derivative of the LHS and integrating by parts), one encounters terms in which ∂ t falls on the weights. Roughly, this leads to the integrals preceded by the factor of −2A on RHS (1.14) (note that t ∈ (0, 1], which explains why the factors −2A are on the RHS).
The key point is that if we choose A to be sufficiently large, then the factor − {2A − C * } on RHS (1.14) is negative, and the corresponding integral has an overall "friction" sign. In particular, it can be discarded, leaving only the less singular error terms ⋯ on RHS (1.14). A careful analysis of ⋯ allows one to conclude, via Gronwall's inequality, the top-order a priori energy estimate Although the above discussion summarizes the main ideas, in our detailed proof, we encounter several hurdles that we overcome using additional ideas. While conceptually simple, these ideas are somewhat technically involved. We close this subsection by highlighting some of the features of this analysis.
• To control error terms, we rely on a variety of norms. In the next point below, we will shed some light on how we use the different kinds of norms; see Subsects. 3.2 and 3.3 for the precise definitions. For example, when bounding Σ t -tangent tensors, we rely on pointwise norms ⋅ F rame that measure the size of components relative to the transported spatial coordinate frame as well as the more geometric norm ⋅ g , which measures the size of tensors using the dynamic spatial metric g. Similarly, our analysis relies on Sobolev norms for the frame components, such as ⋅ H N F rame (Σt) , as well as more geometric Sobolev norms ⋅ H N g (Σt) .
• Although the interpolation estimates described above are sufficient for controlling various error terms at the low derivative levels, interpolation is not quite sufficient, in itself, for controlling some of the error terms in the high-order energy estimates. For example, our proof of the error term estimate (7.13b) relies not only on interpolation, but also on the structure of Einstein's equations in our gauge. We now explain why this is the case. First, the top-order energy estimates can be obtained only in terms of geometric norms such as ⋅ ḢN g (Σt) , since the basic energy identities involve these kinds of norms. In contrast, for M ≤ N, Sobolev interpolation estimates involving the "background differential operators" ∂ i are most naturally stated and derived in terms of the less geometric norms ⋅ ḢM F rame (Σt) . For tensorfields, the discrepancy between the norms ⋅ ḢM F rame (Σt) and ⋅ ḢM g (Σt) is factors of g and g −1 (where the number of factors depends on the order of the tensorfield), and by (1.11), the two norms can differ in strength by (singular) powers of t −q ; in some cases, these powers of t −q are strong enough so that in the energy estimates, certain below-top-order error terms seem, at first glance, to be more singular than the main top-order terms. If this were the case, then our bootstrap argument would not close. Fortunately, this is not the case, but to prove this, we use an argument that involves deriving energy estimates not only at the very highest level, but also at down-to-two derivatives below top. In deriving these below-top-order energy estimates, we use arguments that lose one derivative, by treating the evolution equations like transport equations along the integral curves of ∂ t with derivativelosing source terms. These transport-type estimates lead to better below-top-order estimates compared to the ones that pure interpolation would afford, and we stress that this approach is viable only because of the structure of Einsteins equations in our gauge. To implement this procedure in a consistent fashion, we rely on a hierarchy of energies that features different t weights at different orders and that involves both geometric norms ⋅ ḢM g (Σt) and coordinate frame norms ⋅ ḢM F rame (Σt) . See Defs. 3.5 and 3.6 for the precise definitions of the t-weighted norms that use to control the solution; we prove that the norms from Defs. 3.5 and 3.6 are uniformly bounded up to the singularity.
• In carrying out our bootstrap argument, we find it convenient to derive, as a preliminary step, estimates showing that the lapse n can be controlled in terms of g and k. This is the content of Sect. 5. We have downplayed these estimates in this subsection since they are based on deriving standard estimates for elliptic equations. One feature worth noting is that for most lapse estimates, we rely on the elliptic PDE t 2 g ab ∇ a ∇ b n = (n − 1) + t 2 nR. This is a "good" equation in the sense that it involves source terms that depend only on spatial derivatives of g, which are less singular than time derivatives of g (as represented by k). However, to obtain the top-order lapse estimates, we first use the Hamiltonian constraint (2.9a) to algebraically replace, in the elliptic lapse PDE, R with terms that depend on k; see equation (5.8) and its proof. This algebraic replacement leads to error terms that can be controlled within the scope of our bootstrap approach, both from the point of view of regularity and from the point of view of the structure of the singular error terms that our framework can accommodate.
1.6. Further comparisons with two related works. In this subsection, we compare and contrast the reasons behind the assumed minimum values of D in the present work and in the aforementioned works [25,26]. We start by recalling that, as we described in Subsubsect. 1.4.1, the work [26] provided heuristic evidence for the existence of a large family of monotonic spacelike singularity-forming Einstein-vacuum solutions whenever D ≥ 10, and that such solution families were constructed in [25] (though stability in the sense of the present article was not proved there). There is a substantial gap between the assumption D ≥ 10 and the assumption D ≥ 38 of the present article. In the remainder of this subsection, we provide an overview of how the assumption D ≥ 10 is used in [25,26], and shed some light on how one might approach the problem of extending the results of the present article to apply to a larger range of D values. Actually, we will focus only on the heuristic work [26] since this allows for a simplified presentation of the main ideas.
In [25,26] and the present work, a crucial step is justifying that, in some sense, the spatial derivative terms in the Einstein-vacuum equations are negligible compared to the time derivative terms near the singularity. For example, in the present work (see in particular Subsect. 1.5), this step is embodied by the estimate (1.10), that is, our proof that the components of the Ricci tensor of g obey the following bound: Ric i j ≲ t −p for some constant p < 2. Similarly, the heuristic arguments given in [26] were designed exactly to yield a bound of this type. From this perspective, a primary analytic difference between the present work and the works [25,26] is that the latter works precisely accounted for, in a tensorial fashion, partial cancellations of singular powers of t that can occur in the products featured in the coordinate expression of the Ricci tensor Ric i j of g (i.e., the products on RHS (2.12b)), at least for metrics that can be interpreted as having "spatially dependent Kasner exponents" (see below for more on this point). In contrast, in the present article, we allow for the possibility that all components of g and g −1 are as singular as t −2q (see, for example, the discussion surrounding (1.11)), and we have therefore ignored the possibility of exploiting such cancellations; this is apparent from the proof outline that we gave in Subsect. 1.5. We now further explain the connection between the work [26] and the notion of a spatial metric having "spatially dependent Kasner exponents." The starting point of [26] is the hope that there are singular solutions to the Einstein-vacuum equations that are somehow well-described by a spacetime metric having "spatially dependent Kasner exponents", that is, a metric of the following form, defined for 12 (t, x) ∈ (0, 1] × T D : where {ω (i) (x)} i=1,⋯,D are a linearly independent set of time-independent one-forms on T D (in particular, a (x)dx a relative to local coordinates), and {q i (x)} i=1,⋯,D are "x-dependent" Kasner exponents, subject to the following "spatially dependent Kasner constraints" (i.e., x-dependent analogs of (1.7a)-(1.7b)): (1.16) Given the above assumptions, the authors of [26] then observe (through straightforward but tedious computations) that for metrics of the form (1.15), the Ricci tensor components of the spatial metric verify as long as the following system of inequalities holds: Note that (1.17) is essentially equivalent to the estimate (1.10) that we rigorously obtain in proving our main results; as we described in Subsect. 1.5, (1.10) is a quantified version of the idea that spatial derivative terms should be negligible. The main conclusions of [26] can be summarized as follows.
In view of (1.16), it is not possible to simultaneously satisfy all inequalities (1.18) when D ≤ 9. However, there does exist an open set of {q i } i=1,⋯,D satisfying (1.18) whenever D ≥ 10. We now stress that the metric g featured in (1.15) does not generally solve the Einstein-vacuum equations. However, it does solve a truncated version of the equations in which all spatial derivative terms, including Ric i j , are discarded; in the mathematical general relativity literature, the truncated system is often referred to as the Velocity Term Dominated (VTD) system. For this reason, the condition (1.17), which is supposed to capture the negligibility of the spatial derivative terms, suggests that one can think of the VTD solution g as "asymptotically solving" the Einstein-vacuum equations of Prop. 2.1 as t ↓ 0. That is, for metrics of the form (1.15), the k-involving products in equation (2.10c) are of size t −2 while, by (1.17), the term Ric i j is less singular. Put differently, the metrics g featured in (1.15) solve a PDE system obtained from the Einstein equations by throwing away terms that can be shown to be small in the sense of (1.17). The work [26] can therefore be viewed as providing a kind of "consistency argument" for metrics (1.15) that satisfy (1.18), i.e., an argument based on ignoring terms in the evolution equations that, for the non-solution g, are small compared to the main terms. For this perspective, it is reasonable to speculate that, for metrics g of the form (1.15) that satisfy (1.17), there might be true Einstein-vacuum solutions lying "close" to g.
We now explain some further connections between the heuristic picture painted in [26] and the results of the present article. As we described in Subsect. 1.5, the estimate (1.10) (which is almost the same as the condition (1.17) from [26]) is a main ingredient needed to show that at each fixed x, ∂ t (tk i j (t, x)) is integrable over the time interval (0, 1]. From this integrability condition, one can easily show not only that tk i j (t, x) remains bounded as t ↓ 0 (which is what we prove in our main theorem), but also that the following stronger result holds: tk i j (t, x) converges uniformly to a function K i j (x) as t ↓ 0; see [50,51] for proofs of these kinds of convergent results in a related context, and see also Remark 1.2. Combining these kinds of convergence estimates with related ones (in particular for the lapse n), one could rigorously prove that for the Einstein-vacuum solutions g under study in this article, g is asymptotic to a metric that behaves like the family of metrics that satisfy (1.15)-(1.18); again, readers can consult [50,51] for proofs of these kinds of results in a related context.
The situation can be summarized as follows: the singular solutions that can be shown to be dynamically stable under our framework are asymptotic (as the singularity is approached) to a metric that is well-described by the family of metrics satisfying (1.15)-(1.18). The glaring question, then, is whether or not, under the assumptions (1.16) and (1.17), the "VTD solutions" g defined by (1.15) are always the asymptotic end-state of a true singularity-forming solution of the Einsteinvacuum equations. A natural starting point for addressing this question would be to try to extend the stable blowup results of the present article to apply to perturbations of all Kasner solutions (1.5) whose Kasner exponents {q i } i=1,⋯,D satisfy (1.18), which in particular would extend our results to the cases D ≥ 10. If such a result is in fact true, then its proof would almost certainly involve observing the kinds of tensorial cancellations that the authors of [26] exploited to derive the bound (1.17), a feat that we did not attempt in the present work. To allow for the observation of these kinds of tensorial cancellations, it is likely that an important step in the proof would be anticipating, in some fashion going beyond our work here, that the perturbed solution is converging towards a metric that is similar to ones of the form (1.15). This is a worthy avenue for future investigation, especially since it is intimately tied to the fundamental question of which terms in Einstein's equations are the ones driving the breakdown of solutions. 1.7. Paper outline. The remainder of the paper is organized as follows.
• In Subsect. 1.8, we summarize the notation and conventions that we use in the rest of the article. • In Sect. 2, we set up the ensuing analysis by providing the Einstein-vacuum equations in CMC-transported spatial coordinates and showing that Kasner solutions verifying our assumptions exist when D ≥ 38.
• In Sect. 3, we define the norms that we use to control solutions and formulate the bootstrap assumptions that we use in our analysis.
• In Sect. 4, we provide some preliminary technical estimates, which are standard. • In Sect. 5, we derive elliptic estimates showing that the lapse n can be controlled, in various norms, in terms of g and k.
• In Sect. 6, we derive preliminary estimates for g and k at the low derivative levels.
• In Sect. 7, we derive preliminary estimates for g and k at the top derivative levels, i.e., our main top-order energy estimates.
• In Sect. 8, we derive preliminary energy estimates for g and k at the near-top-order derivative levels; as we explained near the end of Subsect. 1.5, for technical reasons, we need these estimates to close our bootstrap argument.
• In Sect. 9, we combine the results of Sects. 5-8 to obtain the main technical result of the article: a priori estimates showing that appropriately defined solution norms can be uniformly controlled by the initial data, all the way up to the singularity. • In Sect. 10, we derive some results that describe the breakdown as t ↓ 0, i.e., the curvature blows up and past-directed timelike geodesics terminate with a finite length. • In Sect. 11, we synthesize the results of the previous sections and give a relatively short proof of the main stable blowup theorem.

Notation and conventions.
For the reader's convenience, in this subsection, we provide some notation and conventions that we use throughout the article. Some of the concepts referred to here are not defined until later.
1.8.1. Indices. Greek "spacetime" indices α, β, ⋯ take on the values 0, 1, ⋯, D, while Latin "spatial" indices a, b, ⋯ take on the values 1, 2, ⋯, D. Repeated indices are summed over (from 0 to D if they are Greek, and from 1 to D if they are Latin). We use the same conventions for primed indices such as a ′ as we do for their non-primed counterparts. Spatial indices are lowered and raised with the Riemannian metric g ij and its inverse g ij . 1.8.2. Spacetime tensorfields and Σ t -tangent tensorfields. We denote spacetime tensorfields T µ 1 ⋯µ l ν 1 ⋯νm in bold font. We denote Σ t -tangent tensorfields T a 1 ⋯a l b 1 ⋯bm in non-bold font.
1.8.3. Coordinate systems and differential operators. We often work in a fixed standard local coordinate system (x 1 , ⋯, x D ) on T D . The vectorfields ∂ j ∶= ∂ ∂x j are globally well-defined even though the coordinates themselves are not. Hence, in a slight abuse of notation, we use {∂ 1 , ⋯, ∂ D } to denote the globally defined vectorfield frame. We denote the corresponding dual frame by {dx 1 , ⋯, dx D }. In CMC-transported spatial coordinates, the spatial coordinate functions are transported along the unit normal to Σ t , thus producing a local coordinate system (x 0 , x 1 , x 2 , x 3 ) on manifolds-withboundary of the form (T, 1] × T D , and we often write t instead of x 0 . The corresponding vectorfield frame on Relative to this frame, Kasner solutions are of the form (1.5). The symbol ∂ µ denotes the frame derivative ∂ ∂x µ , and we often write ∂ t instead of ∂ 0 and dt instead of dx 0 . Many of our equations and estimates are stated relative to the frame ∂ µ µ=0,1,⋯,D and dual frame dx µ µ=0,1,⋯,D . We use the notation ∂f to denote the spatial coordinate gradient of the function f relative to the transported spatial coordinates. Similarly, if ω is a Σ t -tangent one-form, then ∂ω denotes the Σ t -tangent type 0 2 tensorfield with components ∂ i ω j relative to the frame described above.
is an array comprising D non-negative integers, then we define the spatial multi-indexed differential operator ∂ ⃗ I by ∂ ⃗ I ∶= ∂ n 1 1 ∂ n 2 2 ⋯∂ n D D . The notation ⃗ I ∶= n 1 + n 2 + ⋯ + n D denotes the order of ⃗ I. Throughout, D denotes the Levi-Civita connection of g. We write to denote a component of the covariant derivative of a tensorfield T (with components T The Christoffel symbols of g, which we denote by Γ λ µ ν , are defined by (1.20) We use similar notation to denote the covariant derivative of a Σ t -tangent tensorfield T (with components T a 1 ⋯am b 1 ⋯bn ) with respect to the Levi-Civita connection ∇ of the Riemannian metric g, i.e., the first fundamental form of Σ t . The Christoffel symbols of g, which we denote by Γ i j k , are defined by Above, the notation " ∫ T D f dx" denotes the integral of f over T D with respect to the measure corresponding to the volume form of the standard Euclidean metric E on T D , which has the components diag(1, 1, ⋯, 1) relative to the coordinate frame described in Subsubsect. 1.8.3. Note that dx is not the canonical integration measure associated to the Riemannian metric g.
All of our Sobolev norms are built out of the (spatial) L 2 norms of scalar quantities (which may be the components of a tensorfield). We define the standard L 2 norm ⋅ L 2 (Σt) as follows: (1.23) For integers M ≥ 0, we define the standard H M norm ⋅ H M (Σt) as follows: We also define the following standard homogeneous analog of (1.24): Finally, we define the Lebesgue norm ⋅ L ∞ (Σt) of scalar functions f in the usual way: In Subsects. 3.2 and 3.3, we will introduce additional norms for tensorfields, many of which are built out of the basic ones from this subsubsection.

Parameters.
• ε is a small "bootstrap parameter" that, in our bootstrap argument, bounds the size of the solution norms; see (3.7). The smallness of ε needed to close the estimates is allowed to depend on the numbers N, A, q, σ, and D. • A ≥ 1 denotes an "exponent parameter" that is featured in the high-order solution norms from Def. 3.6. To close our estimates, we will choose A to be large enough to overwhelm various universal constants C * (see Subsubsect. 1.8.6).
• N denotes the maximum number of times that we commute the equations with spatial derivatives (e.g., k ∈ H N (Σ t ) and g, n ∈ H N +1 (Σ t )). To close our estimates, we will choose N to be sufficiently large in a (non-explicit) manner that depends on A, q, σ, and D.
• 0 < q < 1 6 is a constant, fixed throughout the proof, that bounds the magnitude of the background Kasner exponents.
• σ > 0 is a small constant, fixed throughout the argument, that we use to simplify the proofs of various estimates that "have room in them." • q and σ are constrained by (3.1).
• δ > 0 is a small (N, D)-dependent parameter that is allowed to vary from line to line and that is generated by the estimates of Lemma 4.5. In particular, we use the convention that a sum of two δ's is another δ. The only important feature of δ that we exploit in the proof is the following: at fixed D, we have lim N →∞ δ = 0. In particular, if A is also fixed, then lim N →∞ Aδ = 0.
• C denotes a positive constant that is free to vary from line to line. C can depend on N, A, D, q, and σ, but it can be chosen to independent of all ε > 0 that are sufficiently small in the manner described in Subsubsect. 1.8.5.
• C * denotes a positive constant that is free to vary from line to line and that can depend on D. However, unlike C, C * can be chosen to be independent of N, A, and ε, as long as ε is sufficiently small in the manner described in Subsubsect. 1.8.5. For example, 1 + CN!ε = C * while N! = C and N! σ = C, where C and C * are as above.
• We write v = ∏ R r=1 v r to indicate that the Σ t -tangent tensorfield v is a tensor product, possibly involving contractions, of the Σ t -tangent tensorfields v r . We use this notation only when the precise details of the tensor product are not important for our analysis. We sometimes display the indices of v to indicate its order; for example, the where the coefficients in the linear combination are constants ±C, with C as in Subsubsect. 1.8.6. As above, we sometimes display the indices of v to indicate its order.
where the coefficients in the linear combination are constants ±C * , with C * universal constants enjoying the properties described in Subsubsect. 1.8.6. As above, we sometimes display the indices of v to indicate its order. For where the C * ,i are independent of N and A.

Setting up the Analysis
In this section, we start by providing the Einstein-vacuum equations in the gauge that we use to prove our main theorem. Next, we provide some standard expressions for the curvature tensors of the first fundamental form of Σ t . Finally, show that when D ≥ 38, there exist Kasner solutions that satisfy the Kasner exponent assumptions in our main theorem.
2.1. The Einstein-vacuum equations in CMC-transported spatial coordinates. In this subsection, we recall the Einstein-vacuum equations in CMC-transported spatial coordinates gauge; this is the gauge that we use throughout the article.

Basic ingredients in the setup.
In CMC-transported spatial coordinates, the spacetime metric is decomposed into the lapse n and a Riemannian metric g on the constant time hypersurfaces Σ t (known as the first fundamental form of Σ t ) as follows: (2.1) We use g ab to denote the components of the inverse Riemannian metric g −1 . Above and throughout, t is a time function on the spacetime manifold M that we will describe just below and {x a } a=1,⋯,D are standard (locally defined) "spatial coordinates" on the constant-time hypersurfaces Σ t ∶= {(s, x) ∈ M s = t}, which are diffeomorphic to T D . We sometimes use the alternate notation x 0 = t. We use the notation ∂ µ = ∂ ∂x µ to denote the partial derivative vectorfields corresponding to the above coordinates, and we often write ∂ t in place of ∂ 0 . Note that {∂ a } a=1,⋯,D are a globally defined smooth frame on Σ t , even though the coordinate functions {∂ a } a=1,⋯,D are only locally defined. Note also that the dual co-frame {dx a } a=1,⋯,D is also globally defined and smooth. In view of (2.1), we see that the future-directed unit normalN to Σ t can be expressed aŝ Note thatNx a = 0 for a = 1, ⋯, D. Thus, in the gauge under consideration, the spatial coordinates are transported along the flow lines ofN. The second fundamental form k of Σ t is defined by requiring that the following relation holds for all vectorfields X, Y tangent to Σ t : where D is the Levi-Civita connection of g. It is a standard fact that k is symmetric: For such X, Y , the action of the spacetime connection D can be decomposed into the action of the Levi-Civita connection ∇ of g and k as follows: When analyzing the components of k (and in particular when differentiating the components of k), we will always assume that it is written in mixed form (i.e., type 1 1 form) as k i j with the first index upstairs and the second one downstairs. This convention is absolutely essential for some of our analysis; in the problem of interest to us, the evolution and constraint equations verified by the components k i j have a more favorable structure than the corresponding equations verified by k ij .
Throughout the vast majority of our analysis, we normalize the CMC hypersurfaces Σ t as follows: In order for (2.6) to hold, the lapse has to verify the elliptic equation (2.11). We adopt the following sign convention for the Riemann curvature Riem of g: Similarly, we adopt the following sign convention for the Riemann curvature Riem of g: Proposition 2.1 (The Einstein-vacuum equations in CMC-transported spatial coordinates). In CMC-transported spatial coordinates normalized by k a a = −t −1 , the Einstein-vacuum equations (1.1) take the following form.
• The Hamiltonian and momentum constraint equations verified by g and k are respectively: On the nature of Hawking's incompleteness for the Einstein-vacuum equations where ∇ denotes the Levi-Civita connection of g, R = Ric a a denotes the scalar curvature of g, and Ric denotes the Ricci curvature of g (a precise expression is given in (2.12b)).
• The evolution equations verified by g, g −1 , and k are: Standard expressions for curvature tensors of g. For future use, we note the following standard facts: relative to an arbitrary coordinate system on Σ t (and in particular relative to the transported spatial coordinates that we use in our analysis), the components of the type 0 4 Riemann curvature Riem of g and the type 1 1 Ricci curvature Ric of g can be expressed, respectively, as where Γ ijk ∶= g ja Γ a i k and Γ i j k are the Christoffel symbols of g (see (1.21)). To start, we note that in the case D = 36, the following Kasner exponents satisfy (1.7a)-(1.7b) but just barely fail to satisfy (1.7d): Considering now the case D = 38, we let ǫ > 0 be a small parameter, and we perturb the 36 Kasner exponents from (2.13) by ǫ so that they are bounded in magnitude by < 1 6: (2.14) Notice that by (1.7a)-(1.7b), assuming (2.14), any solution (q 37 , q 38 ) to the following system yields, when complemented with the exponents (2.14), a complete set of Kasner exponents: Using (2.15a) to solve for q 38 in terms of q 37 and then substituting into (2.15b), we obtain the equation q 2 37 − 6ǫq 37 − 6ǫ + 36ǫ 2 = 0, which has the solutions We now observe that for any ǫ > 0 sufficiently small, the corresponding solutions q 37 to (2.16) are real and bounded in magnitude by ≲ √ ǫ. From (2.15a), we deduce that the same statement holds for the corresponding exponent q 38 . In view of (2.14), we see that for ǫ > 0 sufficiently small, the exponents {q i } i=1,⋯,38 constructed in this fashion satisfy (1.7a)-(1.7b) and (1.7d). Moreover, in the cases D ≥ 39, we can complement these 38 Kasner exponents with others as follows: q i = 0 for 39 ≤ i ≤ D. In total, we have constructed, for any D ≥ 38, sets of Kasner exponents {q i } i=1,⋯,D that satisfy(1.7a)-(1.7b) and (1.7d). That is, the Kasner solutions whose perturbations we study in our main theorem exist when D ≥ 38.

Norms and Bootstrap Assumptions
In this section, we define the various norms that we use to control the solution. We also state bootstrap assumptions for the solution norms; the bootstrap assumptions are convenient for our analysis in subsequent sections.
3.1. Constants featured in the norms. The norms that we define in Subsect. 3.4 involve the positive numbers q and σ featured in the following definition.
Definition 3.1 (The constants q and σ). Assuming that the Kasner exponents satisfy the condition (1.7d), we fix positive numbers q and σ verifying the following inequalities: Remark 3.1. We consider q and σ to be fixed for the remainder of the article.

Pointwise norms.
To control Σ t -tangent tensorfields, we will rely on two kinds of pointwise norms: one that refers to the transported spatial coordinate frame, and the standard geometric norm that is based on the Riemannian metric g.
3.3. Lebesgue and Sobolev norms. In this subsection, we define the Lebesgue and Sobolev norms that we will use to control the solution. We start by defining the ∂ ⃗ I -derivative of a tensorfield. is a type l m Σ t -tangent tensorfield and ⃗ I is a spatial multi-index, then we define ∂ ⃗ I T to be the tensorfield whose components (∂ ⃗ I T ) a 1 ⋯a l b 1 ⋯bm relative to the CMC-transported spatial coordinate frame are the following: Remark 3.2. The operator ∂ ⃗ I , when acting on Σ t -tangent tensorfields, can be given the following invariant interpretation: it can be viewed as repeated Lie differentiation with respect to the (globally defined) coordinate vectorfields In what follows, ⋅ L 2 (Σt) and ⋅ L ∞ (Σt) denote the standard Lebesgue norms for scalar functions on Σ t ; see Subsubsect. 1.8.4. We now define additional Sobolev and Lebesgue norms that we will use in our analysis.
Definition 3.4 (Sobolev and Lebesgue norms). If T is a type l m Σ t -tangent tensorfield, p ∈ {2, ∞}, and M ≥ 0 is an integer, then we define Remark 3.3 (The omission of norm subscripts when appropriate). If T is a scalar function, then we typically omit the subscripts "F rame" and "g" in the norms since there is no danger of confusion over how to measure the norm of T . For example if f is scalar function, then we write The specific norms that we use to control the solution. Let A ≫ 1 be a large parameter, to be chosen later. We recall that q > 0 and σ > 0 are the real numbers fixed in Subsect. 3.1.
To control the solution variables (g, k, n), we will rely on a combination of norms for the low-order derivatives of the solution and norms for its high-order derivatives. We now define the low-order norms.
Definition 3.5 (Low norms). Letg andk denote the background Kasner solution variables and let q and σ be the fixed constants (which depend on the background Kasner solution (g,k)) satisfying (3.1). We define We will use the following norms to control the high-order derivatives of the solution.
Definition 3.6 (High norms). Let q and σ be the fixed constants (which depend on the background Kasner solution (g,k)) that satisfy (3.1). Let A ≫ 1 be a real parameter (to be chosen later) and let N ≫ 1 be an integer-valued parameter (also to be chosen later). We define 3.5. Bootstrap assumptions. To facilitate our analysis, we find it convenient to rely on bootstrap assumptions. Let T (Boot) ∈ (0, 1) be a "bootstrap time". Until the proof of Theorem 11.1, we assume that the perturbed solution exists classically for (t, x) ∈ (T (Boot) , 1] × T D and that the following bootstrap assumptions hold for the norms from Defs. 3.5 and 3.6: where ε > 0 is a small bootstrap parameter.
Remark 3.4 (The required smallness of ε depends on various parameters). We will continually adjust the required smallness of ε throughout our analysis. The required smallness of ε is allowed to depend on N, A, q, σ, and D, but we will often avoid pointing this out.

Estimates for the Kasner Solution and Preliminary Technical Estimates
In this section, we derive some simple estimates for the background Kasner solution and provide other basic estimates that we will use throughout the paper.

Basic estimates for the Kasner solution.
In controlling various error terms, we will rely on the following simple estimates for the background Kasner solution.
Lemma 4.1 (Basic estimates for the Kasner solution). The following estimates hold for t ∈ (0, 1]: Proof. Recall that the Kasner solution variablesg andk are given, relative to the transported spatial coordinates, by the expressions (1.6) and (1.7c). All estimates stated in the lemma follow as straightforward consequences of these expressions and (3.1).

Norm comparisons.
In controlling various error terms, we will compare the pointwise norms ⋅ F rame and ⋅ g . Our comparisons will often rely on the following lemma.
Lemma 4.2 (Norm comparisons). Let T be a type l m Σ t -tangent tensorfield. Under the bootstrap assumptions (3.7), there exists a universal constant C * > 1 independent of N and A (but depending on m and l) such that if ε ≤ 1, then the following estimates hold for the pointwise norms of Def. 3.2 for t ∈ (T (Boot) , 1]: Proof. Let δ denote the standard Euclidean metric on Σ t , i.e., relative to the transported spatial coordinates, δ has components δ ij equal to diag(1, 1, ⋯, 1), and likewise for the inverse Euclidean metric δ −1 . Then in view of the definition of the norm ⋅ F rame , we have, schematically, T F rame = (δ) l (δ −1 ) m T 2 1 2 , and by g-Cauchy-Schwarz, the RHS of the previous expression is From Def. 3.5, the estimate (4.1a), and the bootstrap assumptions, we deduce that the RHS of the above expression (which we consider from the perspective of the transported coordinate frame) is which yields the second inequality in (4.2). To obtain the first inequality in (4.2), we note that we have, schematically, T g = (g) l (g −1 ) m T 2 1 2 . Writing out the RHS of the previous expression in the transported coordinate frame, we see that it is 3.5, the estimate (4.1a), and the bootstrap assumptions, we deduce that the RHS of the previous expression is ≤ C * t −(m+l)q T F rame , which yields the first inequality in (4.2).

Sobolev interpolation and product inequalities.
In this subsection, we provide some Sobolev interpolation and product inequalities that we will use to control various error terms. We start with the following lemma, which provides basic interpolation estimates.
If in addition M 1 + 1 + ⌊D 2⌋ ≤ M 2 , then the following inequalities hold: Proof. The first inequality in (4.3) follows as special case of Nirenberg's famous interpolation results [39], except that on the RHS, we have replaced the norm ⋅ L 2 (Σt) with ⋅ L ∞ (Σt) ; the replacement is possible because of the estimate v L 2 (Σt) ≲ v L ∞ (Σt) for scalar functions v (which holds because T D is compact). Strictly speaking, Nirenberg stated his results for functions defined on R D , but the arguments given in his paper can be used to derive the same estimates for functions defined on T D . The second inequality in (4.3) follows from the first and Young's inequality. The first inequality in (4.4) is a standard Sobolev embedding result on T D . The second inequality in (4.4) follows from the first and the second inequality in (4.3).
Lemma 4.4 (Sobolev product estimates for tensorfields). Let M be a non-negative integer, let {v r } R r=1 be a finite collection of Σ t -tangent tensorfields, and let v = ∏ R r=1 ∂ ⃗ Ir v r be a (schematically denoted) tensor product, possibly involving contractions. Assume that the bootstrap assumptions (3.7) hold for some ε with ε ≤ 1. Then the following estimate holds for t ∈ (T (Boot) , 1]: Moreover, under the same assumptions as above, and assuming that v is type l m , the following estimate holds for t ∈ (T (Boot) , 1]: Proof. If the {v r } R r=1 are all scalar functions (and hence l = m = 0), then inequality (4.5) is standard; it is proved, for example, as [47,Lemma 6.16]. If one or more of the v r are not scalar functions, then the estimate (4.5) follows a straightforward consequence of the estimate (4.5) for scalar functions, essentially by writing out the definition of LHS (4.5) relative to the transported spatial coordinate frame and estimating the components of all tensorfields.
We then use (4.5) to conclude that the RHS of the previous expression is ≲ RHS (4.6) as desired.
4.4. Sobolev embedding. We will use the following lemma to obtain L ∞ (Σ t ) control over some additional low-order derivatives that are not directly controlled by the low-order norms from Def. 3.5. To achieve the desired control, we borrow, via interpolation with sufficiently large N, a small amount of the high norm. Because the high norms are quite weak near t = 0 (at least when A is large), the interpolation introduces singular behavior into the estimates of the lemma, represented by the factors of t −Aδ on the RHSs of the estimates. However, δ → 0 as N → ∞. Thus, at fixed A, if N is large, then the following basic principle, central to our approach, applies: The singular contribution to the behavior of the low-order norms coming from the high norms is very small. There exists a parameter 13 δ > 0, depending on N and D, such that lim N →∞ δ = 0 and such that if N ≥ 5 + ⌊D 2⌋, then the following estimates hold for t ∈ (T (Boot) , 1]: Proof. To prove (4.7c), we first use (4.4) and the second inequality in (4.3) to deduce n−1 W 4,∞ (Σt) ≲ n − 1 L ∞ (Σt) + n − 1 Ḣ5+⌊D 2⌋ (Σt) . Then using the first inequality in (4.3), we deduce that for N ≥ 5 + ⌊D 2⌋, the RHS of the previous expression is ≲ n − 1 L ∞ (Σt) + n − 1 1−δ L ∞ (Σt) n δḢ N (Σt) , where δ ∶= (5 + ⌊D 2⌋) N. Combining these estimates and appealing to definitions (3.5b) and (3.6b), we find that Young's inequality to deduce where in obtaining the last inequality in (4.8), we have used (3.1), the assumption A ≥ 1, and our running convention that δ is free to vary from line to line, subject only to the restriction lim N →∞ δ = 0. In total, we have derived the desired bound (4.7c). The remaining estimates in the lemma can be proved using similar arguments that take into account Defs. 3.5 and 3.6, and we omit the details.

Control of n in Terms of g and k
Our primary goal in this section is to prove the next proposition, in which we derive the main estimates for n. The proof of the proposition is located in Subsect. 5.3. In Subsects. 5.1-5.2, we derive the identities and estimates that we will use when proving the proposition. The proposition shows that n is controlled, in various norms, by complementary norms of g and k. Achieving such control is possible since n solves an elliptic PDE with source terms that depend on R; here it makes sense to recall that, as we explained at the end of Subsect. 1.5, to obtain suitable estimates for n at the top derivative level, we use the Hamiltonian constraint (2.9a) to algebraically replace, in the elliptic lapse PDE, R with terms that depend on k.
Proposition 5.1 (Control of n in terms of g and k). Recall that L (g,k) (t), L (n) (t), H (g,k) (t), and H (n) (t) are the norms from Defs. 3.5 and 3.6, and assume that the bootstrap assumptions (3.7) hold. There exists a universal constant C * > 0 independent of N and A such that if N is sufficiently large in a manner that depends on A and if ε is sufficiently small, then the following estimates hold for t ∈ (T (Boot) , 1] (where, as we described in Subsect. 1.8.6, constants "C" are allowed to depend on N and other quantities).
Estimates at the lowest order. The following estimate holds: Estimate for ∂ t n. The following estimate holds: Top-order estimates. If ⃗ I = N, then the following estimates hold: Near-top-order estimates. The following estimate holds: 5.1. The equations. In this subsection, we derive the equations that we use to control n in terms of g and k.
5.1.1. The equation at the lowest order. We start with the following lemma, which we will use to derive estimates for n − 1 L ∞ (Σt) .
Lemma 5.2 (A rewriting of the lapse equation). The quantity n − 1 verifies the following elliptic PDE on Σ t : Proof. The lemma is a simple consequence of equation (2.11).

5.1.2.
The equations verified by the time derivative. The next lemma is an analog of Lemma 5.2 for ∂ t n, and we will use it to derive estimates for ∂ t n L ∞ (Σt) . We remark that we need estimates for ∂ t n L ∞ (Σt) only in Subsect. 10.2, when we bound the length of past-directed causal geodesic segments.
Lemma 5.3 (The ∂ t -commuted lapse equation). The quantity ∂ t n verifies the following elliptic PDE on Σ t : where (see Remark 2.1) Proof. Relative to CMC-transported spatial coordinates, the elliptic PDE (5.5) can be expressed as (2.12b) and recall that R = Ric a a ) and g ab Γ c a b ≃ g −2 ∂g. Commuting the elliptic PDE with ∂ t , using equations (2.10a) and (2.10b) to algebraically substitute for ∂ t g and ∂ t g −1 , and carrying out straightforward computations, we conclude the desired equation (5.6). 5.1.3. The equations verified by the high-order derivatives. We will use the equations in the next lemma to control the L 2 (Σ t ) norms of high-order derivatives of n.
Remark 5.1 (Borderline error terms vs. junk error terms). In the rest of the paper, we will denote difficult "borderline" error terms by decorating them with "(Border)", e.g., (Border; ⃗ I) N. These error terms must treated with care since at the top derivative level, there is no room (in the sense of powers of t) in our estimates for such terms. In contrast, the error terms that we decorate with "(Junk)", such as (Junk; ⃗ I) N, are easier to treat in the sense that there is some room in our estimates.

Lemma 5.4 (The ⃗
I-commuted lapse equation). For each spatial multi-index ⃗ I with ⃗ I = N, ∂ ⃗ I n verifies the following equation: where (see Subsubsect. 1.8.7 regarding our use of the notation * ≃ and ≃, and see Remark 2.1) Furthermore, for each spatial multi-index ⃗ I with ⃗ I = N − 1, ∂ ⃗ I n verifies the following equation: Proof. First, we use equation (2.9a) to substitute for R in equation (2.11), use that g ab ∇ a ∇ b n = g ab ∂ a ∂ b n − g ab Γ c a b ∂ c n, and multiply by t A+2 , thereby deriving the equation 12) Commuting equation (5.12) with ∂ ⃗ I and using the schematic identity g ab Γ c a b ∂ c n ≃ g −2 (∂g)∂n, we easily deduce (5.8). We clarify that the terms on RHS (5.9a) arise, respectively, when all ⃗ I derivatives fall on the factor n − 1 in the first product on RHS (5.12) or on one of the factors of k in the next-to-last product t A+2 k a b k b a on RHS (5.12). It is for this reason that the coefficients in the corresponding products do not depend on N, as is indicated by the symbol * ≃ in equation (5.9a). Equation (5.10) follows similarly from multiplying equation (5.5) by t A+1+q and using the schematic identity R ≃ (g −1 ) 2 ∂ 2 g + (g −1 ) 3 (∂g) 2 , which follows from (2.12b) and the fact that R = Ric a a .

5.2.
Control of error the terms in the lapse estimates. In this subsection, we derive estimates for the error terms in the equations of Lemmas 5.3 and 5.4.
Lemma 5.5 (Control of the error terms in the elliptic estimates). Assume that the bootstrap assumptions (3.7) hold. There exists a universal constant C * > 0 independent of N and A such that if N is sufficiently large in a manner that depends on A and if ε is sufficiently small, then the following estimates hold for t ∈ (T (Boot) , 1] (where, as we described in Subsect. 1.8.6, constants "C" are allowed to depend on N and other quantities).
Control of the error term in the equation verified by ∂ t n. The following estimate holds for the error term N ′ defined in (5.7): Borderline top-order error term estimates. If ⃗ I = N, then the following estimate holds for the error term from (5.9a): Non-borderline top-order error term estimates. The following estimate holds for the error term from (5.9b): Just-below-top-order error term estimates. The following estimate holds for the error term from (5.11): Proof. Throughout this proof, we will assume that Aδ is sufficiently small (and in particular that Aδ < σ); in view of the discussion in Subsect. 4.4, we see that at fixed A, this can be achieved by choosing N to be sufficiently large.
Proof of (5.13): This estimate follows in a straightforward fashion from using (3.1), Def. 3.5, Lemma 4.1, and Lemma 4.5 to control the products on RHS (5.7) in the norm ⋅ L ∞ (Σt) , where we bound each factor in the products in the norm ⋅ L ∞ F rame (Σt) . We remark that the most singular (in the sense of powers of t) products on RHS (5.7) are t −3 (n − 1) and the terms in the second sum with p = i 1 = 0.
Proof of (5.15): Let ⃗ I be a spatial multi-index with ⃗ I = N. To bound the first sum on RHS (5.9b), we first consider the top-order case in which ⃗ I 2 = N or ⃗ I 3 = N. Then using g-Cauchy-Schwarz, we see that the terms under consideration are bounded in the norm ⋅ L 2 (Σt) by . From Defs. 3.5 and 3.6 and the bootstrap assumptions, we see that RHS of the previous expression is ≲ t 2−10q−σ H (g,k) (t). In view of (3.1), we see that RHS of the previous expression is ≲ RHS (5.15) as desired. We now consider the remaining cases, in which ⃗ I 1 , ⃗ I 2 , ⃗ I 3 ≤ N − 1. Using (4.3) and (4.5), we bound (using that ⃗ I = N) the terms under consideration as follows: . From Defs. 3.5 and 3.6, the estimates (4.1b) and (4.7c), and the bootstrap assumptions, we deduce that RHS (5.17) ≲ t 3−13q−2σ−Aδ H (g,k) (t) + t 1−q L (g,k) (t). In view of (3.1), we see that the RHS of the previous expression is ≲ RHS (5.15) as desired.
To bound the second sum on RHS (5.9b), we first use (4.5) to deduce that the terms under consideration are bounded in the norm . From Defs. 3.5 and 3.6 and the bootstrap assumptions, we see that the RHS of the previous expression is ≲ t 1−3q−σ H (g,k) (t), which, in view of (3.1), is ≲ RHS (5.15) as desired.
To bound the third sum on RHS (5.9b), we first consider the top-order case in which ⃗ I 3 = N. Using g-Cauchy-Schwarz and the fact that g −1 g ≲ 1, we see that the terms under consideration are bounded in the norm ⋅ L 2 (Σt) by ≲ t A+2 ∂n L ∞ g (Σt) ∂g ḢN g (Σt) . Applying (4.2) (with l = 0 and m = 1) to ∂n L ∞ g (Σt) , we see that the RHS of the previous expression is ≲ t A+2−q n Ẇ 1,∞ (Σt) ∂g ḢN g (Σt) . From Defs. 3.5 and 3.6, the estimate (4.7c), and the bootstrap assumptions, we see that the RHS of the previous expression is ≲ t 3−11q−σ−Aδ H (g,k) (t), which, in view of (3.1), is ≲ RHS (5.15) as desired. We now consider the top-order case in which ⃗ I 4 = N. Arguing as above, we deduce that the terms under consideration are bounded in the norm . From Defs. 3.5 and 3.6, the estimate (4.7a), and the bootstrap assumptions, we see that the RHS of the previous expression is ≲ t 1−5q−Aδ L (g,k) (t)H (n) (t) ≲ t 1−5q−Aδ L (g,k) (t), which, in view of (3.1), is ≲ RHS (5.15) as desired. It remains for us to handle the cases in which ⃗ I 3 , ⃗ I 4 ≤ N − 1. Using (4.3) and (4.5), we bound (using that ⃗ I = N) the terms under consideration as follows: From Defs. 3.5 and 3.6, the estimates (4.1a) and (4.7a)-(4.7c), and the bootstrap assumptions, we deduce that RHS (5.18) ≲ t 4−16q−2σ−Aδ L (g,k) (t) + H (g,k) (t) + t 2−6q−Aδ L (g,k) (t) + H (g,k) (t) . In view of (3.1), we see that the RHS of the previous expression is ≲ RHS (5.15) as desired.
To bound the last sum on RHS (5.9b), we first consider the case in which ⃗ I 2 = N −1. Using (4.2), we bound the terms under consideration in the norm . From Defs. 3.5 and 3.6, the estimate (4.7b), and the bootstrap assumptions, we see that the RHS of the previous expression is ≲ t 1−3q−Aδ L (g,k) (t) + H (g,k) (t) H (n) (t) ≲ t 1−3q−Aδ L (g,k) (t) + H (g,k) (t) , which, in view of (3.1), is ≲ RHS (5.15) as desired. We now consider the remaining cases, in which ⃗ I 2 ≤ N − 2. Using (4.5), we bound (using that ⃗ I = N) the terms under consideration as follows:

ḢN
F rame (Σt) . From Def. 3.6, the estimates (4.7b) and (4.7c), and the bootstrap assumptions, we deduce that In view of (3.1), we see that the RHS of the previous expression is ≲ RHS (5.15) as desired.
Proof of (5.16): Let ⃗ I be a spatial multi-index with ⃗ I = N − 1. To bound the first sum on RHS (5.11), we first consider the case in which ⃗ I 3 = N − 1. Using g-Cauchy-Schwarz and the fact that g −1 g ≲ 1, we see that the terms under consideration are bounded in the norm ⋅ L 2 (Σt) by ≲ t A+1+q ∂ 2 g ḢN−1 g (Σt) . From the definitions of ⋅ ḢM g (Σt) and ⋅ L ∞ F rame , we deduce that the RHS of the previous estimate is ≲ t A+1+q g −1 1 2 L ∞ F rame ∂g ḢN g (Σt) . From Defs. 3.5 and 3.6, the estimate (4.1a), and the bootstrap assumptions, we see that the RHS of the previous expression is ≲ H (g,k) (t), which is ≲ RHS (5.16) as desired. It remains for us to consider the remaining cases, in which ⃗ I 3 ≤ N − 2. Using (4.3) and (4.5), we bound (using that ⃗ I = N − 1) the terms under consideration as follows: . From Defs. 3.5 and 3.6, the estimates (4.1a), (4.7a), and (4.7b), and the bootstrap assumptions, we see that RHS (5.20) ≲ t 1−5q−σ−Aδ L (g,k) (t) + H (g,k) (t) , which, in view of (3.1), is ≲ RHS (5.16) as desired.
To bound the last sum on RHS (5.11), we first use (4.5) and the fact that ⃗ I = N − 1 to deduce that the terms under consideration are bounded as follows: From Defs. 3.5 and 3.6, the estimates (4.7b) and (4.7c), and the bootstrap assumptions, we see that RHS (5.25) ≲ t 1−q−Aδ L (g,k) (t) + H (g,k) (t) + t 4−14q−4σ−Aδ L (g,k) (t) + H (g,k) (t) , which, in view of (3.1), is ≲ RHS (5.16) as desired. This completes the proof of (5.16) and finishes the proof of the lemma. 5.3. Proof of Prop. 5.1. In this subsection, we prove Prop. 5.1. Throughout this proof, we will assume that Aδ is sufficiently small (and in particular that Aδ < σ); in view of the discussion in Subsect. 4.4, we see that at fixed A, this can be achieved by choosing N to be sufficiently large.
Proof of (5.3) and (5.4): In view of Def. 3.6, we see that to obtain (5.3) and (5.4), it suffices to prove the following estimates: To prove (5.28), we let ⃗ I be any spatial derivative multi-index with ⃗ I = N. Multiplying equation (5.8) by t A ∂ ⃗ I n and integrating by parts over Σ t , we deduce Next, we use (3.1), (4.1a), (4.7b), the bootstrap assumptions, and g-Cauchy-Schwarz to deduce that (5.31), and Young's inequality, we deduce that if ε is sufficiently small, then Using (5.14) and (5.15) to bound the last two integrals on RHS (5.32), and soaking (assuming ε is sufficiently small) the first two terms on RHS (5.32) and the first term on RHS (5.14) back into LHS (5.32), we arrive at (5.28). Similarly, to prove (5.29), we let ⃗ I be any spatial derivative multi-index with ⃗ I = N − 1. We multiply equation (5.10) by t A+q−1 ∂ ⃗ I n, use (5.26), use the bound (t∂ b g ab )(t A+q ∂ a ∂∂ ⃗ I n)(t A+q−1 ∂ ⃗ I n) ≲ ε t A+q ∂∂ ⃗ I n g t A+q−1 ∂ ⃗ I n (which follows from essentially the same reasoning we used to prove (5.31)), and argue as in the previous paragraph to deduce the following analog of (5.32): Soaking the first two terms on RHS (5.33) back into LHS and using (5.16) to bound the last integral on RHS (5.33), we arrive at (5.29).

Estimates for the Low-Order Derivatives of g and k
Our main goal in this section is to prove the following proposition, which provides the integral inequality that we use to control the low norm L (g,k) (t). The proof of the proposition is located in Subsect. 6.2. In Subsect. 6.1, we derive the equations that we will use in proving it. Proposition 6.1 (Integral inequality for L (g,k) (t)). Recall that L (g,k) (t) is the low-order norm from Def. 3.5, and assume that the bootstrap assumptions (3.7) hold. There exists a constant C > 0 such that if N is sufficiently large in a manner that depends on A and if ε is sufficiently small, then the following estimate holds for t ∈ (T (Boot) , 1]: L (g,k) (t) ≤ L (g,k) (1) + C 1 s=t s σ−1 L (g,k) (s) + H (g,k) (s) ds.
(6.1) 6.1. The equations. In this subsection, we derive the equations that we use to control g and k at the lowest derivative levels.
Lemma 6.2 (A rewriting of the equations verified by g and k). Letg andk denote the background Kasner solution variables with corresponding Kasner exponents {q i } i=1,⋯,D . Then the following evolution equations hold, where i ≤ j in (6.2a)-(6.2b) and there is no summation over j in (6.2a)-(6.2b) (and note that t −2q jg ij = t 2q j (g −1 ) ij = diag(1, 1, ⋯, 1)): Moreover, the following evolution equation holds (and note that tk i j = −diag(q 1 , ⋯, q D )): Proof. Throughout this proof, there is no Einstein summation over j. To derive equation (6.2a), we first use equation (2.10a) and the fact that t −2q jg ij = δ ij (where δ ij is the standard Kronecker delta) to deduce Sincek = −t −1 diag(q 1 , ⋯, q D ), we can express the first two products on RHS (6.4) as follows: −2t −2q j g ia k a j − 2q j t −2q j −1 g ij = −2t −2q j g ia k a j −k a j . Next, using thatg = diag(t 2q 1 , ⋯, t 2q D ) andk = −t −1 diag(q 1 , ⋯, q D ), we express the RHS of the previous expression as follows: Combining these calculations, we arrive at (6.2a). Equation (6.2b) can be derived by applying similar arguments to equation (2.10b), and we omit the details. Equation (6.3) follows from multiplying both sides of equation (2.10c) by t and using that tk = −diag(q 1 , ⋯, q D ).

6.2.
Proof of Prop. 6.1. In this subsection, we prove Prop. 6.1. In this proof, we will assume that Aδ is sufficiently small (and in particular that Aδ < σ); in view of the discussion in Subsect. 4.4, we see that at fixed A, this can be achieved by choosing N to be sufficiently large.
First, we note that to obtain (6.1), it suffices to show that the following bounds hold: For we can then use the symmetry property k ab = k ba and the fact that (tk a b )(tk b a ) = 1 to derive the identity tk with the bootstrap assumptions (3.7) and the estimates (4.1b) and (6.7), also yields (upon Taylor expanding the square root) the estimate In view of definition (3.5a), from (6.5)-(6.8), we conclude the desired estimate (6.1). It remains for us to prove (6.5)-(6.7). We first prove (6.5) by analyzing equation (6.2a). From Def. 3.5, the estimates (4.1a), (4.1b), and (5.1), and the bootstrap assumptions, we deduce, by bounding each factor on RHS (6.2a) in the norm ⋅ L ∞ F rame (Σt) , that the following estimate holds: From (3.1), we deduce that the last two products on RHS (6.9) are ≤ Ct σ−1−2q−2q j L (g,k) (t) + H (g,k) (t) for t ∈ (T (Boot) , 1]. Hence, with the help of these bounds, we can integrate equation (6.2a) from time t to time 1 and use the initial data bound g −g L ∞ F rame (Σ 1 ) ≤ CL (g,k) (1) to deduce the following estimate for components: Multiplying both sides of (6.10) by t 2q+2q j and using (3.1), we deduce L ∞ F rame (Σs) ds (6.11) We now define G(t) ∶= sup s∈[t,1] s 2q g − s 2qg L ∞ F rame (Σs) . Next, with the help of (3.1), we deduce the following estimate for the first integral on RHS (6.11): From (6.11) and (6.12), we deduce that for (t, x) ∈ (T (Boot) , 1] × T D , we have From (6.13), it follows that For ε sufficiently small, we can absorb the product CεG(t) on RHS (6.14) back into the LHS (at the minor expense of increasing the constants C on the RHS). From these arguments, we conclude the bound G(t) ≤ CL (g,k) (1) + C ∫ The estimate (6.6) can be proved using a similar argument based on the evolution equation (6.2b), and we omit these details.

Estimates for the Top-Order Derivatives of g and k
Our main goal in this section is to prove the following proposition, which provides the main integral inequality for the top-order derivatives of g and k. The proof is located in Subsect. 7.4. In Subsects. 7.1-7.3, we derive the identities and estimates that we will use when proving the proposition.
Proposition 7.1 (Integral inequality for the top-order derivatives of g and k). Let ⃗ I be top-order spatial multi-index, that is, a multi-index with ⃗ I = N. Assume that the bootstrap assumptions (3.7) hold. There exists a universal constant C * > 0 independent of N and A such that if N is sufficiently large in a manner that depends on A and if ε is sufficiently small, then the following integral inequality holds for t ∈ (T (Boot) , 1] (where, as we described in Subsect. 1.8.6, constants "C" are allowed to depend on N and other quantities): 7.1. The equations. In this subsection, we derive the equations that we will use when deriving estimates for g and k at the high derivative levels.
Lemma 7.2 (The equations verified by the high-order derivatives of the metric and second fundamental form). Let ⃗ I be top-order spatial multi-index, that is, a multi-index with ⃗ I = N. Then the following commuted momentum constraint equations hold: where (see Subsubsect. 1.8.7 regarding our use of the notation * ≃ and ≃, and see Remark 2.1) Moreover, for any constant P ≥ 0 and for any spatial multi-index ⃗ I, the following commuted evolution equations hold: Proof. To prove (7.2a), we first write equation (2.9b) relative to the transported spatial coordinates in the schematic form ∂ a k a i = g −1 (∂g)k. Commuting this equation with t A+1 ∂ ⃗ I , we arrive at (7.2a). The proof of (7.2b) is similar, but we start by raising the index i in equation (2.9b) to deduce g ab ∇ a k i b = 0 and then (schematically) writing this equation in coordinates as g ab ∂ a k i b = g −2 (∂g)k. To prove (7.4a), we commute equation (2.10a) with ∂ e ∂ ⃗ I and then with t A+P and carry out straightforward computations.
To prove (7.4b), we first use (2.12b) to decompose where (schematically) Error = g −3 (∂g) 2 . We now use (7.7) to decompose the term Ric i j on RHS (2.10c), commute the evolution equation (2.10c) with ∂ ⃗ I and then with t A+P , and carry out straightforward computations, thereby arriving at (7.4b) .
7.2. Energy currents. When deriving top-order energy estimates (see the proof of Prop. 7.1), we will integrate by parts by applying the divergence theorem on spacetime regions of the form  To each top-order spatial multi-index ⃗ I (i.e., ⃗ I = N), we associate the energy current ( ⃗ I) J, which we define to be the vectorfield with the following components relative to the CMC-transported spatial coordinates: Remark 7.1. Note that the components { ( ⃗ I) J α } α=0,⋯,D are quadratic forms in (∂ ⃗ I k, ∂∂ ⃗ I g, ∂∂ ⃗ I n) with coefficients that depend on g, g −1 , and n. Let ⃗ I be a top-order spatial multi-index, that is, a multi-index with ⃗ I = N. Then for solutions to the commuted equations of Lemma 7.2 with P = 1 in equations (7.4a) and (7.4b), the spacetime vectorfield ( ⃗ I) J from Def. 7.1 verifies the following divergence identity relative to the CMC-transported spatial coordinates: Proof. We view RHSs (7.8a)-(7.8b) as quadratic forms in (∂ ⃗ I k, ∂∂ ⃗ I g, ∂∂ ⃗ I n) with coefficients that depend on g, g −1 , and n. We now consider the expression ∂ α ( ⃗ I) J α . On RHS (7.10b), we place all terms in which spatial derivatives ∂ j fall on the coefficients. In contrast, when ∂ t falls on g or g −1 , we use (2.10a)-(2.10b) to substitute for ∂ t g and ∂ t g −1 and then place the resulting terms on RHS (7.10a) (as the last four products). Next, we consider all of the terms in the expression ∂ α ( ⃗ I) J α in which a derivative falls on one of ∂ ⃗ I k, ∂∂ ⃗ I g, or ∂∂ ⃗ I n. For these terms, we use equations (7.2a)-(7.2b) and (7.4a)-(7.4b) for algebraic substitution, where P = 1 (by assumption) in (7.4a)-(7.4b). More precisely, in the expression ∂ α ( ⃗ I) J α , we use equation (7.4b) to substitute for the factor and equation (7.2b) to substitute for the factor g jc These algebraic substitutions lead to the elimination of all products that depend on a derivative of (∂ ⃗ I k, ∂∂ ⃗ I g, ∂∂ ⃗ I n), at the expense of introducing products that depend on the inhomogeneous terms, for example 2g ac g bd (t A+1 ∂ ⃗ I k a b ) (Border;1; ⃗ I) K c d and 2g ac g bd (t A+1 ∂ ⃗ I k a b ) (Junk;1; ⃗ I) K c d . We place the borderline error term products such as 2g ac g bd (t A+1 ∂ ⃗ I k a b ) (Border;1; ⃗ I) K c d on RHS (7.10a) and the junk error term products such as 2g ac g bd (t A+1 ∂ ⃗ I k a b ) (Junk;1; ⃗ I) K c d on RHS (7.10b). Moreover, as the first product on RHS (7.9), we place 2A t t A+1 ∂ ⃗ I k 2 g , which is generated by the first term on RHS (7.4b) when we use equation (7.4b) (with P = 1) to substitute for the factor ∂ t (t A+1 ∂ ⃗ Finally, as the second product on RHS (7.9), we place A+1 2t t A+1 ∂∂ ⃗ I g 2 g , which is generated by the first term on RHS (7.4a) when we use equation (7.4a) to substitute for the factor ∂ t (t A+1 ∂ ⃗ I ∂ d ∂g ef ) in the expression 1 2 g ad g be g cf (t A+1 ∂ ⃗ I ∂ a ∂g bc )∂ t (t A+1 ∂ ⃗ I ∂ d ∂g ef ). In total, these steps yield the lemma. 7.3. Control of the error terms. In this subsection, we derive estimates for the error terms that will arise when we derive energy estimates for solutions to the equations of Lemma 7.1.

7.3.1.
Pointwise estimates for the error terms in the divergence of the energy current. In this subsubsection, we derive pointwise estimates for the error terms (Border; ⃗ I) J and (Junk; ⃗ I) J that appear in the expression (7.9) for ∂ α Lemma 7.4 (Pointwise estimates for the error terms in the divergence of the energy current). Let ⃗ I be a top-order spatial multi-index, that is, a multi-index with ⃗ I = N. Assume that the bootstrap assumptions (3.7) hold. There exists a universal constant C * > 0 independent of N and A such that if N is sufficiently large in a manner that depends on A and if ε is sufficiently small, then the error terms (Border; ⃗ I) J and (Junk; ⃗ I) J from (7.10a) and (7.10b) verify the following pointwise estimates for t ∈ (T (Boot) , 1] (where, as we described in Subsect. 1.8.6, constants "C" are allowed to depend on N and other quantities): 11a) + C * t (Border;1; ⃗ I) K Proof. Throughout this proof, we will assume that Aδ is sufficiently small (and in particular that Aδ < σ); in view of the discussion in Subsect. 4.4, we see that at fixed A, this can be achieved by choosing N to be sufficiently large. The estimate (7.11a) follows in a straightforward fashion from applying the g-Cauchy-Schwarz inequality and Young's inequality (in the form ab ≤ 1 2 (a 2 + b 2 )) to the products on RHS (7.10a) and using the estimate g g = g −1 g ≤ C * as well as the pointwise estimates n ≤ C * and k g ≤ C * t −1 , which follow from (3.1), Def. 3.5, and the bootstrap assumptions.
7.3.2. L 2 estimates for the error terms in the divergence of the energy current. In this subsubsection, we bound the error terms from Lemma 7.2 in the norm ⋅ L 2 g (Σt) . Lemma 7.5 (L 2 control of the error terms in the top-order energy estimates for g and k). Assume that the bootstrap assumptions (3.7) hold. There exists a universal constant C * > 0 independent of N and A such that if N is sufficiently large and if ε is sufficiently small, then the following estimates hold for t ∈ (T (Boot) , 1] (where, as we described in Subsect. 1.8.6, constants "C" are allowed to depend on N and other quantities).
Borderline top-order error term estimates. For each top-order spatial multi-index ⃗ I (i.e., ⃗ I = N), the following estimates hold for the error terms from (7.3a), (7.3c), (7.5a), and (7.6a), where P = 1 on RHSs (7.5a), and (7.6a): Non-borderline top-order error term estimates. The following estimates hold for the error terms from (7.3b), (7.3d), (7.5b), and (7.6b), where P = 1 on RHSs (7.5b) and (7.6b): Proof. Throughout this proof, we will assume that Aδ is sufficiently small (and in particular that Aδ < σ); in view of the discussion in Subsect. 4.4, we see that at fixed A, this can be achieved by choosing N to be sufficiently large.
Proof of (7.12a)-(7.12d): Let ⃗ I be a spatial multi-index with ⃗ I = N. We first use Def. 3.5, the fact that g −1 g ≤ C * , the bootstrap assumptions, and g-Cauchy-Schwarz to bound the product on RHS (7.3a) as follows: (Border; ⃗ Taking the norm ⋅ L 2 (Σt) of this inequality, we conclude that (Border; ⃗ I) M ≤ C * t A ∂∂ ⃗ I g L 2 g (Σt) as desired. The estimate (7.12b) follows similarly based on equation (7.3c), and we omit the details. The estimates (7.12c) and (7.12d) follow similarly based on equations (7.5a) and (7.6a) (with P = 1 by assumption) and the elliptic estimate (5.3), and we omit the details.
Proof of (7.13a)-(7.13b): We prove only (7.13b); the estimate (7.13a) can be proved by applying a similar argument to the products on RHS (7.3b) and we omit those details. Let ⃗ I be a spatial multi-index with ⃗ I = N. To obtain the desired bound for the sum on RHS (7.3d), we first consider the cases in which ⃗ I 1 = N or ⃗ I 2 = N. Using that g −1 g ≲ 1 and g-Cauchy-Schwarz, we first bound the products under consideration in the norm = 0 and m = 3), we deduce that the RHS of the previous expression is 3.5 and 3.6, the estimate (4.7a), and the bootstrap assumptions, we deduce that the RHS of the previous expression is ≲ t −6q−σ−Aδ H (g,k) (t), which, in view of (3.1), is ≲ RHS (7.13b) as desired. We now consider the case in which ⃗ I 3 = N − 1 on RHS (7.3d). Using that g −1 g ≲ 1 and g-Cauchy-Schwarz, we first bound the products under consideration in the norm From Defs. 3.5 and 3.6, the estimate (4.7b), and the bootstrap assumptions, we deduce that the RHS of the previous expression is ≲ t −6q−σ−Aδ H (g,k) (t), which, in view of (3.1), is ≲ RHS (7.13b) as desired. We now consider the case in which ⃗ I 4 = N on RHS (7.3d). Using that g −1 g ≲ 1 and g-Cauchy-Schwarz, we first bound the products under consideration in the norm ⋅ L 2 g (Σt) by ≲ t A+1 ∂g L ∞ g k ḢN g . Applying (4.2) to ∂g L ∞ g (with l = 0 and m = 3), we deduce that the RHS of the previous expression is ≲ t A+1−3q g Ẇ 1,∞ F rame k ḢN g . From Def. 3.6, the estimate (4.7a), and the bootstrap assumptions, we deduce that the RHS of the previous expression is ≲ t −5q−Aδ H (g,k) (t), which, in view of (3.1), is ≲ RHS (7.13b) as desired. It remains for us to consider the remaining terms on RHS (7.3d), in which ⃗ I 1 , ⃗ I 2 , ⃗ I 4 ≤ N − 1 and ⃗ I 3 ≤ N − 2. We first use (4.2) with l = 1 and m = 0 (since (Junk; ⃗ I) M i is a type 1 0 tensorfield), (4.3), and (4.5) to bound (using that ⃗ I = N) the terms under consideration as follows: . From Defs. 3.5 and 3.6, the estimates (4.1a), (4.1b), (4.7a), and (4.7b), and the bootstrap assumptions, we deduce that RHS (7.14) ≲ t 1−10q−3σ−Aδ L (g,k) (t) + H (g,k) (t) . In view of (3.1), we see that the RHS of the previous expression is ≲ RHS (7.13b) as desired.
Proof of (7.13c): We stress that for this estimate, on RHS (7.5b), we have P = 1 and ⃗ I = N.
To bound the first product on RHS (7.5b), we first use g-Cauchy-Schwarz to deduce that it is bounded in the norm ⋅ L 2 g (Σt) by ≤ t A+1 n − 1 L ∞ (Σt) k L ∞ g (Σt) ∂g ḢN g (Σt) . From Defs. 3.5 and 3.6, (3.1), and the bootstrap assumptions, we see that the RHS of the previous expression To complete the proof of (7.13c), we must bound the three sums on RHS (7.5b) in the norm ⋅ L 2 g (Σt) . To bound the first sum on RHS (7.5b), we first consider the products with ⃗ I 3 = N. Using that g g ≲ 1 and g-Cauchy-Schwarz, we first bound the products under consideration in the norm . Applying (4.2) to ∂n L ∞ g (Σt) (with l = 0 and m = 1), we deduce that the RHS of the previous expression is ≲ t A+1−q n Ẇ 1,∞ (Σt) k ḢN g (Σt) . From Defs. 3.6, the estimate (4.7c), and the bootstrap assumptions, we deduce that the RHS of the previous expression is ≲ t 2−11q−σ−Aδ H (g,k) (t), which, in view of (3.1), is ≲ RHS (7.13c) as desired. It remains for us to bound the products in the first sum on RHS (7.5b) with ⃗ I 1 , ⃗ I 3 ≤ N − 1. We first use (4.6) with l = 0 and m = 3 (since (Junk;P ; ⃗ I) H is a type 0 3 tensorfield) and (4.3) to deduce (using that ⃗ I = N) that the products under consideration are bounded as follows: From Defs. 3.5 and 3.6, the estimates (4.1b), (4.7a), and (4.7c), the elliptic estimates (5.1) and (5.3), and the bootstrap assumptions, we deduce that In view of (3.1), we see that RHS (7.16) ≲ RHS (7.13c) as desired.
To bound the second sum on RHS (7.5b), we first consider the products with ⃗ I 3 = N. Using g-Cauchy-Schwarz, we first bound the products under consideration in the norm . Applying (4.2) to ∂g L ∞ g (Σt) (with l = 0 and m = 3), we deduce that the RHS of the previous expression is ≲ t A+1−3q n L ∞ (Σt) g Ẇ 1,∞ F rame (Σt) k ḢN g (Σt) . From (3.1), Defs. 3.5 and 3.6, the estimate (4.7a), and the bootstrap assumptions, we deduce that the RHS of the previous expression is ≲ t −5q−Aδ H (g,k) (t), which, in view of (3.1), is ≲ RHS (7.13c) as desired. It remains for us to consider the remaining terms in the second sum on RHS (7.5b), in which ⃗ I 2 , ⃗ I 3 ≤ N − 1. We first use (4.6) with l = 0 and m = 3 (since (Junk;P ; ⃗ I) H is a type 0 3 tensorfield) and (4.3) to deduce (using that ⃗ I = N) that the terms under consideration are bounded as follows: . From (3.1), Defs. 3.5 and 3.6, the estimates (4.1b), (4.7a), and (4.7c), and the bootstrap assumptions, we deduce that RHS (7.17 In view of (3.1), we see that the RHS of the previous expression is ≲ RHS (7.13c) as desired.
Proof of (7.13d): We stress that for this estimate, on RHS (7.6b), we have P = 1 and ⃗ I = N. To bound the first sum on RHS (7.6b), we first consider the case in which ⃗ I 2 = N. Then the terms under consideration are bounded in the norm ⋅ L 2 g (Σt) by ≲ t A n−1 L ∞ (Σt) k ḢN g (Σt) . From Defs. 3.5 and 3.6 and (3.1), we see that the RHS of the previous expression is ≲ t 1−10q−σ H (g,k) (t) ≲ t σ−1 H (g,k) (t) as desired. We now consider the remaining cases, in which ⃗ I 1 ≤ N − 1 and ⃗ I 2 ≤ N − 1. First, using (4.6) with l = m = 1 (since (Junk;1; ⃗ I) K is a type 1 1 tensorfield) we bound (using that ⃗ I = N) the products under consideration as follows: . From Defs. 3.5 and 3.6, the estimate (4.7c), and the bootstrap assumptions, we deduce that RHS (7.19) ≲ t 2−15q−2σ−Aδ H (g,k) (t)+t −3q L (g,k) (t) + H (g,k) (t) . In view of (3.1), we see that the RHS of the previous expression is ≲ RHS (7.13d) as desired.
To bound the last sum on RHS (7.6b), we first consider the case in which ⃗ I 3 = N. Using that g −1 g ≲ 1 and g-Cauchy-Schwarz, we deduce that the terms under consideration are bounded in the norm ⋅ L 2 g (Σt) by ≲ t A+1 ∂n L ∞ g (Σt) ∂g ḢN g (Σt) . Next, using (4.2) (with l = 0 and m = 1) to estimate ∂n L ∞ g (Σt) , we see that the RHS of the previous expression is ≲ t A+1−q n Ẇ 1,∞ (Σt) ∂g ḢN g (Σt) . From Defs. 3.5 and 3.6, the estimate (4.7c), and the bootstrap assumptions, we see that the RHS of the previous expression is ≲ t 2−11q−σ−Aδ H (g,k) (t), which, in view of (3.1), is ≲ RHS (7.13d) as desired. We now consider the case in which ⃗ I 4 = N. Using that g −1 g ≲ 1 and g-Cauchy-Schwarz, we deduce that the terms under consideration are bounded in the norm ⋅ L 2 g (Σt) by ≲ t A+1 ∂g L ∞ g (Σt) ∂n ḢN g (Σt) . Next, using (4.2) (with l = 0 and m = 3) to estimate ∂g L ∞ g (Σt) , we see that the RHS of the previous expression is ≲ t A+1−3q g Ẇ 1,∞ F rame (Σt) ∂n ḢN g (Σt) . From Defs. 3.5 and 3.6 and the estimate (4.7a), we see that the RHS of the previous expression is ≲ t −5q−Aδ L (g,k) (t)H (n) (t) ≲ t −5q−Aδ L (g,k) (t). In view of (3.1), we see that the RHS of the previous expression is ≲ RHS (7.13d) as desired. We now consider the remaining cases, in which ⃗ I 3 ≤ N −1 and ⃗ I 4 ≤ N −1. Using (4.6) with l = m = 1 (since (Junk;1; ⃗ I) K is a type 1 1 tensorfield) and (4.3), we deduce (using that ⃗ I = N) that . From (3.1), Defs. 3.5 and 3.6, the estimates (4.1a), (4.7a), (4.7b), and (4.7c), and the bootstrap assumptions, we deduce that RHS (7.22) ≲ t 1−8q−σ−Aδ L (g,k) (t) + H (g,k) (t) . In view of (3.1), we see that the RHS of the previous expression is ≲ RHS (7.13d) as desired. We have therefore proved (7.13d), which completes the proof of the lemma. 7.4. Proof of Prop. 7.1. In this subsection, we prove Prop. 7.1. Throughout this proof, we will assume that Aδ is sufficiently small (and in particular that Aδ < σ); in view of the discussion in Subsect. 4.4, we see that at fixed A, this can be achieved by choosing N to be sufficiently large.
Let ⃗ I be a spatial multi-index with ⃗ I = N and let ( ⃗ I) J be the energy current from Def. 7.1. Applying the divergence theorem on the region [t, 1] × T D , and considering equation (7.8a), we deduce that for t ∈ (T (Boot) , 1], we have We now use equation (7.9) to substitute for the last integral on RHS (7.23), thereby obtaining ds Next, we use the estimates (7.12a)-(7.13d) to bound the terms on lines two to five of RHS (7.25). Also noting that the term ∫ Σ 1 ∂ ⃗ I k 2 g + 1 4 ∂∂ ⃗ I g 2 g dx on RHS (7.24) is ≤ CH 2 (g,k) (1), we arrive at the desired bound (7.1).

Estimates for the Near-Top-Order Derivatives of g and k
Our primary goal in this section is to prove the following proposition, which provides the main integral inequalities that we will use to control the near-top-order derivatives of g and k; recall that, as we explained near the end of Subsect. 1.5, for technical reasons, we need these estimates to close our bootstrap argument. The proof of the proposition is located in Subsect. 8.3. In Subsects. 8.1-8.2, we provide the identities and estimates that we will use when proving the proposition.
Proposition 8.1 (Integral inequalities for the near-top-order derivatives of g and k). Under the bootstrap assumptions (3.7), there exists a universal constant C * > 0 independent of N and A such that if N is sufficiently large in a manner that depends on A and if ε is sufficiently small, then the following integral inequalities hold for t ∈ (T (Boot) , 1] (where, as we described in Subsect. 1.8.6, constants "C" are allowed to depend on N and other quantities): 8.1. The equations. In this subsection, we derive the evolution equations verified by the neartop-order derivatives of g.
Lemma 8.2 (The equations). Let ⃗ I be a spatial multi-index and let P ≥ 0 be a constant. Then the following commuted metric evolution equations hold: Proof. Equation (8.6a) follows in a straightforward fashion from commuting equation (2.10a) first with ∂ ⃗ I and then with t A+P . Equation (8.6b) follows from applying the same procedure to equation (2.10b).

8.2.
Control of the error terms in the near-top-order energy estimates. We now bound various L 2 norms of the error terms from Lemma 8.2 in terms of the solution norms. Lemma 8.3 (L 2 control of the error terms in the near-below-top-order energy estimates for the g and k). Assume that the bootstrap assumptions (3.7) hold. There exists a constant C > 0 such that if N is sufficiently large in a manner that depends on A and if ε is sufficiently small, then the following estimates hold for t ∈ (T (Boot) , 1].
Error term estimates for the near-top-order derivatives of g. The following estimates hold for the error terms from (8.7a)-(8.7b): Furthermore, the following estimates hold for the error terms from (8.7a)-(8.7b): In addition, the following estimates hold for the error terms from (7.5a)-(7.5b): Let ⃗ I be a spatial multi-index with ⃗ I = N − 1 and let T be the type 0 3 Σ t -tangent tensorfield with the following components relative to the transported spatial coordinates: T eij =∶ −2ng ia ∂ e ∂ ⃗ I k a j . Then the following estimate holds: Error term estimates for the just-below-top-order derivatives of k. The following estimates hold for the error terms from (7.6a) and (7.6b): Let ⃗ I be any spatial multi-index with ⃗ I = N − 1 and let T be the type 1 1 Σ t -tangent tensorfield with the following components relative to the transported spatial coordinates: Then the following estimates hold: Proof. Throughout this proof, we will assume that Aδ is sufficiently small (and in particular that Aδ < σ); in view of the discussion in Subsect. 4.4, we see that at fixed A, this can be achieved by choosing N to be sufficiently large.
The estimate (8.8b) can be proved by applying nearly identical arguments to RHS (8.7b), and we omit the details.
Proof of (8.9a)-(8.9b): We first prove (8.9a). We stress that for this estimate, on RHS (8.7a), we have P = q + σ and ⃗ I = N. To bound RHS (8.7a), we first consider the case in which ⃗ I 3 = N. Using that g g ≲ 1 and g-Cauchy-Schwarz, we deduce that the products under consideration are bounded in the norm ⋅ L 2 g (Σt) by ≲ t A+q+σ n L ∞ (Σt) k ḢN g (Σt) . Next, from (3.1), Defs. 3.5 and 3.6, and the bootstrap assumptions, we deduce that the RHS of the previous expression is ≲ t q+σ−1 H (g,k) (t), which is ≲ RHS (8.9a) as desired. We now consider the case in which ⃗ I 1 = N on RHS (8.7a). Using that g g ≲ 1 and g-Cauchy-Schwarz, we deduce that the products under consideration are bounded in the norm ⋅ L 2 g (Σt) by ≲ t A+q+σ k L ∞ g (Σt) n ḢN g (Σt) . Using Def. 3.5, the estimate (4.1b), and the elliptic estimate (5.3), we deduce that the RHS of the previous expression is ≲ t q+σ−1 L (g,k) (t) + H (g,k) (t) , which is ≲ RHS (8.9a) as desired. It remains for us to consider the remaining cases, in which ⃗ I 1 , ⃗ I 2 , ⃗ I 3 ≤ N − 1 on RHS (8.7a). We first use (4.6) with l = 0 and m = 2 (since (q+σ; ⃗ I) G is type 0 2 ) and (4.3) to deduce (using that ⃗ I = N) that the products under consideration are bounded as follows: Since RHS (8.16) is equal to t −3q times LHS (8.15), the arguments surrounding (8.15) imply that RHS (8.16) ≲ t −6q−2σ−Aδ L (g,k) (t) + H (g,k) (t) , which, in view of (3.1), is ≲ RHS (8.9a) as desired.
The estimate (8.9b) can be proved by applying nearly identical arguments to the products on RHS (8.7b) (with P = q + σ and ⃗ I = N) and we omit the details.
Proof of (8.11b): We stress that for this estimate, on RHS (7.5b), we have P = 2q+σ and ⃗ I = N −1. To bound the first product on RHS (7.5b), we first use g-Cauchy-Schwarz to deduce that it is bounded in the norm . From Defs. 3.5 and 3.6 and the bootstrap assumptions, we see that the RHS of the previous expression is ≲ t 1−10q−σ H (g,k) (t), which, in view of (3.1), is ≲ RHS (8.11b) as desired.
Using essentially the same reasoning, we find that the second and third sums on RHS (7.5b) are bounded in the norm ⋅ L 2 g (Σt) by ≲ (8.18) and thus are ≲ RHS (8.11b) as well. We have therefore proved (8.11b).
The estimate (8.13c) can be proved using a nearly identical argument, the key point being that like k L ∞ F rame (Σt) , the term k L ∞ g (Σt) is bounded by ≲ t −1 .
To bound the last sum on RHS (7.6b), we first use (4.5) to bound (using that ⃗ I = N − 1) the terms under consideration as follows: . From Defs. 3.5 and 3.6, the estimates (4.1a), (4.7a), and (4.7c), and the bootstrap assumptions, we see that which, in view of (3.1), is ≲ RHS (8.13b) as desired.
Proof of (8.13d): We stress that for this estimate, on RHS (7.6b), we have P = 3q+σ and ⃗ I = N −1. We claim that we only have to bound (in the norm ⋅ L 2 g (Σt) ) the second sum on RHS (7.6b) in the cases in which ⃗ I 5 = N −1 or ⃗ I 6 = N −1. For by inspecting the proof of (8.13b) given above, and using (3.1), we see that all remaining products on RHS (7.6b) are bounded in the norm ⋅ L 2 F rame by ≲ t 2q+σ−1 L (g,k) (t) + H (g,k) (t) . Hence, using (4.2) with l = m = 1 (since (Junk;3q+σ; ⃗ I K is type 1 1 ), we find that these same products are bounded in the norm ⋅ L 2 g by ≲ t σ−1 L (g,k) (t) + H (g,k) (t) , which is ≲ RHS (8.13d) as desired.
Proof of (8.14b): The arguments given in the proof of (8.14a) yield that t A+3q+σ , which is a better bound than we need. 8.3. Proof of Prop. 8.1. In this subsection, we prove Prop. 8.1. Throughout this proof, we will assume that Aδ is sufficiently small (and in particular that Aδ < σ); in view of the discussion in Subsect. 4.4, we see that at fixed A, this can be achieved by choosing N to be sufficiently large.
To prove (8.1a), we first use the fundamental theorem of calculus and the evolution equation (8.6a) with P = 2q + σ to deduce that for any multi-index ⃗ I with ⃗ I = N, we have (where there is no summation over i, j and we stress that 0 < t ≤ 1): From Def. 3.6, Cauchy-Schwarz, Young's inequality, and the estimate (8.8a), we deduce that the last integral on RHS (8.25) can be bounded as follows:  (1), summing the resulting estimates over 1 ≤ i, j ≤ D and over ⃗ I with ⃗ I = N, and taking ε to be sufficiently small, we arrive at (8.1a).
The estimate (8.1b) can be proved using a similar argument based on the evolution equation (8.6b) with P = 2q + σ and ⃗ I = N and the estimate (8.8b), we omit the details. The estimate (8.4a) can be proved using a similar argument based on the evolution equation (8.6a) with P = 5q + 3σ − 1 and ⃗ I = N − 1 and the estimate (8.10a), and we omit the details. The estimate (8.4b) can be proved using a similar argument based on the evolution equation (8.6b) with P = 5q + 3σ − 1 and ⃗ I = N − 1 and the estimate (8.10b), and we omit the details. The estimate (8.5a) can be proved using a similar argument based on the evolution equation (7.4b) with P = 3q + σ and ⃗ I = N − 1 and the estimates (8.13a), (8.13b), and (8.14a), and we omit the details.
To prove (8.2a), we let ⃗ I be any any multi-index with ⃗ I = N. Using the evolution equation (8.6a) with P = q + σ, the definition of the norm ⋅ g , and equation (2.10b) (to substitute for the factors of ∂ t g −1 that appear when ∂ t falls on the factors of g −1 that are inherent in the definition of ⋅ g ), we deduce that Next, using equation (8.6a) to substitute for the factor ∂ t (t A+2q+σ ∂ ⃗ I g cd ) on RHS (8.27), we obtain Integrating (8.28) over the spacetime slab (t, 1] × T D , using g-Cauchy-Schwarz, and appealing to Def. 3.6, we obtain the following estimate (where we stress that t < 1): dx ds.
The estimate (8.2b) can be proved using a similar argument based on equation (8.6b) with P = q + σ and ⃗ I = N, equation (2.10a) (which one uses to substitute for the factors of ∂ t g that appear when ∂ t falls on the factors of g that are inherent in the definition of ⋅ g ), and the estimate (8.9b); we omit the details.
The estimate (8.3) can be proved using a similar argument based on the evolution equation (7.4a) with P = 2q + σ, and ⃗ I = N − 1, equation (2.10b) (to substitute for the factors of ∂ t g −1 that appear when ∂ t falls on the factors of g −1 that are inherent in the definition of ⋅ g ), and the estimates (8.11a)-(8.11b) and (8.12). We omit the details, noting only that the factors of ∂ t g −1 and the factor −2tk a j in the first braces on RHS (7.4a) lead to the terms 2nk ad g be g cf (t A+q+σ ∂ a ∂ ⃗ I g bc )(t A+q+σ ∂ d ∂ ⃗ I g ef )+ 4(n − 1)g ad k be g cf (t A+q+σ ∂ a ∂ ⃗ I g bc )(t A+q+σ ∂ d ∂ ⃗ I g ef ), which we pointwise bound in magnitude as follows by using g-Cauchy-Schwarz, the fact that g −1 g ≤ C * , (3.1), (4.1b), Def. 3.5, and the bootstrap assumptions: We further remark that these factors of C * lead to the C * -dependent products on RHS (8.3).
The estimate (8.5b) can be proved using a similar argument based on equation (7.4b) with P = 3q + σ and ⃗ I = N − 1 (where we use equations (2.10a)-(2.10b) to substitute for the factors of ∂ t g and ∂ t g −1 that arise when ∂ t falls on the factors of g and g −1 inherent in the definition of ⋅ g ), and the estimates (8.13c), (8.13d), and (8.14b); we omit the details.

The Main A Priori Estimates
In this section, we use the estimates derived in Sects. 5-8 to prove the main technical result of the article: Prop. 9.2, which provides a priori estimates for the solution norms from Defs. 3.5 and 3.6. The proposition in particular yields a strict improvement of the bootstrap assumptions. 9.1. Integral inequality for the high norm. We start with the following lemma, in which we derive an integral inequality for the high norm H (g,k) (t). The lemma is an analog of Prop. 6.1, in which we derived a similar but simpler inequality for the low norm L (g,k) (t).
Lemma 9.1 (Integral inequality for the high norm). Assume that the bootstrap assumptions (3.7) hold. There exists a universal constant C * > 0 independent of N and A such that if N is sufficiently large and if ε is sufficiently small, then the following integral inequality holds for t ∈ (T (Boot) , 1] (where, as we described in Subsect. 1.8.6, constants "C" are allowed to depend on N and other quantities): Proof. We sum the estimates (7.1) over all ⃗ I with ⃗ I = N together with the estimates (8.1a)-(8.5b). In view of definition (3.6a), we conclude the estimate (9.1).

9.2.
A priori estimates for the solution norms. In the next proposition, we provide the main result of Sect. 9. Proposition 9.2 (A priori estimates for the solution norms). Recall that L (g,k) (t), L (n) (t), H (g,k) (t), and H (n) (t) are the norms from Defs. 3.5 and 3.6, and assume that the bootstrap assumptions (3.7) hold. Letǫ be the following norm of the difference between the Kasner initial data and the perturbed initial data: If A is sufficiently large and if N is sufficiently large in a manner that depends on A, q, σ, and D, then there exists a constant C N,A,q,σ,D > 1 such that if ε is sufficiently small in a manner that depends on N, A, q, σ, and D, then the following estimates hold for t ∈ (T (Boot) , 1]: L (g,k) (t) + H (g,k) (t) + L (n) (t) + H (n) (t) ≤ C N,A,q,σ,Dǫ .

Estimates Tied to Curvature Blowup and the Length of Past-Directed Causal Geodesics
In this section, we derive the main estimates needed to show curvature blowup and geodesic incompleteness for the solutions under study.
10.1. Curvature estimates. In the next lemma, we derive a pointwise estimate that shows in particular that the Kretschmann scalar blows up as t ↓ 0.
Proof. Throughout this proof, we will assume that Aδ is sufficiently small (and in particular that Aδ < σ); in view of the discussion in Subsect. 4.4, we see that at fixed A, this can be achieved by choosing N to be sufficiently large.
First, we observe the following identity, which holds relative to CMC-transported spatial coordinates in view of (2.1) and the curvature properties Riem αβγδ = −Riem βαγδ = −Riem αβδγ = Riem γδαβ : t 4 Riem αβγδ Riem αβγδ = t 4 Riem cd ab Riem ab cd + 4t 4 Riem c0 a0 Riem a0 where Riem 0⋅ is defined to be the type 0 3 Σ t -tangent tensorfield with components Riem 0bcd relative to the transported spatial coordinates. Next, from standard calculations based in part on the Gauss and Codazzi equations (see, for example, [52,Equation (4.14)] and [52,Equation (4.18)]), we find that relative to the CMC-transported spatial coordinates, the components Riem µν αβ of the (type 2 2 ) Riemann curvature tensor of g can be decomposed into principal terms and error terms as follows:

3a)
Riem c0 a0 = t −1 k c a + k c e k e a + △ △ △ c0 a0 , (10.3b) where the error terms are defined by △ △ △ cd ab ∶= Riem cd ab , (10.4a) △ △ △ c0 a0 ∶= −t −1 n −1 ∂ t (tk c a ) + t −1 (n −1 − 1)k c a (10.4b) − n −1 g ec ∂ a ∂ e n + n −1 g ec Γ f a e ∂ f n, △ △ △ cd 0b ∶= g ce ∂ e (k d b ) − g de ∂ e (k c b ) (10.4c) In (10.4a), Riem cd ab denotes a component of the (type 2 2 ) Riemann curvature tensor of of g. We now claim that the following estimates hold, wherek is the Kasner second fundamental form:  (6.15), which are needed to help control the terms on RHS (10.4b); we omit the straightforward details. To prove (10.7), we first use the fact that g −1 g ≲ 1 and the g-Cauchy-Schwarz inequality to deduce that the norm ⋅ g of RHS (10.4c) is bounded by ≲ ∂k g + ∂g g k g . Next, using (4.2), we bound the RHS of the previous expression as follows: ∂k g + ∂g g k g ≲ t −3q k W 1,∞ F rame (Σt) + t −3q g

10.2.
Estimates for the length of past-directed causal geodesic segments. In this subsection, we show that for the solutions under consideration, the length of any past-directed causal geodesic segment is uniformly bounded by a constant.
Lemma 10.2 (Estimates for the length of past-directed causal geodesic segments). Under the hypotheses and conclusions of Prop. 9.2, perhaps enlarging N and shrinkingǫ if necessary, the following holds: any past-directed causal geodesic ζ ζ ζ that emanates from Σ 1 and is contained in the region (T (Boot) , 1] × T D has an affine length that is bounded from above by 1 − q , (10.13) where A (t) is the affine parameter along ζ ζ ζ viewed as a function of t along ζ ζ ζ (normalized by A (1) = 0).
Proof. Throughout this proof, we will assume that Aδ is sufficiently small (and in particular that Aδ < σ); in view of the discussion in Subsect. 4.4, we see that at fixed A, this can be achieved by choosing N to be sufficiently large.
Let ζ ζ ζ(A ) be a past-directed affinely parametrized geodesic verifying the hypotheses of the lemma. Note that the component ζ ζ ζ 0 can be identified with the CMC time coordinate. In the rest of the proof, we view the affine parameter A as a function of t = ζ ζ ζ 0 along ζ ζ ζ. We normalize A (t) by setting A (1) = 0. We also defineζ ζ ζ µ ∶= d dA ζ ζ ζ µ and A ′ = d dt A . By the chain rule, we have (10.14) For use below, we note that since ζ ζ ζ is a causal curve, we have (by the definition of a causal curve) that g(ζ ζ ζ,ζ ζ ζ) ≤ 0. Considering also the expression (2.1) for g, we deduce that relative to the CMC-transported coordinates, causal curves satisfy Moreover, from Def. 3.5, (3.1), and the estimate (9.3), we see that the last product on RHS (10.18) obeys the bound Furthermore, also using (4.1a), (4.2), (4.7c), and (5.2), we see that the last term on RHS (10.17) is bounded as follows: Integrating (10.22) from time t to time 1 and using the assumption A (1) = 0, we find that t ∈ (T (Boot) , 1], (10.23) from which the desired estimate (10.13) follows.

The Main Stable Blowup Theorem
We now state and prove our main stable blowup theorem. As we noted in Remark 1.2, it is possible to derive substantial additional information about the solution, going beyond that provided by the theorem.
• N > 0 is sufficiently large, where the required largeness depends on A, q, σ, and D. Here we recall that q > 0 and σ > 0 are the constants fixed in Subsect. 3.1. •ǫ is sufficiently small, where the required smallness depends on N, A, q, σ, and D.
Then the following conclusions hold.
Existence and norm estimates on (0, 1] × T D . The initial data launch a solution (g, k, n) to the Einstein-vacuum equations in CMC-transported spatial coordinates (that is, the equations of Prop. 2.1, where g = −n 2 dt ⊗ dt + g ab dx a ⊗ dx b is the spacetime metric) that exists classically for (t, x) ∈ (0, 1] × T D . Moreover, there exists a constant C > 1 (depending on N, A, q, σ, and D) such that the (N, A, q, σ)-dependent norms from Defs. 3.5 and 3.6 verify the following estimate for t ∈ (0, 1]: L (g,k) (t) + H (g,k) (t) + L (n) (t) + H (n) (t) ≤ Cǫ. In particular, forǫ sufficiently small, Riem αβγδ Riem αβγδ blows up like {C+O(ǫ)}t −4 as t ↓ 0 (where C = 4 ∑ D i=1 (q 2 i − q i ) 2 + ∑ 1≤i<j=D q 2 i q 2 j ). Consequently, the maximal (classical) globally hyperbolic development of the data is ((0, 1] × T D , g), and g cannot be continued as a C 2 Lorentzian metric to the past of the singular hypersurface Σ 0 . That is, the past of Σ 1 in the maximal (classical) globally hyperbolic development of the data is foliated by the family of spacelike hypersurfaces Σ t , along which the CMC condition k a a = −t −1 holds. Bounded length of past-directed causal geodesics. Let ζ ζ ζ be any past-directed causal geodesic that emanates from Σ 1 , and let A = A (t) be an affine parameter along ζ ζ ζ, where A is viewed as a function of t along ζ ζ ζ that is normalized by A (1) = 0. Then ζ ζ ζ crashes into the singular hypersurface Σ 0 in finite affine parameter time Proof. We first fix A and N to be large enough so that all of the estimates proved (under bootstrap assumptions) in the previous sections hold true. By standard local well-posedness, ifǫ is sufficiently small and C ′ is sufficiently large, then there exists a maximal time T (M ax) ∈ [0, 1) such that the solution (g, k, n) exists classically for (t, x) ∈ (T (M ax) , 1] × T D and such that the bootstrap assumptions (3.7) hold with T (Boot) ∶= T (M ax) and ε ∶= C ′ǫ . By enlarging C ′ if necessary, we can assume that C ′ ≥ 2C N,A,q,σ,D , where C N,A,q,σ,D > 1 is the constant from inequality (9.3). Readers can consult [4] for the main ideas behind the proof of local well-posedness, in a similar but distinct gauge for Einstein's equations. Moreover, in view of Defs. 3.5 and 3.6, it is a standard result (again, see [4] for the main ideas) that if ε is sufficiently small, then either T (M ax) = 0 or the bootstrap assumptions are saturated on the time interval (T (M ax) , 1], that is, sup t∈(T (M ax) ,1] L (g,k) (t) + H (g,k) (t) + L (n) (t) + H (n) (t) = C ′ǫ . The latter possibility is ruled out by inequality (9.3). Thus, T (M ax) = 0. In particular, the solution exists classically for (t, x) ∈ (0, 1]×T D , and the estimate (11.3) holds for t ∈ (0, 1]. The remaining aspects of the theorem follow from Lemmas 10.1 and 10.2.