The relativistic Euler equations: ESI notes on their geo-analytic structures and implications for shocks in 1D and multi-dimensions

In this article, we provide notes that complement the lectures on the relativistic Euler equations and shocks that were given by the second author at the program Mathematical Perspectives of Gravitation Beyond the Vacuum Regime, which was hosted by the Erwin Schrödinger International Institute for Mathematics and Physics in Vienna in February 2022. We set the stage by introducing a standard first-order formulation of the relativistic Euler equations and providing a brief overview of local well-posedness in Sobolev spaces. Then, using Riemann invariants, we provide the first detailed construction of a localized subset of the maximal globally hyperbolic developments of an open set of initially smooth, shock-forming isentropic solutions in 1D, with a focus on describing the singular boundary and the Cauchy horizon that emerges from the singularity. Next, we provide an overview of the new second-order formulation of the 3D relativistic Euler equations derived in Disconzi and Speck (2019 Ann. Henri Poincare 20 2173–270), its rich geometric and analytic structures, their implications for the mathematical theory of shock waves, and their connection to the setup we use in our 1D analysis of shocks. We then highlight some key prior results on the study of shock formation and related problems. Furthermore, we provide an overview of how the formulation of the flow derived in Disconzi and Speck (2019 Ann. Henri Poincare 20 2173–270) can be used to study shock formation in multiple spatial dimensions. Finally, we discuss various open problems tied to shocks.

In this article, we provide notes that complement the lectures on the relativistic Euler equations and shocks that were given by the second author at the program Mathematical Perspectives of Gravitation Beyond the Vacuum Regime, which was hosted by the Erwin Schrödinger International Institute for Mathematics and Physics in Vienna in February, 2022. Broadly speaking, our main goals are the following: (1) To briefly introduce the equations. (2) To provide an overview of the new formulation of the equations derived in [41]. (3) To describe the implications of the formulation from [41] for the study of multi-dimensional shocks.
Although our main interest is solutions in 3D without symmetry assumptions, we highlight Sect. 3, in which we rigorously study a family of initially smooth, simple isentropic shock-forming plane-symmetric solutions and give a complete description of a portion of their maximal (classical) globally hyperbolic developments (MGHD for short), up to the boundary. To the best of our knowledge, this is the first article on any quasilinear hyperbolic system in 1D that provides the detailed construction of the MGHD within the vicinity of a shock singularity. Sect. 3 in particular provides an introduction to some of the main ideas behind the study of shock formation in a gentle, semi-explicit 1D setting, where energy estimates are not needed. Roughly, the MGHD is the largest possible classical solution + globally hyperbolic region that is launched by the initial data, and when it exists and is unique (for some solutions to some hyperbolic PDEs, it is not unique -see Sect. 5.5!), it is the holy grail object, at least from the point of view of classical solution.
We stress outright that, although the relativistic Euler equations are posed on an ambient Lorentzian manifold (M, g), from the point of view of the causal structure of the fluid, the correct notion of "globally hyperbolic" is not with respect to the spacetime metric g, but rather with respect to the acoustical metric h, introduced in Def. 2.4. To avoid confusion, we will refer to this as h-global hyperbolicity and the corresponding h-MGHD.
Many of the techniques and geometric insights behind the modern approach to studying shocks have roots in mathematical General Relativity, notably the celebrated proof by Christodoulou-Klainerman [32] of the nonlinear stability of Minkowski spacetime as a solution to Einstein's equations. A key unifying theme between the global existence result [32] and the results on shocks that we discuss here is that nonlinear geometric optics, implemented via an eikonal function U , plays a central role in the study of the flow. We will discuss eikonal functions in detail in Sects. 3 and 6.
We highlight a key advantage of the formulation of the flow derived in [41].
For general solutions (in 3D, i.e., three spatial dimensions), it allows one to implement a sharp version of nonlinear geometric optics. Nonlinear geometric optics is crucial for the study of multi-dimensional shocks, and it has other applications, such as low-regularity wellposedness (see Sect. 2.3.5).
In line with points (1)-(3) above, this article is meant to give an introduction to the mathematical techniques and methods needed for the PDE analysis through the lens of nonlinear geometric optics. A comprehensive topical review of the rich history of shock singularities in the framework of hyperbolic conservation laws would encompass a work of strictly larger order of magnitude. We refer the reader to [18,39,49] and the references therein for a good start.
1.1. Motivation, context, and the structure of the article. For a non-exhaustive yet comprehensive overview of these physical applications, we refer the reader to the texts [78,97,100], and the references therein. There has also been exciting rigorous mathematical progress on the global structure of solutions when the spacetime is expanding [43,44,77,84,90]. We also highlight [29] and [30,Chapter 1] for a modern introduction to the relativistic Euler equations and their connections to the laws of mechanics and thermodynamics.
Despite the above remarks, the rigorous mathematical theory of multi-dimensional solutions is far from complete. Even if one considers the relativistic Euler equations on a given spacetime background (i.e., without coupling to Einstein's equations), there are many mathematically rich phenomena that have not been fully understood. Chief among these is the fundamental issue, going back to Riemann's foundational work [83] on non-relativistic compressible fluids in 1D, that initially smooth solutions can develop shock singularities in finite time. Shocks are singularities such that various fluid variables' gradients blow up in finite time (in a precise, controlled fashion, as it turns out), though the fluid variables themselves (such as the velocity and density) remain bounded. This phenomenon is also known as wave breaking. The relatively mild nature of the singularity gives rise to the hope that one might -at least for short times -be able to uniquely continue the solution weakly past the "first singularity, 1 " subject to suitable selection criteria, typically in the form of Rankine-Hugoniot-type jump conditions across a shock hypersurface (which is not known in advance) and an entropy-type condition. This weak continuation problem is known as the shock development problem, and in full generality, it remains open, for both the 3D relativistic Euler equations and their non-relativistic analog, i.e., the 3D compressible Euler equations; see Sect. 5.6 for further discussion. The long-term goal is certainly the following: Develop a rigorous global-in-time-and-space theory of existence and uniqueness for open sets of multi-dimensional solutions that are allowed to transition from classical to weak due to shock formation, and describe the interactions of different shock hypersurfaces, Cauchy horizons (see Sect. 3.8), and other crucial qualitative features of the flow. This goal is far out of reach as of present. If one could accomplish it even in a single perturbative regime, that would represent truly remarkable progress.
While there have been many works on shock formation for the multi-dimensional non-relativistic compressible Euler equations, such as [1-3, 21-24, 34, 65, 67, 91, 95], there have been relatively few works on multi-dimensional shock formation for the relativistic Euler equations. The most notable work in the relativistic case is Christodoulou's breakthrough monograph [30], in which he proved shock formation and studied some aspects of the h-MGHD for open sets of irrotational and isentropic solutions in 3D that are compactly supported perturbations of non-vacuum constant fluid states. Unlike the non-relativistic case, there are currently no rigorous results that prove multi-dimensional shock formation for the relativistic Euler equations in the presence of vorticity and entropy. However, in Sect. 6, we provide a blueprint for how one can use the equations of [41] and the analytic framework of [1,3,30,65,67,91,95] to prove stable shock formation in the 3D relativistic case for open sets of solutions with vorticity and entropy.
The remainder of the article is organized as follows: • In Sect. 2, we set up the study of the relativistic Euler equations and provide an overview of basic ingredients that play a role in the proof of local well-posedness in Sobolev spaces. In particular, we introduce the acoustical metric h, which is the solution-dependent Lorentzian metric that drives the propagation of sound waves; the acoustical metric is fundamental for implementing nonlinear geometric optics and studying multi-dimensional shocks. While we state the equations in the case of an arbitrary spacetime background (M, g), the vast majority of the article concerns the case in which M = R 1+3 and g = m is the Minkowski metric on M. We again point out that nearly all notions of hyperbolicity will refer to the acoustical metric h and not g or m; see also Sect. 4.2.5 and Remark 6.1. • In Sect. 3, to help readers gain intuition for the subsequent discussion, we study the flow in one spatial dimension, which is equivalent to studying plane-symmetric solutions in 3D. In particular, in Theorem 3.1, we provide a detailed analysis of the h-MGHD for a large set of shock-forming relativistic Euler solutions in 1D. Some aspects of the theorem, notably the results concerning the formation of Cauchy horizons, have not been proved in detail elsewhere in the literature. Although Theorem 3.1 specifically concerns the relativistic Euler equations, the techniques could be applied to a wide variety of strictly hyperbolic genuinely nonlinear PDEs in one spatial dimension. • In Sect. 4, we provide an overview of the new formulation of the 3D relativistic Euler equations derived in [41].
• In Sect. 5, we describe some previous works on shock formation and related problems.
• In Sect. 6, based in part on our prior experience [1,3,65,67,91,95] in studying shock formation in multi-dimensional non-relativistic compressible fluids, we provide an overview of how the new formulation of the flow from [41] can be used to study shock formation for open sets of solutions to the 3D relativistic Euler equations. As of present, rigorous fluid shock formation results in multi-dimensions with vorticity and entropy have been proved only for the non-relativistic compressible Euler equations. Nonetheless, from the perspective of shock formation, the relativistic Euler equations and the non-relativistic compressible Euler equations enjoy many structural commonalities, and we expect that the overview we provide in Sect. 6 could be turned (with substantial effort) into a complete proof. • In Sect. 7, we describe various open problems.

Standard formulations of the relativistic Euler equations and local well-posedness
In this section, we introduce some standard first-order formulations of the relativistic Euler equations. We also introduce the acoustical metric h and set up a corresponding geometric version of the energy method that applies to a first-order formulation of the flow. Finally, we state Prop. 2.3, which is a standard result on local well-posedness, and briefly discuss how its proof is connected to the energy method.
2.1. The basic setup. We start by discussing standard formulations of the relativistic Euler equations on an arbitrary four-dimensional Lorentzian manifold (M, g), where g is the spacetime metric of signature (−, +, +, +). Our discussion in this section is motivated by Christodoulou's presentation in [30,Chapter 1]. We are somewhat terse in our presentation here, and thus we refer readers to [30,Chapter 1] for additional details.
Convention 2.1 (Moving indices with g and m). In this section, we will lower and raise indices with the spacetime metric g and its inverse g −1 respectively, i.e., u α := g αβ u β and ξ α := (g −1 ) αβ ξ β . Similar comments apply to the bulk of the article, once we fix M = R 1+3 and g = m as the Minkowski metric on R 1+3 ; see Sect. 2.2.1. This will become especially important later on, when we introduce the acoustical metric h, i.e., we do not raise or lower indices with h or h −1 .
denotes the enthalpy per particle, and u α denotes the fluid's four-velocity, which is a future-directed vectorfield on M normalized by: In the rest of the paper, denotes an arbitrary fixed, positive value of H; we find it convenient to use H to normalize various constructions. We find it convenient to work with the natural log of the enthalpy.
Definition 2.1 (Logarithmic enthalpy). We define the (dimensionless) logarithmic enthalpy h as follows: Definition 2.2 (The quantity q). We define the quantity q as follows: 2.1.2. A first statement of the relativistic Euler equations. Relative to arbitrary coordinates, the relativistic Euler equations on (M, g) can be expressed as the following system of quasilinear hyperbolic conservation laws: ∇ κ T ακ = 0, (α = 0, 1, 2, 3), (6a) where ∇ is the Levi-Civita connection of the spacetime metric g and: T αβ := (ρ + p)u α u β + p(g −1 ) αβ , (α, β = 0, 1, 2, 3), is the energy-momentum tensor of the fluid. It turns out that under additional assumptions discussed in Sect. 2.1.3, equation (2) can be viewed as a constraint that is preserved by the flow of the PDEs (6a); see equation (15) and the discussion surrounding it. We also emphasize that later on, we will discuss several other formulations of the flow.
2.1.3. Equation of state and thermodynamic relations. The system (2) + (6a)-(6b) + (7) is not closed because there are too many fluid variables and not enough equations. The standard path to closing the system begins with an equation of state, that is, an assumed functional relationship of the form p = p(ρ, s).
The basic hyperbolicity of the equations will require that speed of sound c, defined by: should be real and non-negative, where ∂p ∂ρ | s denotes the partial derivative of p (i.e., of the equation of state) with respect to ρ at fixed s. The vanishing of c causes a severe degeneracy in the system, and we will therefore restrict our attention to solution regimes in which: where we have normalized our setup so that the speed of light is unity (i.e., the second inequality in (9) implies that the speed of sound is less than or equal to the speed of light). We also restrict our attention to solutions such that: ρ > 0, p > 0, n > 0, θ > 0, H > 0.
The laws of thermodynamics demand that the fluid variables satisfy the following functional relations: H = ∂ρ ∂n | s , θ = 1 n ∂ρ ∂s | n , dH = dp n + θds, where ∂ ∂n | s denotes partial differentiation with respect to n at fixed s and ∂ ∂s | n denotes partial differentiation with respect to s at fixed n.

2.2.
Minkowski metric assumption and Minkowski-rectangular coordinates. In the rest of the paper, unless we explicitly state otherwise, it should be understood that the spacetime manifold M is equal to R 1+3 and that g is equal to the Minkowski metric, which for clarity we denote by m. We fix a standard global Minkowski-rectangular coordinate system {x α } α=0,1,2,3 on M, relative to which m = diag(−1, 1, 1, 1). We use the notation {∂ α } α=0,1,2,3 to denote the partial derivative vectorfields in this coordinate system. We also use the alternate notation t := x 0 and ∂ t := ∂ 0 .
2.2.1. Index conventions. From now until the end of the paper, we use the following conventions for indices.
• (Uppercase Latin index conventions) Uppercase Latin spatial indices A, B, etc. correspond to the spatial coordinates (x 2 , x 3 ) and vary over 2, 3. This will be important in Sect. 6, where (x 2 , x 3 ) will correspond to perturbations away from plane-symmetry. • (Einstein summation) We use Einstein's summation convention in that repeated indices are summed, e.g., L A X A := L 2 X 2 + L 3 X 3 .

2.3.
A well-posed first-order formulation. In this section, we provide a first-order formulation of the flow in which the unknowns are taken to be h, u 0 , u 1 , u 2 , u 3 , s . Thanks to the assumptions state in Sect. 2.1.2, the remaining fluid variables can be expressed 2 as functions of the unknowns. With the help of (11), one can compute that for C 1 solutions, (6a)-(6b) is equivalent to the following first-order formulation of the flow: where: denotes projection onto the m-orthogonal complement of u. In particular, Equations (12b)-(12d) + (13) are a first-order quasilinear hyperbolic system and are locally well-posed in suitable Sobolev spaces; see Prop. 2.3. using (12c) We also note by contracting (12b) against u β and using (12c), we deduce -without using (12d) -the following identity: Using (15), one can easily show via Grönwall's inequality that if m(u, u)+1 vanishes at t = 0, then m(u, u)+1 vanishes on any slab [0, T ] × R 3 of classical existence on which the solution is C 1 . Hence, (12d) can be viewed as a constraint on the initial data that is preserved by the flow.
Remark 2.1 (Additional terms in the context of curved spacetime backgrounds). When posed on a general smooth spacetime background (M, g), there are additional terms on RHS (12a)-(12b) of the schematic form Γ · f(V), where Γ schematically denotes Christoffel symbols of g, V denotes the array of fluid variable unknowns (see definition 23), and f schematically denotes a smooth function. In the context of fluid shock formation, the dominant terms in the flow are Riccati-type nonlinearities in the fluid obtained after differentiating the equations one time, i.e., they are quadratic in the first derivatives of the fluid; see Remark 3.5. Note that upon differentiating the term Γ · f(V), one obtains a term that is linear in the derivatives of the fluid variables. This means, in particular, that the for large-gradient solutions, the differentiated terms ∂ α (Γ · f(V)) will be small relative to the shock-driving Riccati-type fluid terms, i.e., the "new terms" coming from the non-flat ambient geometry are relatively small. This observation is a key reason behind our expectation that one can use the methods described later in this article to prove shock formation for the relativistic Euler equations on a general spacetime (M, g), at least for initial data with suitable large gradients.
2.3.1. Equations of variation. As a preliminary step in studying the well-posedness of (12b)-(12d), we will discuss the equations satisfied by the derivatives of solutions. This is important in the sense that proofs of local well-posedness rely on differentiating the equations to obtain energy estimates for the solution's higher derivatives. Specifically, by differentiating the equations, it is straightforward to see that given a solution (h, u α , s) to (12a)-(12d), any derivative (of any order), which we denote by (ḣ,u α ,ṡ), satisfies an inhomogeneous system, known as the equations of variation, of the following form: Although here we do not provide a detailed formula for the structure of the inhomogeneous termsḞ, · · · ,İ on RHSs (16a)-(16d), it is easy to see that they feature terms involving ≤ the number of derivatives represented by the variations (ḣ,u α ,ṡ), i.e., fewer derivatives than the principal terms on LHSs (16a)-(16d).

The acoustical metric.
A basic question about the equations of variation is: what kinds of energy estimates are available for solutions? It turns out that this question is connected to the intrinsic geometry of the flow via the acoustical metric, which plays a fundamental role throughout the article.
Definition 2.4 (The acoustical metric and its inverse). We define: We refer to h as the acoustical metric.
Remark 2.4 (Acoustical metric in a general ambient spacetime (M, g)). When studying the relativistic Euler equations in a general ambient spacetime (M, g) equipped with coordinates (x 0 , x 1 , x 2 , x 3 ), where t := x 0 is a smooth time function (see Sect. 4.2.5), the acoustical metric is (modulo the comments of Remark 2.2): where where the positivity on RHS (22) stems from (9) and the assumption that t is a time function for g, which in particular implies that relative to the coordinates (x 0 , x 1 , x 2 , x 3 ), we have (g −1 ) 00 = (g −1 ) αβ (∂ α t)∂ β t < 0.

2.3.3.
The arrays V andV. For notational convenience, we define the solution array V and the variation arrayV as follows: 2.3.4. Energy currents and energy identity in differential form. We now introduce energy currents, which form the basis for energy identities for the equations of variation.
Definition 2.5 (Energy currents). Given a solution V to (12b)-(12d) + (13), variationsV = (ḣ,u α ,ṡ) that satisfy (16a)-(16d), and scalar functions f 1 and f 2 , we define the corresponding energy current vectorfielḋ J (V,f1,f2) as follows: Note that the components ofJ can be viewed as quadratic forms inV with coefficients that depend on the solution V. Related energy currents were used by Christodoulou in [30]. A general framework for constructing energy currents in the context of regularly hyperbolic PDEs was developed by Christodoulou in [28].
The following lemma provides a differential version of the energy identities afforded by the energy current. The proof is a tedious but straightforward computation, omitted here (c.f. [ where Q (V,f1,f2,∂V,∂f1,∂f2) [V,V] is a quadratic form inV with coefficients that depend on V, f 1 , f 2 , and their first-order partial derivatives; it is straightforward to compute Q in detail, though we do not do so here.
By integrating (26) over a spacetime region and applying the divergence theorem, one obtains an energy identity, where the energies are defined along the boundary of the region. A fundamental question is: when are the energies obtained in this fashion positive definite 3 ? The answer to this question, provided by Lemma 2.2, is tied to the acoustic geometry of the system, captured by the acoustical metric. In short, the energies are positive definite whenever the boundary of the region has h-timelike normals, i.e., whenever the boundary is spacelike with respect to h.
Before proving the lemma, we first introduce a subset of state-space in which the relativistic Euler equations are hyperbolic. Definition 2.6 (Regime of hyperbolicity). We define the regime of hyperbolicity H to be the following subset of solution variable space: . Let ξ be a past-directed h-timelike one-form, i.e., assume that (h −1 ) κλ ξ κ ξ λ < 0 and ξ 0 := ξ(∂ t ) > 0. Let K be a compact subset of the set H defined in (27), and assume that the solution array V defined in (23) satisfies V ∈ K. Then if f 1 and f 2 are sufficiently large (in a manner that depends on ξ and K), there exists a C = C(ξ, K, f 1 , f 2 ) > 0 such that the quadratic form ξ αJ α (V,f1,f2) [V,V] is positive definite in the following sense: We define the energy current vectorfield J to be the vectorfield with the following components relative the standard Minkowski-rectangular coordinates: In (29a), ∇ x ϕ := (∂ 1 ϕ, ϕ 2 , ϕ 3 ) denotes the spatial gradient of ϕ and | · | is the standard Euclidean norm on R 3 . Then straightforward calculations yield that ∂ α J α = 0 and that (∂ t ) α J α = 1 2 (∂ t ϕ) 2 + 1 2 |∇ x ϕ| 2 . These two identities are analogs of (26) and (28) respectively. The upshot is that by integrating ∂ α J α = 0 over [0, t]×R 3 , one obtains the celebrated conservation of energy formula for solutions to the linear wave equation: The approach is called the "multiplier method" because it is a geometric version of the more standard proof of (30), in which one multiplies the wave equation by ∂ t ϕ and integrates by parts over [0, t] × R 3 .
Proof. Throughout, C = C(ξ, K, f 1 , f 2 ) denotes a positive constant (often smaller than 1) that can vary from line to line. We start by defining the scalar function ξ ∥ by: Then, using (12d), we decompose: where the one-form ξ ⊥ α satisfies ξ ⊥ κ u κ = 0. In particular, ξ ⊥ α is m-orthogonal to the m-timelike vectorfield u α and thus ξ ⊥ α is m-spacelike whenever it is non-zero, i.e., (m) −1 ξ ⊥ κ ξ ⊥ λ ≥ 0. Since (h −1 ) κλ ξ κ ξ λ < 0 and h κλ u κ u λ < 0 (see (20)), and since ξ α is past-directed while u α is future-directed, it follows that there is a C with 0 < C < 1 such that: We can similarly decomposeu α into a variation that is parallel to u α and one that vanishes upon contraction with u α . We claim that prove (28), it suffices to show that when f 1 = f 2 = 0,ṡ = 0, and u κu κ = 0, we have: The estimate (34) would yield control overḣ 2 and variations such thatṡ = u κu κ = 0 (note thatu α is m-spacelike for such variations and thus the termu κu κ on RHS (34) is positive-definite in such variations). Then to handle general variations, since the first and last terms on RHS (25) yield the two coercive (in view of (33)) terms f 1 ξ ∥ṡ2 and f 2 ξ ∥ (u κu κ ) 2 , we can choose f 1 and f 2 to be sufficiently large (for example, we can choose them to be large, positive constants) to obtain control overṡ 2 and (u κu κ ) 2 , and large enough so that these two coercive terms and the coercivity (34) can collectively be used to absorb (via Young's inequality) all the remaining terms on RHS (25). In total, this would yield (28).
We make the following remarks.
• In Sects. 3 and 6, we will study the formation of shock singularities, in which the solution remains in the interior of the regime of hyperbolicity and the breakdown is precisely due to the scenario (39). • The techniques of [50], which rely only on energy estimates, Sobolev embedding, and the fractional Sobolev-Moser calculus, can be used to save half a derivative compared to Prop. 2.3, that is, to prove local well-posedness for the 3D relativistic Euler equations wheneverh,ů 1 ,ů 2 ,ů 3 ,s ∈ H 5/2 + (R 3 ). • In [102], a low-regularity local well-posedness result for the 3D relativistic Euler equations was proved. The low-regularity assumption was on the "wave-part" of the initial data, which was assumed to belong to H 2 + (R 3 ), while the "transport-part" (including the vorticity and entropy) was assumed to satisfy some additional smoothness assumptions, which were shown to be propagated by the flow; the terminology "wave-part" and "transport" is tied to the formulation of the flow provided by Theorem 4.1 and is explained in [102]. Lindblad's work [64] implies that the assumptions on the wave-part of the initial data are optimal in the scale of Sobolev spaces due to the phenomenon of instantaneous shock formation for initial data in H 2 (R 3 ). The analysis in [102] relies on the new formulation of the flow presented in Theorem 4.1, which allows one to employ analytical techniques based on nonlinear geometric optics and Strichartz estimates. The techniques used in [102] have roots in the works [40,99] (see also [103][104][105]) on low-regularity well-posedness for the non-relativistic compressible Euler equations, which in turn have roots in the works [55-60, 86, 98] on low-regularity well-posedness for quasilinear wave equations.
Discussion of proof of Prop. 2.3. This is a standard result that can be proved by linearizing the equations, deriving a priori estimates for solutions to the linearized equations, and then invoking a standard iteration scheme. We refer to [87][88][89] for the main ideas, where an analog of Prop. 2.3 was proved for the coupling of the relativistic Euler equations to Nordström's theory of gravity. The main idea behind the a priori estimates (in H N ) is to construct energy currents for the linearized systems such that analogs of (26) and (28) hold, and to control RHS (26) in L 2 using the Sobolev-Moser calculus (which is standard in 3D when N ≥ 3). □

A study of the maximal development in simple isentropic plane-symmetry
Our goal in this section is to prove Theorem 3.1, which yields a complete description of a localized portion of the boundary of the maximal (classical) h-globally hyperbolic development (h-MGHD for short) for open sets of shock-forming solutions to the relativistic Euler equations (12a)-(12c) in simple isentropic plane-symmetry; see Fig. 2. Roughly, the h-MGHD is the largest possible classical solution + h-globally hyperbolic region that is launched by the initial data. By a h-globally hyperbolic region R of classical existence launched by the initial data, we mean that every inextendible h-causal curve in R intersects Σ 0 := {t = 0}, where h is the acoustical metric from Def. 2.4. In particular, the h-MGHD cannot contain any points lying in the h-causal future of any singularity, where the h-causal future of a point is effectively determined by integral curves of the two transversal null vectorfields L and L from Def. 3.1.
The portion of the boundary that we study includes a singular boundary, where the fluid's gradient blows up in a shock singularity, and a Cauchy horizon, along which the solution remains smooth. For the solutions under study, the past boundary of the Cauchy horizon intersects the past boundary of the singular boundary in a distinguished point that we refer to as "the crease;" see Fig. 3. The crease plays the role of the true initial singularity in a shock-forming relativistic Euler solution.
Sect. 3 is intended to provide a gentle introduction to the main ideas behind the study of the h-MGHD in a simplified setting in which the geometry is easy to study and one can avoid the burden of energy estimates. It will also help set the stage for Sect. 6, in which we outline some of the main ideas needed to study shock formation in 3D.

3.1.
A description of the symmetry class. As we mentioned, we will study simple isentropic planesymmetric solutions, which we now describe. First, by "isentropic solutions," we mean that s ≡ 0; for any constant s 0 , we could treat the case s ≡ s 0 using the same arguments. From (12c), we infer that for C 1 isentropic solutions, the condition s = s 0 is preserved by the flow if it is satisfied at t = 0. By "planesymmetry," we mean that u 2 = u 3 ≡ 0 and that h = h(t, x 1 ), u 1 = u 1 (t, x 1 ), and u 0 = u 0 (t, x 1 ). Such solutions can be viewed as solutions in 1D in which s ≡ 0, and throughout the rest of the paper, we make that identification. By "simple," we mean that only one Riemann invariant (see Sect. 3.4) is non-vanishing.
3.2. The basic setup and conventions in plane-symmetry. In Sect. 3, since we are studying isentropic plane-symmetric solutions, we will slightly abuse notation and use h to denote the restriction of the acoustical metric (17a)-(17b), to the (t, x 1 )-plane. That is, in Sect. 3, relative to the (t, x 1 )-coordinates, we have: where n is defined by (17c), and we recall that (19) holds. For future use, we note that for C 1 isentropic plane-symmetric solutions, the relativistic Euler equations (12a)-(12c) in Minkowski spacetime are equivalent to: For future use, we also note the following algebraic identity, which holds in plane-symmetry due to (17c) and (41c): To ensure ensure that shocks can form in solutions near constant fluid states with constant enthalpy H > 0, we assume the following non-degeneracy assumption: where in (43), we are viewing c = c(H), c ′ := d dH c, and LHS (43) is defined to be the constant obtained by evaluating 1−c 2 + c ′ Hc at H = H. The condition (43) implies in particular that the isentropic plane-symmetric relativistic Euler equations are genuinely nonlinear near the constant-state solution (H, u 1 ) ≡ (H, 0); see also Sect. 3.5. Christodoulou showed [30, Chapter 1] that irrotational relativistic fluids (of which isentropic plane-symmetric fluids are a special case) admit a Lagrangian formulation for a potential function, and the only equation of state leading to 1−c 2 + c ′ Hc ≡ 0 is such that the potential function solves the timelike minimal surface equation. 4 For that equation, which satisfies a nonlinear version of the null condition that is stronger than the classic null condition identified by Klainerman [54], the first author proved that perturbations of plane-symmetric data (which are similar to the data we study here in Sect. 3) lead to global, nonlinearly stable C 2 solutions in both R 1+1 and R 1+3 [4,5].
3.3. Null frame. In 1D, it is possible to explicitly write down a pair of transversal null vectorfields that dictate the geometry of sound waves. As we explain in Sect. 6.3, in more than one spatial dimension, it is not possible to explicitly write down null vectorfields in terms of the fluid variables. Rather, in that context, the null vectorfields depend on the gradient of an eikonal function U . Definition 3.1 (Null vectorfields in 1D). Relative to the Minkowski-rectangular coordinates (t, x 1 ), we define the vectorfields L and L as follows: Straightforward calculations based on (40a) and (44a)-(44b) yield the following identities: In particular, (45) shows that L and L are h-null. The pair {L, L} is known as a h-null frame. From (45) and straightforward calculations, we deduce the following lemma.

Riemann invariants.
In this section, we introduce the Riemann invariants, which are a pair of fluid variables, one of which is constant along the integral curves of L and the other along the integral curves of L. For isentropic plane-symmetric solutions to the relativistic Euler equations, the Riemann invariants can be used as state-space variables for the fluid, i.e., the other fluid variables are determined by them. Riemann invariants were discovered by Riemann in his proof of shock formation for the non-relativistic compressible Euler equations in 1D. In [96], Taub discovered Riemann invariants for the isentropic relativistic Euler equations in 1D. Studying the flow with respect to the Riemann invariants allows one to study the h-MGHD using techniques based on transport equations and ODE theory. In particular, energy estimates are not needed in 1D. Although Riemann invariants are not available in multi-dimensions, approximate Riemann invariants are available, and they are convenient for studying perturbations of isentropic plane-symmetric solutions; see Sect. 6.2.
To proceed, we let F = F (H) be the solution to the initial value problem: where H > 0 is the constant (3), and in (47), we are viewing c as a function of H. Note that since c(H) ̸ = 0 by assumption, it follows from (47) that the function F has a local inverse F −1 defined near the origin such that F −1 (0) = H. Definition 3.2 (Riemann invariants). We define the Riemann invariants R PS (+) and R PS (−) as follows: In the next proposition, for isentropic plane-symmetric solutions, we express the relativistic Euler equations in terms of (R PS (+) , R PS (−) ), and we provide explicit formulae for various quantities in terms of (R PS (+) , R PS (−) ). Proposition 3.2 (Isentropic plane-symmetric flow in terms of the Riemann invariants). Let (h, u 0 , u 1 ) be a C 1 solution to (41a)-(41c), and let L and L be the vectorfields from Def. 3.1. Then the Riemann invariants of Def. 3.2 satisfy the following quasilinear transport system: Moreover, the components of the fluid velocity in the Minkowski-rectangular coordinate system (t, x 1 ) can be expressed in terms of the Riemann invariants as follows: In addition, with F −1 denoting the inverse function of the map F from Def. 3.2, we have: Finally, the following identities hold: Proof. Equations (49a)-(49b) follow from straightforward but tedious calculations involving Def. 3.2, equations (41a)-(41c), the identity u 0 ∂ α u 0 = u 1 ∂ α u 1 (which follows from differentiating (41c) with ∂ α ), and the identity ∂ α F (H) = 1 c ∂ α h (which follows from (4) and (47)). To prove the identities in (50a), we first subtract (48b) from (48a) and use (41c) to compute that = ln(u 0 + u 1 ). Rewriting (41c) as (u 0 + u 1 )(u 0 − u 1 ) = 1 and carrying out straightforward calculations based on the definitions of cosh and sinh, we conclude (50a). (50b) follows from adding (48a) to (48b), dividing by 2, and then applying F −1 to each side of the resulting identity.
Finally, (51a)-(51b) follow from definitions (44a)-(44b) and the already proven identity (50a). □ Convention 3.1 (The fluid variables as functions of the Riemann invariants). In view of (50a)-(50b), throughout Sect. 3, we often view all the fluid variables as smooth functions of the Riemann invariants (or of only R PS (+) when we are studying simple isentropic plane-symmetric solutions), e.g., u 0 , u 1 , h, H = f(R PS (+) , R PS (−) ), where throughout the paper, f schematically denotes a smooth function that can vary from line to line. For example, we use (47), (50b), the inverse function theorem, and the chain rule to compute the following identity, which holds for isentropic plane-symmetric solutions: where on RHS (52), c ′ := dc dH . 3.5. Brief remarks on Riccati-type blowup and genuine nonlinearity. Consider a simple isentropic plane-symmetric solution R PS (+) to (49a) (with R PS (−) = 0). By differentiating the equation with respect to ∂ 1 , one obtains a Riccati-type equation for L 1 ̸ = 0 holds, then for solutions R PS (+) launched by smooth compactly supported initial data, one can easily show that ∂ 1 R PS (+) typically blows up along the integral curves of L, much like in the case of Burgers' equation ∂ t Ψ + Ψ∂ 1 Ψ = 0; we refer to [39, Section 7.5] for a detailed discussion of the notion of genuine nonlinearity in the context of 1D hyperbolic conservation laws. One can check that the non-degeneracy assumption (43) , that equation (49a) is genuinely nonlinear near (R PS (+) , R PS (−) ) = (0, 0). Our study of the blowup in Theorem 3.1 does not rely on differentiating the equations with ∂ 1 to obtain a Riccati-type PDE. Instead, it relies on a much sharper approach tied to nonlinear geometric optics, which we set up in Sect. 3.6. In 3D, such a sharpened approach seems essential, as there is no known simple approach for studying gradient-blowup based only on differentiating the equations with respect to a standard, solutionindependent partial derivative operator. In particular, in our discussion of shock formation in 3D in Sect. 6, nonlinear geometric optics will play a fundamental role.
3.6. Simple isentropic plane symmetry and nonlinear geometric optics. Our study of the h-MGHD in Theorem 3.1 intimately relies on nonlinear geometric optics, implemented via an acoustic eikonal function U , and a careful study of the corresponding acoustic geometry. The main idea of the proof of Theorem 3.1 is to show that the solution remains smooth relative to the "geometric coordinates" (t, U ) all the way up to the shock and to recover the formation of the shock as a degeneracy between the geometric coordinates and the Minkowski-rectangular coordinates (t, x 1 ).
In this section, we construct the eikonal function and corresponding acoustic geometry for the isentropic relativistic Euler equations in 1D, i.e., to the following initial value problem: More precisely, we restrict our attention to simple isentropic plane-symmetric solutions, defined here to be those solutions such that R PS (−) ≡ 0. From (53) and (44b), it follows that such solutions arise from initial data withR PS (−) ≡ 0. In this case, (50a) yields the identity u 0 = cosh( 1 2 R PS (+) ). As we alluded to in Convention 3.1, (50a)-(50b) also yield similar formulas for u 1 , H as functions of R PS (+) . We start by defining the acoustic eikonal function, which plays a central role in all of the forthcoming analysis.
Definition 3.3 (The acoustic eikonal function in 1D and the characteristics). Let h −1 be the inverse acoustical metric defined in (46). Then we define the acoustic eikonal function (eikonal function for short) to be the solution U to the following fully nonlinear hyperbolic PDE with Cauchy data: We refer to the level sets of U as "the characteristics." Remark 3.1 (Geometric optics is easier in 1D). Although understanding the fully nonlinear structure of (54a) is necessary in multiple space dimensions, in 1D, (54a) is equivalent to the linear (in U ) transport equation of U derived in the following lemma. This makes the analysis of the acoustic geometry in this section considerably simpler than in multi-D. It is also unique to 1D that the h-null vectorfield L can be explicitly expressed in terms of the fluid via the formula (44a). The correct analog of (44a) in higher dimensions is precisely Proof. Suppose U solves (54a)-(54b). Then by (46), we find that 0 = (LU ) · LU . If LU = 0, there is nothing to prove. Hence, we assume that LU ̸ = 0, and we will show that this leads to a contradiction. Then by (46), we have LU = 0. By (9), (51b), (54c), (55b), and our smallness assumption on R PS (+) , we have: (56) contradicts (54b). Conversely, if (55a)-(55b) hold, then since (46) implies that (h −1 ) αβ ∂ α U ∂ β U = −(LU ) · LU , it immediately follows that (54a) and (54c) hold. To prove (54b), we use (9), (51a), (55b), and the smallness of R PS (+) , arguing as in the proof of (56) to deduce that ∂ t U > 0, as is desired. □ Some remarks are in order.
• The minus sign in (55b) is just a convention imposed for consistency with related results that have appeared in the literature, such as [1,65,67,95]. • The precise initial condition in (55b) is not important; small perturbations of this initial condition could also be used to study the h-MGHD. • We have adapted our eikonal function to L because in Theorem 3.1, we will study initial data that lead to shock-forming solutions such that the singularity is tied to the infinite density of the Lcharacteristics and the blowup of ∂ 1 R PS (+) . There is nothing special about L compared to L, and it is only for convenience that we have adapted our constructions to it. That is, using the methods of this section, one could also construct an eikonal function U adapted to L and study different initial data that lead to a singularity tied to the infinite density of the L-characteristics and the blowup of One could also try to study richer problems in which singularities occur along both families of characteristics. We now define µ, the inverse foliation density of the characteristics. The name is justified by the fact that 1 µ is a measure of the density of the stacking of the characteristics relative to the level-sets Σ t ′ := {t = t ′ }. Roughly speaking, µ is initially positive, while shock formation corresponds to µ ↓ 0, which signifies the infinite density of the integral curves of L (or equivalently, the level sets of U in (t, x 1 )-space); we refer to 5 The solutions from Theorem 3.1 satisfy the needed smallness. The key point is that if R PS (+) is large, then L and L could both be right-pointing or both be left-pointing, in which case the initial condition (55b) would not be sufficient for identifying the "correct root of the eikonal equation," namely the "root" LU = 0. Fig. 3, a context in which µ vanishes along the curve portion denoted by "Υ(B [−U ,0] )" (we will explain these issues in much more detail later on).
Definition 3.4 (Inverse foliation density). We define the inverse foliation density of the characteristics, denoted by µ, as follows: In Sect. 6, which concerns solutions in 3D, we will define µ to be − 1 (h −1 ) κλ ∂κt∂ λ U ; see (157). Under the initial condition (156) that we will choose for U , that definition of µ can be shown to agree with the one for plane-symmetric solutions stated in (57).
We also note the following identities, which follow from (55b) and (57): We now define a collection of vectorfields associated to the eikonal function.
Definition 3.5 (The vectorfields X,X,L). Let L and L be the vectorfields from Def. 3.1 (they are h-null by (45)). We define: Using (44a)-(44b), (45), and (59), it is straightforward to see that X is Σ t -tangent (i.e., Xt = 0) and satisfies h(X, X) = 1, that h(X,X) = µ 2 , and that h(L,L) = 0. We also note that in plane-symmetry, using Def. 3.1, (42), and Def. 3.5, one can compute the following identity: Our main results concern solutions in which the factor c n on RHS (60) satisfies c n ≈ 1. Definition 3.6 (Geometric coordinates in 1D and the corresponding partial derivative vectorfields). We define (t, U ) to be the geometric coordinates, and we denote the corresponding geometric coordinate partial derivative vectorfields by ∂ ∂t , ∂ ∂U (which are not to be confused with the partial derivatives {∂ t , ∂ 1 } in the (t, x 1 )-coordinate system).
The following sets play a fundamental role in our analysis of solutions.
Definition 3.7 (Portions of submanifolds of geometric coordinate space R t × R U ). Given intervals I, J ⊂ R, relative to the geometric coordinates, we define: We often write P U instead of P . We often refer to the P U as the characteristics.
The following lemma shows that the geometric coordinates degenerate with respect to the Minkowskirectangular coordinates precisely when µ = 0. However, Theorem 3.1 shows that for appropriate initial data, Υ will still be a homeomorphism on the subset of the boundary of the h-MGHD where µ vanishes.
Lemma 3.4 (The Jacobian determinant of (t, U ) → (t, x 1 )). Let: denote the change of variables map from geometric coordinates to Minkowski-rectangular coordinates. Then the following identity holds: Proof. (67) follows from straightforward calculations based on the identity ∂ ∂U x 1 = − c n µ, which is implied by (59)-(60) and (62).  (62), we see that in simple isentropic planesymmetry, the transport equation (49a) takes the following form in geometric coordinates: From (68) and (53), we see that in geometric coordinates, the solution to (68) is: 3.6.2. The evolution equation for µ, the non-degeneracy condition, and explicit expressions for various solution variables. In this section, we derive various identities involving µ. In particular, we will later use the identities to show that suitable assumptions on the initial data of R PS (+) lead to the vanishing of µ in finite time, i.e., the formation of a shock.
We begin with the following lemma, which provides the evolution equation for µ. Moreover, Proof. First, using (51a), (52), (60), the quotient rule, the identity cosh 2 (z) − sinh 2 (z) = 1, and straightforward calculations, we deduce the following commutator identity, valid for simple isentropic plane-symmetric solutions: Equation (70) then follows from multiplying (57) by c/n, differentiating both sides with respect to L, and using the following commutator identity, which is a consequence the eikonal equation LU = 0: Equation (72) then follows from differentiating (70) with respect to L and using (49a), (60), (63), and the assumption of simple isentropic plane-symmetry. □ We find it convenient to work with an anti-derivative of the nonlinearity G defined in (71). To this end, viewing the quantity − as follows: From (62), (73), the fundamental theorem of calculus, and the chain rule, it follows that: We highlight the important fact that in simple isentropic plane-symmetry, the non-degeneracy assumption (43), the initial condition of (47), (73), and (74b) collectively imply that: From (73) and (75), it follows that the function R PS (+) → A[R PS (+) ] is a diffeomorphism from a neighborhood of the origin (i.e., near R PS (+) = 0) onto a neighborhood of the origin; we will silently use this basic fact in the rest of Sect. 3. As we mentioned in Sect. 3.5, (75) is fundamental for ensuring the formation of shocks, for it implies that near (R PS (+) , R PS (−) ) = (0, 0), equation (49a) is genuinely nonlinear. Next, we note that in simple isentropic plane-symmetry, in view of (69), the quantities G, c, and n can be viewed as functions of U alone. We sometimes emphasize this fact by using notation such as G(U ), c(U ), etc. The following corollary is a straightforward application of this fact and the above discussion; we omit the straightforward proof.  (73). Then for simple isentropic plane-symmetric solutions, the following identities hold relative to the geometric coordinates, where 3.7. Initial data that lead to admissible shock-forming solutions. In this section, we construct examples of initial data that fall under the scope of Theorem 3.1, i.e., initial data that launch shock-forming solutions such that we can derive the structure of a portion of the maximal classical development that includes a neighborhood of its boundary. Our construction could be substantially generalized to exhibit a much larger class of initial data for which the results of Theorem 3.1 hold; here we have striven for simplicity. Figure 1. The graph of a representative "seed profile," with U increasing from right to left.
3.7.1. Assumptions on the "seed" profile. To start, we fix any scalar "seed profile"φ with the following properties (it is straightforward to show that such functions exist): has a unique, non-degenerate minimum at U = 0 (our choice U = 0 is only for convenience).
In particular, d 3 dU 3φ (0) > 0. • Modifyingφ by multiplying it by a constant and composing it with a linear map of the form U → zU for some constant z ∈ R if necessary, we assume that (for the modifiedφ), there is a constantδ * satisfying:δ * > 0 (80) such that: such that: and such that: and that there is a constant p satisfying: such that: See Fig. 1 for the graph of a representative seed profile.

3.7.2.
The initial data for R PS (+) and amplitude smallness assumption. We now construct the initial data for R PS (+) . Letφ be as in Sect. 3.7.1. We define the initial data functionR PS (+) := R PS (+) | t=0 to be the following function of U :R where (73).
Amplitude smallness assumption on the initial data By (75) and continuity, we can further modifyφ by multiplying it by a constant to shrink the amplitude of R PS (+) and ensure that the following two conditions hold: (1) The isentropic plane-symmetric initial data arrayV : where H is the regime of hyperbolicity defined in (27).
(2) The following non-degeneracy condition holds: In the rest of Sect. 3, we will use the above two assumptions without always explicitly mentioning them. Taking into account (43) and (74b), we deduce from Taylor expansions and standard calculus that the following conclusions hold, where C > 0 depends onφ, and for functions and satisfies the following bounds: ] has a unique, negative, non-degenerate minimum at U = 0. In particular, d 3 ) ≤ C and such that: • For j = 1, 2, 3, there exists a continuous function Λ j : and such that for M = 0, 1, 2, we have: where p ∈ [0, 1) is the non-negative constant in (86), and in (91) and throughout, A ≈ B means that there is a C ≥ 1 such that 1 C B ≤ A ≤ CB. • The constantδ * appearing in (81) satisfies: where [z] − := max{−z, 0} is the negative part of z.

The time of first blowup.
Definition 3.8 (Time of first blowup). Given any initial data functionR PS (+) (U ) constructed above, with δ * > 0 as in (81), we define: Note that (77) and (93) imply that T Shock is the Minkowskian time at which µ first vanishes, which by (79) is the positive time of first blowup for |∂ 1 R PS (+) |. 3.7.4. Sharp estimates for the inverse foliation density. In Lemma 3.7, we provide various estimates for µ and its derivatives in a neighborhood of [0, T Shock ] × {0} ⊂ R t × R U . These estimates will be useful in the description of the boundary of the maximal development.
Note that by (9) and (17c), we have: Then there exists 7 a constant U depending only on the seed profileφ and satisfying 0 < U < 1 such that the following hold.
Behavior of µ on [−U , U ]. Let T Shock > 0 be as defined in (94). Then the following estimates hold for t ∈ [0, T Shock ]: Moreover, for (t, U ) ∈ [0, T Shock ] × [−U , U ], the following estimates hold, where here and throughout, A = O pos (B) means that there exist constants C 1 , C 2 (depending onφ) 8 such that 0 < C 1 ≤ C 2 and C 1 B ≤ A ≤ C 2 B, and A = O(B) means that there exists a C > 0 (depending onφ) such that |A| ≤ C|B|: In addition, we have: 7 With minor additional effort, one could explicitly compute such a constant U in terms of the constants associated toφ such asδ * , b, etc. 8 For example, the implicit constants present on RHS (99) can be taken to be C 1 = C −1 and C 2 = c −1 ; see (97).
µ PS is uniformly positive away from the interesting region. With p ∈ [0, 1) denoting the non-negative constant on RHS (92b), we have the following estimate: Properties of G in the interesting region. Let G be the source term from (71), viewed as a function of U , and let G ′ := d dU G. Then for U ∈ [−U , U ], the following estimates hold: In particular: Proof. All of the results except for those concerning G are straightforward consequences of Cor. 3.6, the properties of the initial dataR PS (+) (U ) from Sect. 3.7.2, and Taylor expansions; we omit the details. To prove the results for G, we first use (69) and (74a) Using this identity and (90), we conclude (107c), provided U > 0 is taken to be sufficiently small. □ 3.8. The crease, singular curve, and Cauchy horizon in 1D simple isentropic plane-symmetry.
In this section, we provide a description of a localized portion of the boundary of the h-MGHD. We emphasize the following key point, which follows from (69), (77), and (79): For the initial data from Sect. 3.7.2, the corresponding fluid solution R PS (+) and the inverse foliation density µ are explicit smooth functions on all 9 of R 2 in the (t, U ) differential structure, even though the Minkowskian partial derivative ∂ 1 R PS (+) blows up in finite time. Therefore, when we speak of the h-MGHD, we are referring to the standard differential structure corresponding to the (t, x 1 ) coordinates. Nonetheless, we prefer to carry out most our analysis of the h-MGHD in (t, U ) coordinates, thanks to the explicit structure of the equations and solutions in geometric coordinates revealed in Cor. 3.6. Hence, to reach conclusions about the behavior of solution in the (t, x 1 )-differential structure, we have to control the change of variables map Υ(t, U ) := (t, x 1 ); indeed, in Theorem 3.1, we reveal various diffeomorphism and homeomorphism properties of Υ.
For the solutions and localized region under study, when described in geometric coordinates, the boundary of the h-MGHD consists of two pieces, depicted in . We highlight again that the singular behavior is only visible relative to the differential structure of Minkowski space in (t, x 1 ) coordinates. The discrepancy between the behavior of the solution in the two coordinate systems is tied to the breakdown in the diffeomorphism property of Υ when µ vanishes; see (67).

The Relativistic Euler Equations
3.8.1. Singular curve, crease, Cauchy horizon, and developments. We now define some subsets of geometric coordinate space R t × R U that, in Theorem 3.1, will play an important role in describing the shape of the h-MGHD. Many of our definitions are localized to the subset {|U | ≤ U }, where U > 0 is the constant from Lemma 3.7. This is mainly for convenience, since Lemma 3.7 provides sharp control over the solution when |U | ≤ U .  (71), viewed as a function of U . We define: . (108) We define the truncated singular curve to be the following subset of geometric coordinate space R t × R U : whereX is the vectorfield defined in (59), and the second identity in (109) follows from (77), (96), and (108).
We define the crease, which we denote by ∂ − B, to be the following subset of geometric coordinate space We therefore have the following alternate characterization of B I : We make the following remarks.
• One might wonder why in the singular curve definition (109), there is a restriction to regions in whichXµ ≤ 0. This is because the portion of {µ = 0} on which we formally haveXµ > 0 never actually dynamically develops in the h-MGHD; that "fictitious portion" is cut off by the Cauchy horizon before it has a chance to dynamically form; see also Remark 3.2.  • To further flesh out the previous point, we note that for the solutions we study in Theorem 3.1, the smallest value of t on B [−U ,0] will occur precisely at U = 0, and we will have t Sing (0) = T Shock .
Definition 3.10 (Truncated Cauchy horizon). Assume that ∂ − B = (T Shock , 0), where T Shock is defined by (94) (Theorem 3.1 shows that this assumption is satisfied for the solutions under study here). LetL be the vectorfield defined in (59). Letγ be the integral curve of 1 2L in geometric coordinate space R t × R U parametrized by U (note that 1 2L U = 1 by (64)), with the initial conditionγ(0) = ∂ − B. That is,γ solves the following initial value problem: We use the notation t CH to denote the t-component of the curveγ. That is, by (64), we have: where t CH solves the following initial value problem: Then for an interval I ⊂ [0, ∞) of non-negative 10 U -values, we define the truncated Cauchy horizon to be: We now define M * , which is a subset of R t × R U that is tailored to the shape of the h-MGHD.  (108) and (113). We define the following subsets of geometric coordinate space R t × R U : The set M * defined by (118) is the region of geometric coordinate space lying between Σ where we can justify the use of (t, U ) coordinates to reach conclusions about the behavior of the solution in (t, x 1 ) coordinates. 11 3.8.2. The main theorem. We now prove Theorem 3.1, which is the main result of Sect. 3. On the region M * , the theorem provides a sharp description of the solution in the differential structure of the geometric coordinates (t, U ) as well as the behavior of the solution on Υ(M * ) in differential structure of the Minkowskirectangular coordinates (t, x 1 ), where we recall that Υ(t, U ) = (t, x 1 ) is the change of variables map. In particular, the theorem reveals various diffeomorphism and homeomorphism properties of Υ and exhibits the blowup of |∂ 1 R PS (+) | as Υ(B [−U ,0] ) is approached by points in Υ(M * ). In Fig. 2, we have depicted the region in geometric coordinate space R×R U . In Fig. 3, we have depicted the region in Minkowski-rectangular coordinate space R × R x 1 , i.e., the image under Υ of the region from geometric coordinate space.  10 In Theorem 3.1, the Cauchy horizon will correspond to non-negative values of U . 11 With some additional effort, given the results of Theorem 3.1, we could extend the region of spacetime on which Υ(t, U ) = (t, x 1 ) is a homeomorphism to include the infinite region M Large * ∞)). We refrain from giving the proof because it would require some technical, but uninteresting, adjustments of the arguments presented in the proof of the theorem. The details are uninteresting in the sense that µ is uniformly positive on M Large   Classical existence with respect to the geometric coordinates. With respect to the geometric coordinates (t, U ), the following results hold.
• R PS (+) , the null vectorfields L andL, the vectorfieldX, and the inverse foliation density µ exist classically on • The following estimate holds: In particular, with M * denoting the closure of the region M * defined in (118) • The crease, which is defined in (110), is a single point: • The singular curve function t Sing = t Sing (U ) from (108) and (111) satisfies the following estimates: • The rectangular spatial coordinate x 1 satisfies the following estimate on B [−U ,0] : • The integral curve initial value problem forγ(U ) given by (112a)-(112b) has a unique solution on the interval U ∈ [0, U ]. Hence, the truncated Cauchy horizon C [0,U ] , defined in (115), can be parameterized as follows: as in (115). Moreover, t CH (U ) satisfies the following estimates: • The inverse foliation density µ, defined in (57), satisfies the following estimate along the Cauchy horizon portion C [0,U ] : • The rectangular spatial coordinate x 1 satisfies the following estimate along the Cauchy horizon portion • µ is positive on M * but vanishes on B [−U ,0] . In particular, (126) implies that µ is positive on C (0,U ] but vanishes on the crease ∂ − B ⊂ C [0,U ] , which corresponds to U = 0 in (126).
Description of the singularity formation and the h-MGHD in Minkowski-rectangular coordinates.
• Let Υ(t, U ) := (t, x 1 ) be the change of variables map from geometric coordinates to Minkowskirectangular coordinates. Then Υ is a diffeomorphism from M * onto its image. Moreover, Υ is a homeomorphism from M * onto its image.
• The solution R PS (+) exists classically relative to the Minkowski-rectangular coordinates (t, x 1 ) on the subset Υ(M * ) of R t × R x 1 .
With respect to the (t, The causal structure of the Cauchy horizon.
• C [0,U ] is a h-null curve portion, i.e., its tangent vector is null with respect to h.
The causal structure of the singular curve.
• Relative to the geometric coordinates (t, U ), the vectorfield Q defined by: is well-defined on the singular curve portion B [−U ,0) , and Q is tangent to B [−U ,0) . • On B [−U ,0) , we have Qt > 0 and QU < 0. In particular, Q is future-directed and transversal to the characteristics P U . • For every q ∈ B [−U ,0) , the pushforward vectorfield [dΥ(q)] · Q(q), which is tangent to the curve Since there exist integral curves of the h-null vectorfield L κ ∂ κ that foliate and are h-orthogonal to Υ(B [−U ,0) ), it follows that Υ(B [−U ,0) ) is a h-null curve in the Minkowskian coordinate differential structure.
Before proving the theorem, we make some remarks.
• All the important results of the theorem are stable under perturbations of the initial data of R PS (+) . That is, even though the theorem was stated only for initial data that satisfy the assumptions of Sect. 3.7.2, the proof shows that if one perturbs any of those data by an arbitrary (small) C 4 function, then the corresponding perturbed simple isentropic plane-wave solution has an h-MGHD/exhibits singularity formation that is quantitatively and qualitatively close to the un-perturbed solution.
• In a similar vein as the previous point, with only modest additional effort, the results of Theorem 3.1 could be extended to a large, open set of (non-simple) isentropic plane-symmetric solutions in which R PS (−) is sufficiently small compared to R PS (+) (or, alternatively, in which R PS (+) is sufficiently small compared to R PS (−) ). • The emergence of the Cauchy horizon from the crease is fundamentally tied to the fact that in 1D, the isentropic compressible Euler equations have two characteristic directions, L and L, as is evident from equations (49a)-(49b). In particular, even within the class of simple isentropic planesymmetric solutions in which R PS (−) (which is transported by L) vanishes, when investigating the structure of the h-MGHD, one cannot "ignore" the L direction. This is in contrast to Burgers' equation, in which there is only one characteristic direction. Relatedly, for Burgers' equation, there does not exist a notion of "globally hyperbolic" that enjoys the same fundamental significance as the notion of "globally hyperbolic" in the context of Lorentzian geometry, relativistic Euler flow, or wave equations.
• From the point of view of determinism, the Cauchy horizon is an ally in that it "blocks" the nonuniqueness of classical solutions (in h-globally hyperbolic regions). In particular, in Fig. 3, one can see that the Cauchy horizon prevents the characteristics in the left of the figure from entering into the region on the right, where the fluid's gradient blows up. That is, the Cauchy horizon prevents crossing of the characteristics Υ(M * ), which in turn prevents multi-valuedness/non-uniqueness of classical solutions in this region. This should be contrasted against Burgers' equation, which exhibits nonuniqueness of classical solutions due to multi-valuedness stemming from the crossing of characteristics emanating from widely separated regions in the initial data hypersurface. ), which in particular shows that L α does not enjoy the Lipschitz regularity that lies behind standard ODE uniqueness theorems for integral curves. • One can show that in 1D, the acoustical metric takes the following form in geometric coordinates: h = −µdt⊗dU −µdU ⊗dt+µ 2 dU ⊗dU . In particular, in geometric coordinates, h vanishes at points such that µ = 0. Hence, h is not capable of "measuring" causality along the singular curve portion B [−U ,0] , along which µ vanishes. We note, however, that in more than one space dimension, even when µ = 0, there are some remaining non-degenerate components of h in geometric coordinates. Relatedly, in our 3D paper [1] on the non-relativistic compressible Euler equations, we prove that the ∂ − B is a co-dimension 2 h-spacelike submanifold of geometric coordinate space.
Proof of Theorem 3.1. Throughout the proof, we will silently shrink the size of U as necessary; this will in particular guarantee that some important terms arising in various Taylor expansions have a sign.
Classical existence and properties of the solution relative to the geometric coordinates: The fact that R PS (+) (t, U ) =R PS (+) (U ) is smooth on all of R t × R U was proved in (69). In particular, in geometric coordinates, R PS (+) smoothly extends to the closure M * . The fact that R PS (+) vanishes on the complement of {(t, U ) ∈ R 2 | U ∈ [−U 1 , U 2 ]} then follows from our assumptions on the initial data, which in particular imply thatR PS (+) is supported in [−U 1 , U 2 ]. Next, using the identity (77), the estimate (96), and the fact that R PS (+) is supported in {(t, U ) ∈ R 2 | U ∈ [−U 1 , U 2 ]}, we conclude the asserted properties of µ in geometric coordinates, except for the positivity of µ on M * , which we prove below. The asserted properties of L,L, Lx 1 , and Lx 1 in geometric coordinates then follow from (51a)-(51b), (62), (63), and the fact that The estimate (119) follows from the first equality in (76) and the estimate (107a). The estimates (122a)-(122b) follow from definitions (94) and (108) and the estimates (107b)-(107c).
Proof of the properties of C [0,U ] relative to the geometric coordinates: We first prove that the integral curve initial value problem forγ(U ) given by (112a)-(112b) (or equivalently, (114)) has a unique solution on the interval U ∈ [0, U ]. Since µ is a smooth function on all of (t, U ) ∈ R 2 , the standard existence and uniqueness theorem for ODEs yields a unique solution to (114) for U ∈ [0, U ], provided we shrink the size of U if necessary.
Finally, C [0,U ] is a h-null curve portion because its tangent vector is everywhere proportional to the vectorfieldL, which is h-null by (45) and (59).
Proof that µ > 0 on M * : We have already shown that µ > 0 on M Sing . Hence, in view of definition 118, we see that to prove µ > 0 on M * , it suffices to prove µ > 0 on M Reg . To this end, we first note that the same arguments given above imply that the curve U → (t Sing (U ), U ) is smooth for U ∈ [−U , U ] and that t Sing (U ) satisfies the estimate stated in (122b) for U ∈ [−U , U ] (not just for U ∈ [−U , 0]). The portion of this curve corresponding to U ∈ (0, U ] is shown as a dotted curve in Fig. 2; for purposes of this proof, we denote this dotted curve, which by definition does not contain the crease, by "B Fictitious " ("fictitious" because the image of B Fictitious under Υ in (t, x 1 )-space is not part of the h-MGHD). The same arguments that we used to prove that µ > 0 on M Sing (which relied on (77), (96), and (108)) imply that µ > 0 on D := (t, U ) ∈ R 2 | 0 ≤ t < t Sing (U ), U ∈ (0, U ] , where we note that the top boundary of this region is B Fictitious . The key point is that (113), (115), and the estimates (122b) (valid also for U ∈ (0, U ]) and (125b) imply that t Sing (U ) < t CH (U ) for U ∈ (0, U ] and thus C (0,U ] lies below B Fictitious . In view of definition (117), we see that M * ⊂ D and hence µ > 0 on M Reg , as desired.
Proof of the diffeomorphism and homeomorphism properties of Υ: From Def. 3.11, we see that the top boundary of M * is the following U -parameterized curve (note that the curve is well-defined since t Sing (0) = t CH (0) = T Shock ): Using (123), (127), and the mean value theorem, we see that x 1 (t top (U ), U ) is strictly decreasing for U ∈ [−U , U ], i.e., the x 1 coordinate function strictly decreases as we move from right to left along the top boundary of M * in Fig. 2. Next, we note that (59) and (60) imply that ∂ ∂U x 1 = − c n µ. Using this identity, the fact that c n > 0, and the fact that we have already shown that µ > 0 on M * and µ vanishes precisely along the curve portion B [−U ,0] in M * , we see that for any (t, U ) ∈ M * , the map U → x 1 (t, U ) is strictly decreasing. From these two monotonicity properties of x 1 , we conclude, given the shape of M * (see Fig. 2), that the map Υ(t, U ) = (t, x 1 (t, U )) is injective on the compact set M * and thus is a homeomorphism from M * onto Υ(M * ). Also considering the identity det dΥ = − c n µ proved in (67), we further conclude that Υ is a diffeomorphism on M * , where µ > 0.
Proof of (128) and the structure of the h-MGHD in Minkowski-rectangular coordinates (t, x 1 ): We deduce from (59), (60), (69), (76), the chain rule, the non-degeneracy assumption (88), and (119) that for (t, U ) ∈ M * , the following estimate holds: From (134), we conclude (128). The continuity of R PS (+) , L M R PS (+) and h αβ = h αβ (R PS (+) ) all the way up to Υ(M * ) in Minkowski-rectangular coordinates (t, x 1 ) follows from the fact that R PS (+) (t, U ) is smooth on R t × R U and the already proven fact that Υ is a homeomorphism from M * onto Υ(M * ). The non-degeneracy of the 2 × 2 matrix h αβ relative to the (t, x 1 )-coordinates follows from the expression (40a) for h αβ = h αβ (R PS (+) ), our assumed smallness on the amplitude of the initial data (which in particular implies that 0 < c ≤ 1), and the identities stated in Prop. 3.2.
Proof of the causal structure of the singular curve in Minkowski-rectangular coordinates (t, x 1 ): Using (62), (76)- (77), and (109), we deduce that along B [−U ,0) , the vectorfield Q defined by (129) can be expressed as follows, where G is the source term from (71), viewed as a function of U , and G ′ is its derivative: Hence, using (107a) and (107c), we see that Q is well-defined on B [−U ,0) , i.e., the denominator term G ′ (U ) on RHS (135) is non-zero for U ∈ [−U , 0). Next, we deduce from (129) that Q c n µ = 0, i.e., that Q is tangent to the level sets of c n µ. Using the fact that c n > 0 (see (96)), we see that the zero level set of c n µ coincides with the zero level set of µ. Hence, in view of definition (109), we see that Q is tangent to B [−U ,0) . The inequalities Qt > 0 and QU < 0 follow easily from (135), (107a), and (107c).

The new formulation of the flow in three spatial dimensions
In Sect. 6, we provide an overview of how to extend some aspects of the simple isentropic plane-symmetric shock formation results yielded Theorem 3.1 to the much harder case of three spatial dimensions in the presence of vorticity and entropy. The strategy we present in Sect. 6 is based on the approach we used in studying shock formation for the non-relativistic multi-dimensional compressible Euler equations [1,65,67]. Fundamental to that strategy is the availability of a new formulation of the flow as as system of geometric wave-transport-div-curl equations that exhibit remarkable regularity properties and geometric null structures. As we explain in Sect. 6, such a formulation allows one to implement multi-dimensional nonlinear geometric optics, which is important for the study of shocks (as is already evident from Theorem 3.1). In the non-relativistic case, such a formulation was derived in [66,94], while in the relativistic case, such a formulation was derived in [41]. In Sect. 4, in the relativistic case, we set up the machinery needed to state the formulation from [41] and present it in a simplified, schematic form as Theorem 4.1. In Sect. 4.2.5, we make some brief comments on how the new formulation can be extended to more general spacetimes (M, g).

4.1.
Connection to quasilinear wave equations. Although our study of isentropic plane-symmetric solutions in Sect. 4 was based on analyzing the first-order Riemann invariant system (49a)-(49b), as we will see in Sect. 6, in multi-dimensions, it is advantageous to work with the formulation of the flow provided by Theorem 4.1, which features geometric wave equations. The main advantage is that for geometric wave equations, there is an advanced geo-analytic machinery available, tied to the vectorfield method, for studying the global behavior of solutions. To motivate the relevance of wave equations for the study of relativistic fluids, we now establish a simple connection between quasilinear wave equations and the Riemann invariant formulation of isentropic plane-symmetric solutions that we provided in Sect. 4. More precisely, upon differentiating (49a)-(49b), we deduce that for isentropic plane-symmetric solutions, the Riemann invariants (R PS (+) , R PS (−) ) satisfy the following coupled system of quasilinear wave equations: Definition 4.1 (The u-orthogonal vorticity of a one-form). Let ξ be a one-form. We define its u-orthogonal vorticity, denoted by vort α (ξ), to be the following vectorfield, where ϵ is the fully antisymmetric symbol normalized by ϵ 0123 = 1: and we recall that throughout the paper, we use the Minkowski metric m -as opposed to h -to lower and raise indices. In particular, ϵ 0123 = −1.
Definition 4.2 (Vorticity vectorfield). We define the vorticity ϖ α to be the following vectorfield: Definition 4.3 (Entropy gradient one-form). We define S α to be the following one-form: The fluid variables C and D from the next definition play a fundamental role in the setup of the new formulation of the flow. In particular, Theorem 4.1 shows that they satisfy transport equations with source terms that enjoy remarkable regularity and null structures. . We define C α and D to respectively be the following vectorfield and scalar function: 4.2.2. Null forms. In Def. 4.6, we recall the definition of standard null forms relative to h. The key point is that in Theorem 4.1, all derivative-quadratic terms on the RHS of the equations are null forms relative to the acoustical metric. In Sect. 6.8.3, we will provide a detailed explanation as to why, in the context of shock formation, such null forms are harmless error terms. However, at this point in the paper, we merely state that h-null forms are derivative-quadratic nonlinearities that, when a shock forms, are strictly weaker than the Riccati-type terms that drive the blowup (see Sect. 3.5). As we will further explain in Remark 6.4, the dangerous Riccati-type terms are "hiding" in the definition of the covariant wave operator on LHS (147a).
Definition 4.6 (Standard null forms relative to h). We define the standard null forms relative to h (which we refer to as "standard h-null forms" for short) to be the following derivative-quadratic terms, where ϕ and ψ are scalar functions and 0 ≤ µ < ν ≤ 3:

Two vectorfields.
The two vectorfields featured in the following definition will play a role in our forthcoming discussion.
Definition 4.7 (B and N). We define B and N to respectively be the vectorfields with the following components relative to the Minkowski-rectangular coordinates (t, x 1 , x 2 , x 3 ): Note that Bt = 1, while by (19), N is the future-directed h-unit normal to Σ t . In particular, We also compute that: and, with the help of (12d) and (17a), that: Note that RHS (146) < 0 because n > 0, B is future-directed and h-timelike, L is future-directed and h-null, and the h-null cones are inside the m-null cones (see Remark 2.3).

4.2.4.
The geometric wave-transport-div-curl formulation of the flow. We now state -in a somewhat abbreviated form -the formulation of relativistic Euler flow derived in [41]. The proof of theorem is based on differentiating the equations (12a)-(12d) in well-chosen directions and finding a myriad of cancellations and special structures.
4.2.5. Comments on the case of more general ambient spacetimes (M, g). Here we make some brief comments on some of the adjustments that would be needed to extend the results of Sect. 4 to general g-globally hyperbolic spacetimes (M, g). By g-globally hyperbolic, we mean that there exists a g-Cauchy hypersurface Σ 0 ⊂ M. By a "g-Cauchy hypersurface," we mean a co-dimension one submanifold Σ 0 such that every inextendible g-causal curve in M intersects Σ 0 exactly once. It is well known [17,45] that g-globally hyperbolic spacetimes admit a smooth time function t : M → R. By a "time function," we mean that t −1 (0) = Σ 0 , that the level sets Σ t ′ := {t ≡ t ′ } are g-spacelike Cauchy hypersurfaces, that M = ∪ t Σ t , and that t has a past-directed timelike gradient ∇ # t := g −1 · dt (which is everywhere normal to Σ t ). Here, ∇ # t denotes the vectorfield equal to the dual of ∇t with respect to g. First, we note that the extra terms described in Remark 2.1 would lead to the presence of extra terms in the system (147a)-(147d). However, as we described in Remark 2.1, in the context of the study of shock formation, in which fluid gradients are large, the extra terms would be small relative to the shock-driving Riccati-type fluid terms (which are "hidden" in □ h( ⃗ Ψ) Ψ; see Remark 6.4). Next, we note that under the above assumptions, we can define a lapse function Φ by Φ := (−g(∇t, ∇t)) −1/2 . We can then define the vectorfield T := −Φ 2 ∇ # t, which is g-normal to Σ t and satisfies T t = 1. Hence, given a coordinate system (x 1 , x 2 , x 3 ) on Σ 0 , it can be propagated to any Σ t by the flow of T , thereby yielding a coordinate system (t, x 1 , x 2 , x 3 ) on spacetime. In this setup, the ambient spacetime metric can be decomposed as follows: where g = g(t) is the Riemannian metric on Σ t induced by g, i.e., the first fundamental form of Σ t with respect to g. Given g and the relativistic fluid, the corresponding acoustical metric h is defined by (21). The upshot of introducing a time function t is that it allows for many of the constructions from this section to easily be generalized. For example, using (2) and (148), one sees that u 0 := u · ∇t := u α ∂ α t > 0 and hence one could naturally define the analog of (143a) by B := 1 u 0 u and the analog of (143b) by N := −h −1 · dt.

Some prior works on shocks
In this section, we discuss some prior works on shocks. There is a vast literature, and we cannot hope to mention all of it here. We have aimed to discuss works that provide further context and motivation for Sects. 2 and 4 and to help prepare the reader for Sects. 6 and 7.

5.1.
Comments on the 1D theory. In Riemann's foundational paper [83], he developed the method of Riemann invariants and combined it with the method of characteristics to prove that for the 1D nonrelativistic compressible Euler equations, there exist large sets of initial data such that a shock forms in the sense that the fluid's gradient blows up in finite time, though the fluid variables themselves remain bounded. Once one has formulated the flow in terms of Riemann invariants, the proof of the blowup is based on differentiating the equations to obtain a coupled system of Riccati-type equations along the characteristics. This can be viewed as a multi-directional (i.e., there are two characteristic directions in 1D isentropic compressible Euler flow) version of the blowup of ∂ x Ψ that often occurs in solutions to Burgers' equation ∂ t Ψ(t, x) + Ψ∂ x Ψ(t, x) = 0; in Sect. 3.5, we discussed this approach in more detail in the case of simple isentropic plane-symmetric solutions to the relativistic Euler equations. Lax later generalized [62] Riemann's results to 2 × 2 genuinely nonlinear hyperbolic systems in 1D. John later extended [53] Lax's blowup-results to apply to some systems in 1D with more than two unknowns. Christodoulou and Raoul-Perez used a sharpened, more geometric version (in the style of our proof of Theorem 3.1) of Lax's approach to give a detailed proof of shock formation in 1D for electromagnetic waves in nonlinear crystals [35]. Related results have recently been proved for plane-symmetric solutions to the equations of elasticity [12] and the equations of non-relativistic magnetohydrodynamics [13], and in these works, the plane-symmetric blowup was also used to prove ill-posedness results in 3D; see also the survey article [14]. For the relativistic Euler equations in 1D, stable blowup-results have been proved for large sets of initial data (which are allowed to be large) for isentropic flows [16] as well as flows with dynamic entropy [15]. Theorem 3.1 provides a sharpened version of these kind of results for simple isentropic plane-symmetric solutions to the relativistic Euler equations; by "sharpened," we mean that, different from the works mentioned above, Theorem 3.1 follows the solution all the way up to the singular boundary and Cauchy horizon, rather than just up to the time of first blowup.
Even though shocks can form for large classes of hyperbolic PDEs, signifying the end of the classical evolution, ideally, one would like to develop a theory of weak solutions that can accommodate the formation of shocks and their subsequent interactions. In 1D, a robust such theory exists, at least for strictly hyperbolic systems and initial data that are small in a suitable bounded variation (BV) function space. We refer to [39] for a comprehensive discussion of the 1D theory. In multi-dimensions, much less is known. A key obstacle is that it is known, thanks to the important paper [80] by Rauch, that quasilinear hyperbolic PDEs are typically ill-posed in BV spaces. In fact, for most multi-dimensional quasilinear hyperbolic systems of interest, even local well-posedness is known only in L 2 -type Sobolev spaces. In practice, this means that to prove shock formation in multi-dimensions, one must commute the equations and derive energy estimates up the singularity; this turns out to be a difficult task, for reasons we discuss in detail in Sect. 6.

5.2.
Blowup-results in multi-dimensions without symmetry: proofs by contradiction. In light of the previous paragraph, it is perhaps not surprising that the first blowup-results without symmetry assumptions in fluid mechanics were non-constructive. Specifically, in [85], Sideris studied the 3D compressible Euler equations under equations of state that satisfy a convexity assumption. For open sets of initial data (including near-constant-state-data) without symmetry or irrotationality assumptions, he used arguments based on integrated quantities to give a proof by contradiction that a smooth solution cannot exist for all time. For the 3D relativistic Euler equations and relativistic Euler-Maxwell equations, Guo and Tahvildar-Zadeh proved similar results [48] for open sets of initial data that are large perturbations of constant-states with positive density. In contrast to the present work, [48,85] relied on convexity assumptions on the equation of state.

5.3.
Proofs of shock formation in multi-dimensions via nonlinear geometric optics. Alinhac [8][9][10][11] was the first to prove stable shock formation without symmetry assumptions for scalar quasilinear wave equations in 2D and 3D of the form: whenever the nonlinearities in (149) fail to satisfy the null condition. In (149), h is a Lorentzian metric whose components h αβ in the standard coordinate system are prescribed functions of ∂ ∂ ∂Φ, i.e., functions of the spacetime-gradient of Φ in the standard coordinate system. Alinhac's proof applied to open sets of initial data that satisfy a non-degeneracy assumption that guaranteed that within the constant-time hypersurface of first blowup, the singularity is an isolated point. His non-degeneracy assumption guaranteed that the singular boundary B and its past boundary ∂ − B are strictly convex, as in Fig. 4, and in the context of the figure, his approach allowed him to follow the solution up to the flat constant-time hypersurface that contains the point b lowest ; in Sect. 5.5, we discuss B and ∂ − B in more detail. Alinhac's main results showed that at b lowest , a shock forms in the sense that ∂ ∂ ∂ 2 Φ blows up while ∂ ∂ ∂ ≤1 Φ remains bounded. His proof relied on nonlinear geometric optics (i.e., eikonal functions and the method of characteristics, as in Sects. 3.6 and 6.3), and his proof showed that relative to a geometric coordinate system constructed out of the eikonal function, Φ and its partial derivatives remain bounded, except possibly at the high derivative levels. To close the energy estimates in the geometric coordinate system without derivative loss, he relied on a Nash-Moser iteration scheme that necessarily terminated at the time of first blowup. In Christodoulou's breakthrough monograph [30], he proved a substantially sharper result for the 3D relativistic Euler equations in regions where the solution is irrotational and isentropic. The dynamics in such regions is described by a scalar quasilinear wave equation for a potential function Φ (i.e., ∂ ∂ ∂Φ determines the four-velocity and proper energy density) that is of the form (149) and that is invariant under the Poincaré group. Like Alinhac, Christodoulou used an eikonal function in his proof, that is, a solution U to the eikonal equation: The level sets of U are characteristic for equation (149), i.e., hypersurfaces that are null with respect to h. Also like Alinhac, Christodoulou constructed a geometric coordinate system tied to U and proved that relative to the geometric coordinates, the solution remains smooth, except possibly at the high derivative levels. He showed that for open sets of initial data, the shock singularity -i.e., the blowing-up of the secondorder Minkowski-rectangular derivatives ∂ ∂ ∂ 2 Φ while ∂ ∂ ∂ ≤1 Φ remains bounded -is tied to the vanishing of the inverse foliation density µ, defined by: where t is the standard time coordinate on Minkowski spacetime. The vanishing of µ signifies the infinite density of the level-sets of U in Minkowski-spacetime, i.e., the infinite density of the characteristics and the blowup of ∂ ∂ ∂U , as in Theorem 3.1. Motivated by techniques used in his joint proof of the stability of Minkowski spacetime with Klainerman [32], Christodoulou was able to avoid Nash-Moser estimates. Rather, to avoid derivative loss in the eikonal function and to close the proof, he used a geometric energy method for wave equations and U ; see Sect. 6 for further discussion of a related geometric energy method in the context of the 3D relativistic Euler equations with vorticity and entropy. The results of [30] go beyond Alinhac's in several crucial ways. First, for smooth compactly supported perturbations of non-vacuum constant-state initial data, Christodoulou proved that the solution remains global unless µ vanishes, i.e., he proved a conditional global existence result. Second, his geometric energy method allowed for commutator and multiplier vectorfields with time and radial weights, which allowed him to simultaneously study the interaction between the dispersive tendency of waves and the fact that nonlinearities can focus waves and cause singularities. Third, for open sets of initial data that do not have to satisfy Alinhac's non-degeneracy condition, Christodoulou proved that µ vanishes in finite time, leading to the blowup of ∂ ∂ ∂ 2 Φ. In particular, even without strict convexity of the type shown in Fig. 4, Christodoulou was able show that there is at least one singular point. Finally, for open sets of initial data that satisfy a non-degeneracy assumption that is weaker than Alinhac's, Christodoulou was able to follow the solution beyond the time of first blowup and to reveal a portion of the h-MGHD, up to the boundary. The portion of the boundary of the h-MGHD revealed by [30] was not described explicitly and thus we cannot make straightforward comparisons with Theorem 3.1. The weaker notion of non-degeneracy used by Christodoulou coincides with our notion of transversal convexity, which we describe in Sect. 5.5.
Christodoulou's results [30] have been extended and applied to other equations, including: • Under the assumption of an irrotational and isentropic flow, the (non-relativistic) 3D compressible Euler equations were handled in [34]. • A larger class of wave equations (e.g., without the assumption of Poincaré invariance) was treated in [75,91]. See also the survey article [49], which is centered around the works [30,91]. • Solution regimes that are different than the small, compactly supported data regime were treated in [49,75]. • Solutions that exist classically precisely on a past-infinite half-slab were studied in [74].
• Other multiple speed systems in multi-dimensions were treated in [92,93].
• 2D compressible Euler solutions with vorticity were handled in [65] with the help of the new formulation of the flow derived in [66], which yields a non-relativistic analog of Theorem 4.1 for barotropic equations of state, which by definition are such that the pressure is a function of the density. • 3D compressible Euler solutions with vorticity and entropy were treated in [67] with the help of the new formulation of the flow derived in [94], which extends the results of [66] to general equations of state in which the pressure is a function of the density and entropy. • For the 3D compressible Euler equations with vorticity and entropy, a complete description of open sets of solutions up to a neighborhood of the boundary of the h-MGHD is provided by our series [1][2][3], which collectively yield a multi-dimensional version of Theorem 3.1 for the non-relativistic equations. See also Sect. 5.5.

5.4.
A different approach to multi-dimensional singularity formation. In the works [22,23], the authors developed a new approach for proving gradient-blowup for for solutions the 3D compressible Euler equations under adiabatic equations of state with vorticity (and entropy in [23]). Instead of using nonlinear geometric optics, the authors used modulation parameters to show that for open sets of initial data with large gradients, the solution forms a gradient singularity in finite time that is a perturbation of a self-similar Burgers-type shock. The approach allows one to follow the solution to the time of first blowup, but not further. It applies to non-degenerate initial data such that within the constant-time hypersurface of first blowup, the singularity occurs at an isolated point. That is, the approach applies when the singular boundary B and its past boundary ∂ − B are strictly convex, as in Fig. 4, and in the context of the figure, it allows one to follow the solution up to the flat constant-time hypersurface that contains the point b lowest , where the singularity occurs. Such singularities can be viewed as fluid analogs of the ones studied by Alinhac in his aforementioned works [8][9][10][11] on quasilinear wave equations. See also the precursor work [24], in which the authors studied the same problem for the 2D compressible Euler equations in azimuthal symmetry, and the work [21], which, in the same symmetry class, constructs unstable self-similar solutions whose cusp-like behavior at the singularity is non-generic. Self-similar blowup is a phenomenon that occurs in many other PDEs besides those of fluid mechanics. In particular, singularity formation modeled on a self-similar Burgers'-type shock has been proved for various non-hyperbolic PDEs [36-38, 61, 76, 79, 101].

5.5.
Maximal classical h-globally hyperbolic developments. It is of mathematical and physical interest to study shock-forming solutions in a region of classical existence that goes beyond the constant-time hypersurface of first blowup, e.g., in a region larger than the one studied by Alinhac in [8][9][10][11], as we described in Sect. 5.3. The holy grail object in this vein is the maximal (classical) h-globally hyperbolic development (h-MGHD), which is what we studied in Theorem 3.1 within the class of simple isentropic plane-symmetric solutions to the relativistic Euler equations. Recall that, roughly speaking, the h-MGHD the largest possible classical solution + h-globally hyperbolic region that is uniquely determined by the initial data. In particular, the h-MGHD can have a complicated boundary that includes points lying to the future of the constant-time hypersurface of first blowup, as in Figs. 3 and 6. Here, "h-globally hyperbolic" means that the development contains a Cauchy hypersurface, i.e., a surface such that every inextendible h-causal curve (where h is the acoustical metric) in the development intersects it. In the important work [42], the authors showed that for general quasilinear hyperbolic PDEs, one cannot ensure uniqueness of the MGHD until one constructs it and shows that it enjoys some crucial structural properties, which are global in nature.
In our series [1][2][3], for open sets of initial data for the 3D compressible Euler equations leading to shockforming solutions, we constructed a large (though bounded) portion of the h-MGHD, up to the boundary, i.e., in the non-relativistic case without symmetry assumptions, we proved an analog of Theorem 3.1. In particular, these are the first results to fully justify Fig. 6 for open sets of solutions to the 3D compressible Euler equations with vorticity and dynamic entropy. Although these works are of mathematical interest in themselves, they are also fundamental for setting up the shock development problem, as we describe in Sect. 5.6.
Our results [1][2][3] exhibit a localized version of the crucial property that the h-MGHD "lies on one side of its boundary," as is shown in Fig. 6. In [42], the authors showed, roughly, that if this "one-sidedness" property holds globally, then the h-MGHD is unique. They also showed, for a specific hyperbolic PDE and specific initial data, that uniqueness of the MGHD fails and the property fails! As in Theorem 3.1, the boundary portion we construct in [1][2][3] has two kinds of components: • A singular boundary B, constructed in [3], along which the fluid's gradient blows up, the past boundary of which is a 2-dimensional acoustically spacelike 12 submanifold known as the crease and denoted by "∂ − B" in Fig. 6. • A Cauchy horizon, constructed in [2] and denoted by "C" in Fig. 6, which is an acoustically null hypersurface that emanates from the crease and along which the solution extends smoothly (except at the crease). The crease plays the role of the true initial singularity in shock-forming solutions. Our analysis in [1,2] relies on an assumption of transversal convexity, which is a weak, stable form of convexity that can be read off the initial data and propagated by the flow. Roughly, transversal convexity means that when the inverse foliation density µ (which in the context of [1,2] is an analog of (151) or (157) for the 3D compressible Euler equations) is small, µ has a positive second derivative in a direction transversal to the level sets of the eikonal function U . For the simple isentropic plane-symmetric solutions we studied in Sect. 3, transversal convexity is ensured by the estimate (102) and is manifested in Fig. 3 by the upwards-bending nature of Υ(B [−U ,0] ). We highlight the importance of transversal convexity with the following remarks.
Even for simple isentropic plane-symmetric solutions, in the absence of transversal convexity, the qualitative structure of the singularity can be radically altered. For example, if instead of the estimate (102) there heldXXµ ≡ 0 at the first singularity, then in the (t, x 1 )-coordinates picture, an entire continuum of characteristics would fold into a single point, and the change of variables map Υ from Theorem 3.1 would dramatically fail to be an injection. This would be a serious obstacle to even the local well-posedness of the shock development problem, which we describe in Sect. 5.6. This is in contrast to the "favorable" situation shown in Fig. 3, in which, thanks to transversal convexity, the characteristics graze the singular boundary Υ(B [−U ,0] ), but distinct characteristics do not actually intersect along Υ(B [−U ,0] ). In the context of the 3D solutions depicted in Fig. 6, the transversal convexity is manifested by B being "convex in the x 1 -direction" (i.e., upwards bending in the x 1 direction) but not necessarily convex in the (x 2 , x 3 )-directions. Note that if we view plane-symmetric solutions as solutions in 3D with symmetry, then the corresponding singular boundary B is not strictly convex because of the "symmetric directions" (x 2 , x 3 ).
Of the works on shock formation cited in Sect. 5.3, Christodoulou's monograph [30] on 3D irrotational and isentropic relativistic fluids is the only one that follows the solution beyond the time of first blowup. He used eikonal functions and foliations of spacetime by h-spacelike planes Σ (more precisely, the Σ were coordinate planes with respect to the standard Minkowski-rectangular coordinates) to reveal an implicit portion of the boundary of the h-MGHD, where, for example, the portion of the crease ∂ − B that was revealed is such that the crease had to lie to the future of any Σ that is tangent to it. For initial data such that the crease is a strictly convex subset of Minkowski-rectangular coordinate space (as in Fig. 4), his approach can be used to study the entire crease and a neighborhood of the singular boundary that emerges from it. However, strict convexity does not hold for all shock-forming solutions. In particular, as we noted in the previous paragraph, strict convexity fails for symmetric solutions (roughly, strict convexity will fail in the directions of symmetry) and general small perturbations of them. In [1], under the weaker assumption of transversal convexity, we construct the entire crease and a neighborhood of the singular boundary that contains it for solutions with vorticity and dynamic entropy. To handle the lack of strict convexity, we dynamically construct a new foliation of spacetime by "rough acoustically spacelike hypersurfaces," depicted in Fig. 5, that are precisely adapted to the shape of the crease and singular boundary. We refer to the foliations as "rough" because they are less differentiable than the fluid solution, which is a key difficulty that has to be overcome in the PDE analysis. To handle the presence of vorticity and dynamic entropy, we combine the special structure of the geometric wave-transport-div-curl formulation of the flow derived in [94] with upgraded versions of the integral identities derived in [3], which in total yield a regularity theory for the vorticity and entropy that is adapted to the rough foliations (and hence the shape of the singular boundary) and that allows one to treat the fluid equations as a perturbation of a wave equation system. In particular, as in [65,67], the framework yields (and requires) one extra degree of differentiability for the vorticity and entropy compared to standard estimates. See also Sect. 6.9.2, in which we discuss the extra differentiability of the vorticity and entropy in the context of the 3D relativistic Euler equations. Figure 5. Rough foliations and the singular boundary from [1], depicted in Cartesian coordinate space 5.6. Shock development problem. In Majda's celebrated works [70,71], he solved the shock front problem for open sets of solutions to the non-relativistic 3D compressible Euler equations. Roughly, this means that he considered piecewise smooth initial data that have a jump discontinuity across a smooth hypersurface, and he then constructed a local weak solution to the equations and a shock hypersurface such that the solution is piecewise smooth on either side of shock hypersurface but has a jump discontinuity across it and satisfies the Rankine-Hugoniot condition and an entropy-type condition.
The shock development problem is closely related to the shock front problem but is laden with additional technical difficulties. It is the problem of describing how initially smooth solutions form a "first singularity" (which is the crease, as described in Sect. 5.5), from which emerges a shock hypersurface (which is not known in advance, as in the shock front problem), and describing the transition of the solution from classical to weak such that in the "weak solution region," the solution jumps across the shock hypersurface and satisfies the Rankine-Hugoniot jump condition and an entropy-type condition. In particular, in the shock development problem, there is no jump discontinuity in the initial data; rather, the jump emerges dynamically, starting from the crease. In problems that are effectively one-dimensional (and thus energy estimates can be avoided), the shock development problem has been solved for various hyperbolic PDEs [20,27,33,51,52,63].
The only solution to a multi-dimensional shock development problem without symmetry assumptions was provided by Christodoulou's recent breakthrough monograph [31], in which he used nonlinear geometric optics (i.e., a pair of acoustic eikonal functions, one ingoing and one outgoing) to solve the restricted shock development problem in an arbitrary number of spatial dimensions for the compressible Euler equations and the relativistic Euler equations. Roughly, the word "restricted" means that he solved an idealized problem in which he ignored the jump in entropy and vorticity across the shock hypersurface, thereby producing a weak solution to a hyperbolic PDE system that approximates the true one; in the real problem, which remains unsolved, the vorticity and entropy must jump across the shock hypersurface. Compared to the 1D case, the main new difficulty in multiple spatial dimensions is that one must derive energy estimates, which are degenerate near the first singularity for reasons related to the energy degeneracies that occur in the shock formation problem, as we discuss in Sect. 6.9.3.
We highlight the following key issue, which is a primary motivating factor for our construction of large portions of the h-MGHD in [1][2][3]: The setup of the shock development problem (see [31,33]) is such that the data for it are an h-MGHD 13 of a shock-forming solution launched by smooth initial data, with a precise description of the gradient-blowup along a singular boundary, a precise description of the classical solution's regular behavior along a Cauchy horizon C, and a sharp description of the structure of the "first singularity," which in Sect. 5.5 (also in Theorem 3.1), we referred to as "the crease" and which we denote by "∂ − B" in Fig. 6. Crucially, under the "transversal convexity" assumption described in Sect. 5.5, ∂ − B is a co-dimension 2, h-spacelike submanifold equal to the intersection of the Cauchy horizon and the singular boundary; without this structure, it is not even clear whether the shock development problem is well-posed. These structures, as well as the full state of the fluid up to ∂ − B, C, were posited as data in [31,33]. In 1D, Theorem 3.1 provides the first rigorous construction of these qualitative and quantitative aspects of the data for the shock development problem, starting from smooth initial conditions on a spacelike Cauchy hypersurface. For the non-relativistic Euler equations in 3D and without symmetry assumptions, our works [1][2][3] collectively provide an analog of Theorem 3.1.

Rarefaction waves.
In his foundational papers [6,7], for a large class of multi-dimensional hyperbolic systems that includes scalar conservation laws and the compressible Euler equations as special cases, Alinhac proved local existence and uniqueness for rarefaction wave solutions. Multi-dimensional rarefaction waves are analogs of a class of solutions to the well-known Riemann problem in 1D, in which the initial data are piecewise smooth and discontinuous, and the initial discontinuity is immediately smoothed out by the flow. His approach relied on Nash-Moser estimates to overcome derivative loss in linearized versions of the equations. Alinhac proved [6,7] that, as in the 1D case, in the corresponding multi-dimensional region + solution, the initial discontinuity is immediately smoothed out by the flow. In the important recent works [68,69], the authors study the isentropic 2D compressible Euler equations with a family of irrotational discontinuous initial data that are (asymmetric) perturbations of plane-symmetric data for a corresponding 1D Riemann problem. For irrotational data, their main result provides a sharpened version of the local existence and uniqueness results of Alinhac [6,7]. Their work shows in particular that the 2D irrotational rarefaction solution is a perturbation of the standard 1D rarefaction wave solution to the Riemann problem. Compared to [6,7], the works [68,69] yield two important advances. First, the techniques of [68,69] avoid a loss of derivatives in a corresponding linearized problem and consequently, the authors were able to close the energy estimates without Nash-Moser estimates. Second, [68,69] provide a complete description of the characteristic geometry in the problem. 6. Some ideas behind the proof of shock formation in 3D without symmetry In this section, we provide an outline (without complete proofs) of how to extend some aspects of Theorem 3.1 to apply to open sets of 3D relativistic Euler solutions without symmetry, isentropicity, or irrotationality assumptions. We anticipate that all aspects of Theorem 3.1 can be extended. A key reason behind our expectation is the availability of the new formulation of the flow provided by Theorem 4.1, which is qualitatively similar to the formulation of the non-relativistic flow from [94], which was used in prior related works [1-3, 65, 67] on shocks for the (non-relativistic) 3D compressible Euler equations. To keep the discussion short, we only sketch some key ideas behind the proof of Conjecture 1, which essentially concerns studying 3D solutions until their time of first blowup. Our discussion here mirrors the strategy we used in our work [67], in which we proved an analog of Conjecture 1 for the non-relativistic 3D compressible Euler equations.
6.1. Conjectures for the 3D relativistic Euler equations. To set the stage, we state three conjectures for the 3D relativistic Euler equations that are tied to the following figure: For the 3D compressible Euler equations, Fig. 6 has been justified [1][2][3] for open sets of initial data.

Conjecture 1 (Time of first blowup for the 3D relativistic Euler equations with vorticity and entropy).
Consider the simple isentropic plane-symmetric solutions of Theorem 3.1 as "background solutions" to the 3D relativistic Euler equations with symmetry. Consider smooth initial data in 3D -without symmetry, irrotationality, or isentropicity assumptions -that are close (in a sufficiently high-order Sobolev space) to the data of one of the background solutions. Then the perturbed solution forms a shock at a time that is a perturbation of T Shock , i.e., a perturbation of the shock-formation-time of the background solution; the conjectured perturbed first blowup-time is the smallest value of t along the crease ∂ − B depicted in Fig. 6.
Conjecture 2 (The structure of the crease and singular boundary for the 3D relativistic Euler equations with vorticity and entropy). Under the assumptions of Conjecture 1, the perturbed solution has a crease and singular boundary that are perturbations of the crease singular boundary of the background solution; the background crease and singular boundary are depicted in Fig. 3, while the conjectured perturbed singular boundary and crease are depicted in Fig. 6.
Conjecture 3 (The structure of the Cauchy horizon with vorticity and entropy). Under the assumptions of Conjecture 1, the perturbed solution has a Cauchy horizon that is a perturbation of the Cauchy horizon of the background solution; the background Cauchy horizon is depicted in Fig. 3, while the conjectured perturbed Cauchy horizon is depicted in Fig. 6. Conjecture 2 is strictly harder than Conjecture 1, while Conjecture 3 is independent (though closely related).

Almost Riemann invariants.
To study general (asymmetric) perturbations of simple isentropic plane-symmetric solutions, it is convenient to use analogs of the Riemann invariants from Def. 3.2.
where H > 0 is the constant (3), and in (152), we are viewing the speed of sound c as a function of H and s. We define the almost Riemann invariants R (+) and R (−) as follows: Note that for plane-symmetric solutions with s ≡ 0, the quantities defined by (153a)-(153b) coincide with the ones defined in (48a)-(48b).
We view R (+) and R (−) as replacements for H and u 1 that are convenient for studying perturbations of simple isentropic plane-symmetric solutions. In particular, for values of (H, s) near (H, 0), the factor 1 Hc(H,s) in (152) is positive. The implicit function theorem then allows one to solve for H as a smooth function of F and s. We can then use the relations (153a)-(153b) to express (H, u 1 ) as smooth functions of (R (+) , R (−) , s). Then, under the algebraic relation (12d), a complete set of state-space variables for the 3D relativistic Euler equations is given by: Although we have previously used the symbol ⃗ Ψ to denote the array h, u 0 , u 1 , u 2 , u 3 , s , in the rest of Sect. 6, we use ⃗ Ψ to denote the array in (154). Away from isentropic plane-symmetry, R (+) and R (−) no longer solve the homogeneous transport equations (49a)-(49b). However, from (147a), the chain rule, and Def. 6.1, it follows that for general solutions, R (+) and R (−) satisfy geometric wave equations that have the same schematic form as (147a). In the rest of Sect. 6, we use this fact without always explicitly mentioning it. 6.3. Nonlinear geometric optics. As in the 1D case treated in Sect. 3, the study of shock formation in 3D relies on nonlinear geometric optics, i.e., multi-dimensional analogs of the constructions from Sect. 3.6. The main object behind the construction is an eikonal function. The use of eikonal functions to study the global properties of multi-dimensional quasilinear hyperbolic PDE solutions originated Christodoulou-Klainerman's proof [32] of the stability of Minkowski spacetime. See also Sect. 5.3 for a discussion of nonlinear geometric optics in the context of proofs of multi-dimensional shock formation.
As we discussed at the beginning of Sect. 3.6, for simple isentropic plane-symmetric solutions, the role of nonlinear geometric optics/geometric coordinates is to yield a framework/differential structure in which the solution remains smooth all the way up to the shock, thereby allowing one to reduce the problem of shock formation to a more traditional problem in which one studies long-time existence. The same remarks hold in the present context of three spatial dimensions, although many additional geometric constructions are required.
In particular, all of the geometric quantities we introduce in Sects. 6.3 and 6.8, including the acoustic eikonal function, the geometric coordinates, various geometric vectorfields, and various geometric energies are introduced in order regularize the problem. However, unlike in the simple isentropic plane-symmetric case, in the present multi-dimensional context, one needs to derive energy estimates, and it turns out that a remnant of the shock singularity survives, even in the "good" geometric coordinates. The remnant manifests as singular estimates for the high-order geometric energies, though it is crucial for our approach that the mid-order-and-below geometric energies remain bounded; see Sect. 6.9.3 for further discussion.
Definition 6.2 (Acoustic eikonal function). An acoustic eikonal function (eikonal function for short) is a solution to the eikonal equation, which is the following fully nonlinear hyperbolic PDE, where h is the acoustical metric from Def. 2.4: When studying shocks close to plane-symmetry, it is convenient to assume the following initial condition for U , which is the exact same initial condition (55b) that we used in our study of isentropic plane-symmetric solutions: We refer to the level sets of U as the characteristics. One can show that for isentropic plane-symmetric solutions, the solution "U " to (55a)-(55b) coincides with the solution to (155)-(156). Definition 6.3 (Inverse foliation density). We define the inverse foliation density of the characteristics to be the following scalar function: where the second equality in (157) follows from (143b).
One can show that for isentropic plane-symmetric solutions, the quantity "µ" defined by (57) coincides with the one defined in (157). As in Sect. 3, the formation of a shock occurs when the density of the characteristics becomes infinite, that is, when µ ↓ 0. 6.4. Geometric coordinates and vectorfield frames tied the eikonal function. In this section, we introduce geometric coordinates and related vectorfields that are analogs of quantities we used in our study isentropic plane-symmetric solutions.
6.4.1. Geometric coordinates. Definition 6.4 (Geometric coordinates and the corresponding partial derivative vectorfields). We refer to {t, U, x 2 , x 3 } as the geometric coordinates on spacetime, where U is the eikonal function and t, x 2 , and x 3 are the Minkowski-rectangular coordinates, with t := x 0 . We denote the corresponding partial derivatives by { ∂ ∂t , ∂ ∂U , ∂ ∂x 2 , ∂ ∂x 3 } (which are not to be confused with the partial derivatives ∂ α in the (t, x 1 , x 2 , x 3 )coordinate system). Definition 6.5 (Some important submanifolds of geometric coordinate space).
• We denote the level set {U = U ′ } by P U ′ , and we refer to these level sets as the characteristics; see Fig. 7. • We define Σ t ′ to be the hypersurface {t = t ′ }, where t is the Minkowski time function.
The surfaces P U are characteristic for the wave operator □ h on LHS (147a), or equivalently, null with respect to the acoustical metric h (h-null for short).
As we mentioned in Sect. 3 and just above Def. 6.2, a key strategy behind the proof of shock formation is to show that the solution remains smooth with respect to the geometric coordinates. The blowup of the solution's gradient in (t, x 1 , x 2 , x 3 )-coordinates, i.e., the blowup of ∂ ∂ ∂ ⃗ Ψ, occurs when µ vanishes. One can prove (see e.g. [67, Lemmas 2.22 and 2.23]) the following relationship between the partial derivatives in the two coordinate systems, where the a β α are smooth (solution-dependent) functions of the geometric coordinates: The first product on RHS (158) shows in particular why ∂ ∂ ∂ ⃗ Ψ can blow up as µ ↓ 0, even when the derivatives of ⃗ Ψ with respect to elements of { ∂ ∂t , ∂ ∂U , ∂ ∂x 2 , ∂ ∂x 3 } remain bounded. The use of geometric coordinates allows one to study the problem of shock formation by using ideas that have traditionally been used in long-time existence problems. However, different from the 1D case, in 3D, some difficult degeneracies can occur in the energy high-order energy estimates; we discuss this difficulty in more detail in Sect. 6.9.3. 6.4.2. Vectorfield frames. Experience has shown (e.g., [30] and the related works cited in Sect. 5.3) that in the study of shocks, to avoid loss of derivatives in commutator estimates and other difficulties, it is advantageous to use vectorfield frames adapted to the eikonal function. In particular, once one has constructed U , one can construct the following vectorfield frames, depicted in Fig. 7: which we describe below. The frame Z spans the tangent space of geometric coordinate space, the subset P spans the tangent space of the characteristics P U , and the subset Y spans the tangent space of the two-dimensional surfaces ℓ t,U , i.e., the span of Y is equal to that of Figure 7. The frame Z at two distinct points on P U with the x 3 -direction suppressed We now briefly describe how to construct the vectorfields. First, L is defined by L α := −µ(h −1 ) ακ ( ⃗ Ψ)∂ κ U . L is clearly h-orthogonal to the characteristics P U , and by (155) and (157), we see that h(L, L) = 0 and Lt = 1. In particular, since L is h-orthogonal to P U and h-null, P U is a h-null hypersurface. One can show that for isentropic plane-symmetric solutions, the quantity "L" defined in this fashion coincides with the one defined in (44a). However, for general solutions in 3D, a big difference with the plane-symmetric case is that there does not exist any explicit formula for L in the spirit of (44a). That is, in general, L depends on the gradient of U , which in turns depends implicitly on the fluid solution via the coefficients in the eikonal equation (155).
Next, we define X to be the Σ t -tangent vectorfield that is h-orthogonal to the ℓ t,U , normalized by h(X, X) = 1, and is such that U increases along the integral curves of X, i.e., XU > 0. We then de-fineX := µX, where µ is defined in (157). It follows that h(X,X) = µ 2 . With the help of (19), one can show thatXU = 1. That is,X is a geometric replacement 14 for ∂ ∂U , and among the elements of Z , it is the 14 (63) shows that in plane-symmetry,X = ∂ ∂U . However, in the general 3D case,X is equal to ∂ ∂U plus corrections belonging to span one vectorfield that is transversal to the characteristics. One can show that for isentropic plane-symmetric solutions, the quantities "X" and "X" defined in this fashion coincide with the ones defined in (60). Next, for A = 2, 3, we define Y (A) to be the h-orthogonal projection of ∂ A onto ℓ t,U , where ∂ A is the standard Minkowski-rectangular partial derivative vectorfield. In our study of isentropic plane-symmetric solutions, we did not rely on any analogs of the Y (A) . 6.4.3. The first fundamental form of ℓ t,U and related decompositions. Standard calculations (see, e.g., [95,Equation (2.40b)]) imply that the inverse acoustical metric from (17b) can be decomposed as follows: where g / is the Riemannian metric on ℓ t,U induced by h, i.e., the first fundamental form of ℓ t,U . We view g / to be a spacetime tensor that vanishes on contraction with L or X and that agrees with h on ℓ t,U -tangent vectors. Moreover, in (160), (g / −1 ) αβ := (h −1 ) αγ (h −1 ) βδ g / γδ is the corresponding inverse first fundamental form.
In what follows, ∇ / denotes the Levi-Civita connection of g / and | · | g / denotes the pointwise norms of tensors with respect to g /, e.g., for scalar functions f , we have |∇ / f | 2 g / = (g / −1 ) αβ ∇ / α f ∇ / β f . For future use, we note that for scalar functions f , we have the following decomposition of µ□ h( ⃗ Ψ) f (see, e.g., [95,Proposition 2.16]), where ∆ / denotes the covariant Laplacian associated to g /: On RHS (161), · · · denotes terms that depend on the first derivatives of f and that we will not discuss in detail here. (of the characteristics, with respect to Σ t ). III) X could be defined to be a Σ t -tangent vectorfield that is h-orthogonal to ℓ t,U and normalized by h(X, X) = 1 (this would define X up to an overall minus sign, where in Sect. 3, we defined X such that is is "left-pointing"). We do, however, highlight the following point: general spacetime manifolds have no symmetries and thus do not generally allow for the existence of the simple plane-symmetric solutions constructed in Sect. 3. Nonetheless, we expect the robust geometric framework described above and many of the key ideas presented in the subsequent sections to be useful for proving shock formation results for initial data with large gradients. We also note that the initial data (156) for the eikonal function generally needs to be adjusted to fit the regime one is studying. 6.5. Description of the compactly supported initial data of perturbations of simple isentropic plane-symmetric solutions. We now describe the initial data in Conjecture 1 in more detail. One simply takes any of the "background" simple isentropic plane-symmetric shock-forming solutions from Theorem 3.1 and then considers a general small (asymmetric) perturbation of its initial data on Σ 0 . For the proof to close, the perturbed initial data need to belong to a sufficiently high order Sobolev space; we will discuss the issue of regularity in more detail later on. For convenience, it is easiest to consider perturbed initial data that, like the background data, are compactly supported in Σ Then by finite speed of propagation, the corresponding perturbed solution remains spatially compactly supported for all time and in particular vanishes along the characteristic P −U1 (which therefore appears flat in Fig. 6). The compact support allows one to avoid difficult boundary terms in the elliptic estimates for the vorticity and entropy; see Sect. 6.9.2.
For such initial data, all of the wave variables in the array ⃗ Ψ defined in (154) -except for R (+) -will initially have a negligible effect on the dynamics. Similarly, the fluid variables ϖ α , S α , C α , and D from Sect. 4.2 will also be small initially. The challenge is to propagate suitable versions of this smallness all the way up to the shock. 6.6. Bootstrap assumptions and the regions M [0,t ′ ),[−U1,U ′ ] . To prove that a shock forms, one commutes the equations up to N top times and derives energy estimates up to top-order and L ∞ estimates at the lower orders. For reasons described later on, the proof requires N top to be rather large (i.e., substantially larger than a proof of local well-posedness requires) and thus the initial data need to belong to a high order Sobolev space. To close the proof, it is convenient to make L ∞ bootstrap assumptions that capture the expectation that at the lower derivative levels, the solution behaves like a perturbation of a simple isentropic plane-symmetric solution (in which only R (+) is non-vanishing). That is, one aims to propagate various aspects of the smallness enjoyed by the initial data described in Sect. 6.5. The bootstrap assumptions are convenient because explicit formulae in the spirit of Cor. 3.6 are not available in 3D. They capture the expectation that nonetheless, the 3D solution should obey estimates that are similar to the ones proved in Lemma 3.7 and Theorem 3.1 in the simple isentropic plane-symmetric case.
Before introducing the bootstrap assumptions, we first introduce some notation. In what follows, P N denotes an arbitrary order N string of elements of the P U -tangent vectorfields P defined in (159). We refer to such vectorfields as "tangential" (to the characteristics P U ). In contrast, the vectorfieldX from (159) The proof will show that µ first vanishes at a time that is a perturbation of the plane-symmetric blowup-time T Shock from Theorem 3.1. That is, in all estimates and arguments, one can safely assume that 0 < T Boot ≤ 2T Shock . From this assumption, the assumption that the data are supported in Σ , and finite speed of propagation, one can compute/choose U 3 > U 2 such that for t ∈ [0, T Boot ), the solution is trivial 15 at points in Σ t belonging to the complement of Σ t ∩ {U ∈ [−U 1 , U 3 ]}.
Specifically, for N up to mid-order, i.e., roughly, for N ≤ Ntop 2 , one assumes L ∞ smallness on M [0,TBoot),[−U1,U3] for the following quantities, which, aside from Lµ, vanish in the case of simple isentropic plane-symmetric solutions: • The pure P u -tangential derivatives P N R (+) , P N ϖ α , P N S α , P N C α , P N D, P N L i , and P N µ, except Lµ does not have to be small (cf. (119)). • Mixed tangential-transversal derivatives of all these quantities except µ, e.g., LXR (+) .
One also assumes L ∞ boundedness (not smallness) on M [0,TBoot),[−U1,U3] for the following quantities, which are non-zero for the simple isentropic plane-symmetric solutions we studied in Theorem 3.1: One also assumes that: The assumption (163) captures that no shock has occurred on M [0,TBoot),[−U1,U3] , though it leaves open the possibility that the shock happens exactly at time T Boot . The difficult part of the proof is proving energy estimates that allow one to derive, through a combination of Sobolev embedding and transport estimates, improvements of the bootstrap assumptions. Then through a standard continuation argument, one can extend the solution all the way up until the first time that µ vanishes. 6.7. The behavior of µ and the formation of a shock. Given the framework established above and the bootstrap assumptions, it is not difficult to prove that shocks can form for open sets of initial data that are perturbations of simple isentropic plane-symmetric data. We now sketch the argument. Using the eikonal equation (155) and the bootstrap assumptions of Sect. 6.6 one can show that for perturbations of simple isentropic plane-symmetric solutions, the inverse foliation density defined in (157) where · · · denotes small error terms. In addition, using that LXR (+) is small by the bootstrap assumptions, one can show that on M [0,TBoot),[−U1,U3] , we have: that is,XR (+) stays close to its initial condition (cf. (69), which shows that in simple isentropic planesymmetry R (+) depends only on U ). Combining (164)-(165), we see that on M [0,TBoot),[−U1,U3] , we have: Since Lt = 1 (i.e., L = d dt along the integral curves of L), and since the ℓ t,U -tangential derivative ∂ ∂x A µ is small by the bootstrap assumptions, we deduce from (166) and the fundamental theorem of calculus that on M [0,TBoot),[−U1,U3] , we have: Thus, for initial data such that there are points whereXR (+) (0, U, x 2 , x 3 ) is negative and large enough to dominate the terms · · · on RHS (167), we infer from (167) that min Σt µ will vanish in finite time, and that the time of first vanishing can be controlled by the initial data. Moreover, since this argument implies that |XR (+) | ̸ = 0 at the points where µ vanishes, and since: it follows that |XR (+) | blows up like C µ at the points where µ vanishes (cf. 128 and (158)). In total, these arguments show why µ can vanish in finite time and reveal how the vanishing is connected to the blowup of |XR (+) |.
We also note that the above arguments yield the following: which is an analog of the estimate (119) in simple isentropic plane-symmetry. In particular, at fixed (U, x 2 , x 3 ), µ vanishes linearly in t. These facts have important implications for the energy estimates, described below; see in particular Remark 6.3.
6.8. Vectorfield multiplier method and energies. In this section, we describe some basic ingredients that are needed to set up the energy estimates. To control the solution up to the shock, one needs geometric machinery that is much more advanced than that of Sect. 2.3.4 and that takes advantage of the structure of the equations provided by Theorem 4.1.
6.8.1. Energies and null fluxes for wave equations. In this section, we set up the vectorfield multiplier method for the wave equations (147a).
Given a scalar function f , we define the energy-momentum tensor associated to it to be the following symmetric type 0 2 tensor: In (170) and throughout, D denotes the Levi-Civita connection of h.
Straightforward calculations, based on the symmetry property D α D β f = D β D α f and the Leibniz rule, yield that given any multiplier vectorfield Z, we have the following identity: where: is the deformation tensor of Z (with respect to h). The identity (171) is convenient for bookkeeping in the divergence theorem when one integrates by parts. To derive energy-null flux identities for wave equations, we will use (171) with the multiplier vectorfield Z equal toT , which is defined as follows, where L andX are as in Sect. 6.4.2: One can compute that h(T,T ) = −4µ(1 + µ) and thusT is h-timelike whenever µ > 0. The h-timelike nature ofT is crucial for generating coercive energies. The µ-weights in (173) have been carefully placed.
To close the L 2 estimates, one needs energies on: and null fluxes on: Recall that N is the future-directed h-unit normal to Σ t defined in (143b). The strength of our energies on the Σ is determined by the following identity (see, e.g., [95,Lemma 3.4]): Similarly, the strength of our null fluxes on the P [0,t) U is determined by the following identity (again, see, e.g., [95,Lemma 3.4]): These identities (176)-(177) motivate the following definition.
where g / and |∇ / f | 2 g / are as in Sect. 6.4.3, and on each ℓ t ′ ,U ′ , is the canonical volume form induced by g /. In particular, the determinant in (179) is taken relative to the coordinates (x 2 , x 3 ) on ℓ t ′ ,U ′ .
Remark 6.2 (µ weights and error terms). Note that one L-derivative involving term on RHS (178b) lacks a µ-weight. This means that F (Wave) [f ](t, U ) is strongly coercive in the L-derivatives even in regions where µ is small. This is crucial for controlling L-derivative-involving energy estimate error terms that lack µ-weights.
In contrast, all terms on RHS (178a)-(178b) involving the ∇ / -derivatives of f have a µ-weight. Hence, to control dangerous ∇ / -derivative-involving energy estimate error terms that lack µ-weights, one must rely on the special structure afforded by the coercive spacetime integral described in Remark 6.3.
Using (171) with Z :=T , one can prove the following lemma (see, e.g., [95,Proposition 3.5]), which forms the starting point for the L 2 -analysis of solutions to the wave equations (147a).
Remark 6.3 (The coercive spacetime integral K[f ](t, U )). Upon decomposing the Q-involving integral on RHS (180) relative to the vectorfields in Z , one finds that it contains the following term, which, after moving it to LHS (180), takes the following form: In (181), [z] − := max{0, −z} denotes the negative part of z. Importantly, the discussion in Sect. 6.7 shows that Lµ is negative in regions close to the shock (i.e., in regions where µ is small) and thus in such regions, K[f ](t, U ), which we again emphasize has the favorable sign (181) when moved to LHS (180), yields non-µweighted spacetime L 2 control of |∇ / f | 2 g / . As we mentioned above, the control afforded by this term is crucial for controlling error terms in the energy estimates that depend on ∇ / f but lack µ-weights.
Note that the discussion below (146) implies that the factor −nm(B, L) on RHS (182b) is positive, which of course is important for the coercivity of F (Transport) [f ](t, U ). For the same reasons given in Sect. 6.8.1, it is crucially important that RHS (182b) lacks a µ-weight.
By applying the divergence theorem to the vectorfield f 2 B on M [0,t),[−U1,U ] , using (145)-(146), and using that in geometric coordinates we have (see, e.g., [91,Corollary 3.47]): one can prove the following lemma, which forms the starting point for the L 2 -analysis of solutions to the transport equations (147b)-(147c).
6.8.3. Commutator method and the shock-driving terms. To derive suitable energy estimates, one commutes µ-weighted 17 versions of the equations of Theorem 4.1 with the elements of the P U -tangent subset P defined in (159) and then applies the energy-null flux identities of Lemmas 6.1 and 6.2 and the elliptic estimates described in Sect. 6.9.2. In particular, it has been understood since [67] that to close the energy estimates in a shock formation problem, one does not need to commute the equations with the vectorfieldX, which is transversal to the characteristics. We stress that it does not seem possible to close the estimates by commuting the equations with the Minkowski-rectangular partial derivatives ∂ α ; such an approach would generate error terms that we do not know how to control. Although one must commute with all the elements of P to close the energy estimates, the most difficult terms arise from commuting the wave equations (147a) with a string of elements of the ℓ t,U subset Y defined in (159). Schematically, for Ψ ∈ ⃗ Ψ, with Y N denoting an order N string of elements of Y , we have: where χ is the null second fundamental form of P U , tr g / χ denotes its trace with respect to g /, and · · · denotes easier error terms.
The Y N C-and Y N D-involving terms on RHS (185) clearly arise from the first term on RHS (147a), and we will discuss how to control them in Sect. 6.9.2. We refer to the proofs of [95, Lemmas 2.18 and 4.2] for details on the origin of the commutator term (XΨ) · Y N tr g / χ, which we describe how to control just below; this term would lead to the loss of a derivative at the top order if handled in a naive fashion.
The key point is that there are not any terms on RHS (186) that are proportional toXψ ·Xϕ; signature considerations imply that such terms, if present, would be multiplied by a dangerous factor of 1 µ . Such a factor would become singular as µ ↓ 0, which would spoil the philosophy of proving that the solution remains quite regular in geometric coordinates, even as µ vanishes. A decomposition similar to (186) also holds for the up-to-top-order derivatives of the µ-weighted h-null forms. We stress that the absence of terms proportional toXψ ·Xϕ on RHS (186) is a statement about the full nonlinear structure of the h-null forms from Def. 4.6. In particular, this notion of null form is tied to the acoustical metric h and is stronger than the "classic null condition" established by Klainerman in [54], which is tied to the Minkowski metric and is indifferent to the structure of most cubic nonlinearities.
Remark 6.4 (The shock-driving terms are hidden in □ h( ⃗ Ψ) Ψ). The upshot is that all terms in (147a) that drive the formation of the shock are hidden in the covariant wave operator term □ h( ⃗ Ψ) Ψ on the LHS. More precisely, the shock-driving terms are semilinear Riccati-type terms of type XΨ · XΨ (cf. Sect. 3.5) that become visible if one first expands LHS (147a) relative to the Minkowski-rectangular coordinates and then decomposes the semilinear terms with respect to the frame Z from (159).
In view of the above discussion, it follows that the goal is to show that all terms on RHS (185) are error terms that do not interfere with the shock formation processes. As we have alluded to earlier, this is indeed possible except at the high derivative levels, where the difficult regularity theory of the eikonal function (in particular, the difficult regularity theory of χ) leads to singular high-order estimates; see Sect. 6.9.3.
We now explain how to control the first product on RHS (185). One starts with the Raychaudhuri equation [81], which plays an important role in mathematical General Relativity, and which takes the following form in the present context: 17 The µ-weights allow one to avoid problematic error terms, and they are compatible with the µ-weights in Lemmas 6.1 and 6.2.
In (187), Ric is the Ricci curvature of the acoustical metric h and Ric LL := Ric αβ L α L β . The difficulty is that since h = h( ⃗ Ψ), Ric depends on the second derivatives of ⃗ Ψ. This suggests that control of the term Y N tr g / χ on RHS (185) requires control over N + 2 derivatives of ⃗ Ψ, which is one too many derivatives to be compatible with the regularity for ⃗ Ψ afforded by energy estimates for the wave equation (185), i.e, one has to worry about losing a derivative at the top-order.
Fortunately, based on ideas originating in [32,55], there is a way to avoid the loss of derivatives by combining the special structure of Ric LL with the special structures of the equations of Theorem 4.1. We will briefly explain the main ideas of the argument. First, one decomposes Ric LL to deduce an algebraic identity that can schematically be expressed as: where f schematically denotes a smooth function that is allowed to vary from line to line and · · · denotes lower order terms that depend on the first derivatives of ⃗ Ψ. Next, one uses (161) and (164)  ⃗ Ψ with the inhomogeneous terms on RHS (147a). Combining with (187) and bringing all perfect L-derivative terms to the left, we deduce a transport equation that can schematically be expressed as: where · · · denotes easier error terms that involve ϖ, S, and the first derivatives of ⃗ Ψ. In particular, RHS (189) does not involve the second derivatives of ⃗ Ψ, which is its main advantage compared to equation (187). Similar results hold for the commuted equations, i.e., given any order N string of ℓ t,U -tangent vectorfields Y N , with (Y N ) X denoting the "fully modified" order N version of tr g / χ defined by: one can derive a transport equation of the following schematic form (see, e.g., the proof of [91, Proposition 11.10] for the main ideas behind the proof): Equation (191) is an order N version of equation (189) that allows one to control (Y N ) X without derivative loss, assuming that one can adequately control RHS (191). In Sect. 6.9.2, we will explain how to control the Y N C-and Y N D-involving terms on RHS (191). One remaining difficulty is that the term µY N (|χ| 2 g / ) on RHS (191) involves the full tensor χ, rather than just the trace part featured on LHS (191). However, one can control this term using a strategy that goes back to [32,55]. Specifically, one can use elliptic estimates on the two-dimensional surfaces ℓ t,U , based on the Codazzi equations from geometry, to control the ℓ t,U -tangent derivatives of |χ| 2 g / in terms of the ℓ t,U -tangent derivatives of tr g / χ plus error terms with an admissible amount of regularity; similar results hold up to top-order, and we refer, for example, to the proof [91,Lemma 20.20] for the main ideas behind the analysis. In total, the approach described above allows one to control the product (XΨ) · Y N tr g / χ on RHS (185) in appropriate L 2 -based Sobolev spaces. However, the argument introduces a difficult factor of 1 µ into the top-order estimates, as we explain in Sect. 6.9.
6.9. Energy estimates. In this section, we will describe some features of the proof of the energy estimates, which are the most technical and difficult aspect of studying shock formation.
6.9.1. Transport variable energy estimates. We now explain how to derive energy estimates for (ϖ, S), which solve the transport equation (147b). These estimates are relatively easy to derive. Let P N denote an arbitrary order N string of elements of the P U -tangent vectorfields P defined in (159). One multiplies equation (147b) by µ, commutes with P N for N ≤ N top , and uses the bootstrap assumptions and the transport equation energy identity (184) to deduce that for (t, U ) ∈ [0, T Boot ) × [−U 1 , U 3 ], we have: From (192) and Grönwall's inequality in U , we deduce that: where here and throughout, "data" denotes a small term that depends on the initial data and that measures the perturbation of the data away from simple isentropic plane-symmetric data. The terms "µ-regular ⃗ Ψ terms" in (193) are easily controllable by the wave energies and null fluxes, and in particular, no singular factor of 1 µ is present in these terms. Hence, the estimate (193) shows that the behavior of P ≤Ntop (ϖ, S) is effectively determined by the behavior of the wave variables ⃗ Ψ (defined in (154)). We clarify that even though there is no singular factor of 1 µ present in the µ-regular ⃗ Ψ terms, these terms can still blow up at the high derivative levels as µ ↓ 0, for reasons we describe in Sect. 6.9.3. Hence, in view of the coupling between the transport and wave energies shown by (193), we see that the blowup of the ⃗ Ψ-energies can cause the blowup of the (ϖ, S)-energies.
6.9.2. Elliptic-hyperbolic estimates for the vorticity and entropy. Recall that the terms Y N C and Y N D appear on RHSs (185) and (191). To close the energy estimates, we need to control P ≤Ntop (C, D) in L 2 . These terms would lead to a loss of derivatives at the top-order if handled in a naive fashion. There are several reasons why one cannot control these terms by using only the same transport equation energy methods we used to derive (193). One reason is that Def. 4.4 shows that from the point of view of regularity, we have (C, D) ∼ ∂ ∂ ∂(ϖ, S) + · · · ; since the inhomogeneous terms in the transport equation (147b) satisfied by (ϖ, S) suggest (incorrectly, as it fortunately turns out) that ∂ ∂ ∂(ϖ, S) have at best the same regularity as ∂ ∂ ∂ 2 ⃗ Ψ, this is formally inconsistent (from the point of view of regularity) with having ∂ ∂ ∂(ϖ, S) as a source term in the wave equations (147a) for ⃗ Ψ. A second reason is that the transport equations (147c) satisfied by (C, D) feature the inhomogeneous terms Q(∂ ∂ ∂ϖ, ∂ ∂ ∂ ⃗ Ψ) and Q(∂ ∂ ∂S, ∂ ∂ ∂ ⃗ Ψ), which depend on the general first derivatives of (ϖ, S), rather than the special combinations of first derivatives of (ϖ, S) present in the definitions (140a)-(140b) of (C, D). That is, due to these inhomogeneous terms, (147c) cannot be treated as a pure transport equation in (C, D).
To overcome the difficulties highlighted in the previous paragraph, one can treat (147b)-(147d) as a coupled transport-div-curl system that yields L 2 -control of (C, D) and (∂ ∂ ∂ϖ, ∂ ∂ ∂S), where (C, D) are controlled with hyperbolic transport energy estimates in the spirit of (193), and the elliptic estimates are used only to handle the inhomogeneous terms Q(∂ ∂ ∂ϖ, ∂ ∂ ∂ ⃗ Ψ) and Q(∂ ∂ ∂S, ∂ ∂ ∂ ⃗ Ψ) on RHS (147c). Similar results hold up to top-order and yield L 2 control over ∂ ∂ ∂P ≤Ntop (ϖ, S) and P ≤Ntop (C, D). One key difficulty in this argument is that to obtain elliptic estimates on the spacelike hypersurfaces Σ t , one needs to extract a spatial div-curl subsystem from (147b)-(147d); the difficulty is that equations (147b)-(147d) appear to involve spacetime divcurl equations. Nevertheless, by splitting various derivative operators into a B-parallel part and a Σ t -tangent part, one can extract the desired spatial div-curl subsystem; we refer to the proof of [41, Lemma 9.21] for details on this extraction in the context of a local well-posedness argument. We also refer to [67, Sections 11.2 and 11.3] for similar but simpler elliptic estimates in the context of shock formation for the non-relativistic 3D compressible Euler equations.
The argument described in the previous paragraph applies when the initial data of the vorticity and entropy are compactly supported, which is the case for the initial data described in Sect. 6.5. The key point is that under compact support, one can avoid boundary terms in the elliptic estimates (the boundary terms vanish thanks to the compact support). To treat solutions that are not (spatially) compactly supported, one would need to handle the boundary terms that arise in the elliptic estimates. In the case of nonrelativistic 3D compressible Euler equations, based on the special structures found in the formulation of the flow derived in [94], localized spacetime "elliptic-hyperbolic" integral identities for the vorticity and entropy were derived in [3]. The identities of [3] allow one, in the non-relativistic case, to handle the boundary terms. In particular, in [1,2], we used specialized versions of those identities to study the structure of h-MGHD for the non-relativistic 3D compressible Euler equations. To extend these results to 3D relativistic Euler solutions without compact support, one would need to derive relativistic analogs of the integral identities from [3]. 6.9.3. The wave energy estimate hierarchy. We now discuss the energy estimates for the wave variables ⃗ Ψ, defined in (154). As before, let P N denote an arbitrary order N string of elements of the P U -tangent vectorfields P defined in (159). The main difficulty in closing the energy estimates in multi-dimensions is that the best estimates we know how to derive allow for the possibility that the top-order energies blow up as µ ↓ 0. Before proceeding, for notational convenience, we define: where E (Wave) and F (Wave) are as in Def. 6.6 and K is as in Remark 6.3. The estimates for the wave energies and null fluxes take the following hierarchical form on (t, U ) ∈ [0, T Boot )×[−U 1 , U 3 ] (see, e.g., [65,Proposition 14.1] for complete proofs in the case of the 2D non-relativistic compressible Euler equations with vorticity): . . .
where as before, "data" denotes a small term that depends on the initial data and that measures the perturbation of the data away from simple isentropic plane-symmetric data, µ ⋆ (t, U ) := min{1, min denotes an arbitrary string of elements of the P U -tangent vectorfields P of order in between 18 1 and A 2 , and for reasons described below, A ≫ 1 is a universal constant (independent of N top , the initial data, and the equation of state). Energy estimates in the spirit of (195a)-(195d), which are singular at the high derivative levels, are the only kinds of energy estimates that are known in the context of multi-dimensional shock formation. One might be concerned that the high-order energies are allowed to blow up when µ vanishes, as this seems to be in conflict with the philosophy described in Sect. 6.3, namely that the solution should look regular in geometric coordinates, all the way up to the shock. However, the full hierarchy (195a)-(195d) shows that the wave energies become less singular by two powers of µ −1 ⋆ with each level of descent below the top, 18 One can close the proof without deriving energy estimates in the case N = 0. This is convenient because the order 0 wave energies can initially be large, stemming from the initial largeness of the data forXR (+) (largeness that is present even for the simple isentropic plane-symmetric solutions treated in Theorem 3.1). In contrast, the energies in (195a)-(195d) are initially small (in fact, they vanish for the solutions in Theorem 3.1).
until one reaches the level (195d) at which the energies remain bounded. In particular, the boundedness of the mid-order and below geometric energies, as shown by the estimate (195d), capture the sense in which the solution remains regular relative to the geometric coordinates. Due to coupling (see, e.g., (193)), these estimates imply that the fluid variables C, D, ϖ, S obey a similar energy hierarchy (omitted here for brevity, but see [67] for a detailed statement and proof of the hierarchy in the case of the non-relativistic 3D compressible Euler equations) featuring related but distinct blowup-rates, which are nonetheless compatible with proving (195a)-(195d). We highlight the following key point: The non-singular estimates (195d) are what allow one to improve, through Sobolev embedding, a smallness assumption on "data", and derivative-losing transport-equation-type estimates, the L ∞ bootstrap assumptions described in Sect. 6.6. The singular top-order wave equation energy estimate (195a) stems from the following integral inequality, whose proof we describe in Sect. 6.9.4: In (197), A ≫ 1 is the universal positive constant mentioned above and · · · denotes similar or easier 19 error terms. From (197) and Grönwall's inequality, we deduce that: An argument based on precise, refined versions of (166)-(167) yields the following crucial bound, which can be proved using the same arguments given in [95, Section 10] (the proof is trivial if one ignores the terms · · · in (166)-(167)): where µ ⋆ is defined in (196). Combining (198) and (199), we conclude (195a). The less singular below-top-order estimate (195b) stems from the following integral inequality, whose proof we describe in Sect. 6.9.5: Note that RHS (200) involves a coupling between the top-order energies and the just-below-top-order energies. The actual proof of (195b) involves a delicate Grönwall argument that is coupled to the proof of the top-order estimate (195a) as well as the following estimates, which hold for numbers B > 1 and which can be proved by using arguments similar to the ones needed to prove (199): From (178a), it follows that the first integral on RHS (206) is bounded by the A-multiplied integral on RHS (197) (this integral contributes a "portion of" the "A"). The (Y N top ) X -involving integral on RHS (206) can be bounded by a similar error term involving a double time-integral, where one of the time integrations comes from integrating (191); we refer to the proof of [95,Proposition 14.2] for the details. The part of the first integral on RHS (205) that is generated by the (1 + 2µ)L term in (173) can be handled through an argument involving integration by parts in L, which leads to difficult critical-strength boundary terms that make a similar contribution to the blowup of the top-order energies as the A-multiplied integral on RHS (197); the same arguments given in the proof of [95,Proposition 14.2] can be used to handle these terms. The Y Ntop C-and Y Ntop D-involving integrals on RHS (205) are easy to handle, thanks to the helpful factor of µ in the integrand and the strategy for bounding Y Ntop C and Y Ntop D via transport-div-curl estimates that we described in Sect. 6.9.2. We now consider the remaining cases, namely Ψ ∈ {R (−) , s, u 2 , u 3 }. The identity (205) holds in these cases too. However, the factor ofXΨ in the identity is now very small in L ∞ in these cases, since we are considering perturbations of simple isentropic plane symmetric solutions (for which these quantities all vanish). Hence, due to the smallness, all error terms on RHS (205) are much easier to treat and can in fact be relegated to the error terms "· · · " on RHS (197); see [65, page 154] for further details. 6.9.5. Discussion of the proof of the just-below-top-order inequality (200). We now sketch some key ideas behind the proof of the just-below-top-order inequality (200). First, one uses the identity (205), which, for any Ψ ∈ ⃗ Ψ, holds with P Ntop−1 in the role of Y Ntop . The key step, which is the one that is different compared to the proof of (197), is that one bounds the error integral on LHS (206), namely (T P Ntop−1 Ψ) · (XΨ) · P Ntop−1 tr g / χ dϖ g / dU ′ dt ′ , in a different way. First, one uses the bootstrap assumptions to bound the factor |XΨ| in (207) by ≲ 1. Next, one uses (173), the Cauchy-Schwarz inequality, and (178a) to bound the integral in (207) by: where we have used that the P U -tangent derivative terms in the energy (178a) contain a µ-weight, while there is an L factor in the definition ofT that does not contain a µ weight. Next, to bound the factor ∥Y Ntop−1 tr g / χ∥ L 2 (Σ [−U 1 ,U ] t ′ ) , one commutes the Raychaudhuri equation (187) with P Ntop−1 and uses the schematic relation Ric LL = P 2 ⃗ Ψ + · · · to deduce the following evolution equation, schematically depicted: LY Ntop−1 tr g / χ = P Ntop+1 ⃗ Ψ+· · · , where · · · denotes lower order error terms. This equation represents a loss of one derivative in the estimates and leads to the coupling between the different order energies described below (200). However, such a loss is permissible below the top-order. Since Lt = 1 and LU = 0, we can integrate this evolution equation and use (178a) to deduce that: where to obtain the last line of (209), we have again used that the P U -tangent derivative terms in the energy (178a) contain a µ-weight. Inserting (209) into (208), we conclude (200).

Open problems
In this section, we describe various open problems tied to shocks.
(1) Prove Conjecture 1. In Sect. 6, we outlined how to achieve this. We are confident that the conjecture can be proved, especially since the formulation of relativistic Euler flow provided by Theorem 4.1 is qualitatively similar to the formulation of 3D non-relativistic compressible Euler flow derived in [94], a system for which shock formation results have been derived (see Sect. 5).
Nonetheless, there are many non-trivial details that have to be checked, and this would be a good project for someone who wants to learn the field. (2) Prove Conjecture 2. This problem is much more technically demanding than proving Conjecture 1.
However, for the same reasons mentioned above, we expect that this can be achieved by adapting the methods used in the non-relativistic work [1] to the equations of Theorem 4.1. (3) Prove Conjecture 3. This problem is also much more technically demanding than proving Conjecture 1. For the same reasons mentioned above, we expect that this can be achieved by adapting the methods that we are using in our forthcoming non-relativistic work [2] to the equations of Theorem 4.1. (4) The shock development problem. Recall that in [31], Christodoulou solved the restricted shock development problem for the relativistic Euler equations and the non-relativistic compressible Euler equations in an arbitrary number of spatial dimensions; see Sect. 5.6. The general problem, i.e., the shock development problem with vorticity and entropy, remains an outstanding open problem for both systems in two or more spatial dimensions. (5) The global behavior of solutions with shocks. After the shock development problem (which is local-in-time) is solved, a particularly compelling problem will be to understand the global -in-timeand-space behavior of the corresponding weak solutions, at least in a perturbative regime. Among the many outstanding challenges in this problem is that of understanding the long-time behavior of the vorticity and entropy; even in the context of smooth incompressible non-relativistic flows, the long-time behavior of the vorticity is not well understood. We do, however, mention that in the breakthrough work [26], for axisymmetric incompressible non-relativistic flows in 3D with a boundary, smooth initial conditions were identified such that the solution (including the vorticity) blows up in finite time in an approximately self-similar fashion; such vorticity blowup, if present in multi-dimensional relativistic Euler flow, would be a serious obstacle to finding a meaningful way to continue the solution weakly past the singularity. (6) Extending the results to the coupled Einstein-Euler system. It is of great interest to extend the above results to the Einstein-Euler system in multiple spatial dimensions. The expectation is that while the fluid will exhibit shock wave phenomena (as in the uncoupled problem), the gravitational metric will exhibit less singular behavior As of present, the only known singularity formation result for the coupled system is in the plane symmetric case [82], in which the dynamics are described by 1 + 1-dimensional hyperbolic PDEs. In [82], it was shown that for many equations of state and a large class of plane-symmetric initial data, the Einstein-Euler solution breaks down in finite time. Although [82] provided heuristic arguments suggesting that the singularity is of shock-type, a precise description of the singularity was not given. The multi-dimensional problem is a difficult PDE problem because surfaces that are null or "barely spacelike" with respect to the acoustical metric are in fact timelike with respect to the gravitational metric; the reason is that speed of sound is slower than the speed of propagation of gravitational waves (at least when the speed of sound is less than unity). A corresponding key difficulty that arises in the context of the PDE estimates is that one must control the spacetime metric on various gravitationally timelike surfaces on which the fluid is singular; this is difficult at the top derivative level because generally, energy estimates for the spacetime metric are not available on surfaces that are timelike with respect to the spacetime metric. In contrast, in multiple wave speed systems such that the shock forms in the fastest wave, it is possible to close the top-order energy estimates and thus prove stable shock formation for the coupled problem [92]. (7) Extending the results to more complicated multiple speed systems. It is also of great physical and mathematical interest to extend the above results to other multiple speed systems, such as the equations of compressible magnetohydrodynamics, the GRMHD equations (i.e., the coupling of Einstein's equations to the equations of relativistic magnetohydrodynamics), the equations of elasticity, the equations of nonlinear electromagnetism, and the equations of crystal optics. A key difficulty of these systems is that their principal symbols are more complicated than wave operators and transport operators. Relatedly, their corresponding geometry is more complicated than Lorentzian geometry, and the characteristics comprise multiple sheets, which can be singular even at the tangent space level. These difficulties are a serious obstacle to implementing the kind of sharp version of nonlinear geometric optics (i.e., eikonal functions) that has been so successfully employed to study shocks for wave equations and fluid equations. (8) Implosion singularities. In the breakthrough works [72,73], the authors proved implosion singularity formation for some initially C ∞ spherically symmetric solutions to the compressible Euler equations and Navier-Stokes equations under adiabatic equations of state p = ρ γ for γ > 1, outside of a countable set of γ-values. In [19], the results were extended to allow for all γ > 1. Implosion singularities are much more severe than shocks in the sense that the density and velocity themselves blow up (in shock singularities, it is their gradient that blows up) in finite time at the center of symmetry. The proof depends on the detailed structure of the equations and in particular relies on a careful analysis of a phase portrait for self-similar solutions. Hence, it is of interest to decide whether a similar result holds for the relativistic Euler equations. We highlight the important related work [47], in which the authors studied spherically symmetric solutions the Einstein-Euler equations under adiabatic equations of state p = ερ for 0 < ε ≪ 1 sufficiently small and proved the existence of self-similar solutions that form a naked singularity. See also [46], in which the authors studied solutions to the non-relativistic Euler-Poisson system under adiabatic equations of state p = ρ γ for γ ∈ (1, 4 3 ) and showed the existence of spherically symmetric self-similar imploding solutions modeling gravitational collapse, i.e., initially smooth solutions such that the density blows up in finite time. (9) Inviscid limits. In [25], the authors studied the inviscid limit of solutions to the 1D viscous Burgers' equation all the way up to the time of first shock formation in the inviscid solution. They decomposed the viscous solution into a singular piece and a smoother piece and proved that the viscous solution converges to the singular piece in L ∞ as the viscosity vanishes, where the L ∞ norm is taken over the entire slab of classical existence of the inviscid solution. This is the first result of its type that extends all the way to the time of first singularity formation. Important open problems include extending this result to the 1D compressible Euler equations (where the corresponding viscous equations are the Navier-Stokes equations) and, after that, to multi-dimensions. It would also be of interest to extend the result to the relativistic Euler equations, but as of present, it is not clear if there exist any well-posed relativistic viscous fluid models that suppress the singularity formation while retaining physically desirable features such as causality. (10) Rarefaction waves. In Sect. 5.7, we described the recent works [68,69] on irrotational rarefaction wave solutions to the 2D compressible Euler equations. It is of interest to eliminate the irrotationality assumption and to extend these results to the relativistic Euler equations, and perhaps even the Einstein-Euler equations.