The emergence of the singular boundary from the crease in 3D compressible Euler flow
LA,JS]Leonardo Abbrescia$^{* \dagger}$ and Jared Speck$^{** \dagger\dagger}$

$^{*}$Vanderbilt University, Nashville, TN, USA.

$^{**}$Vanderbilt University, Nashville, TN, USA.

$^\dagger$ LA gratefully acknowledges support from an NSF Postdoctoral Fellowship.

$^{\dagger\dagger}$JS gratefully acknowledges support from NSF grant # DMS-2054184 and NSF CAREER grant # DMS-1914537.

We study the Cauchy problem for the $3D$ compressible Euler equations under an
arbitrary equation of state with positive speed of sound,
aside from that of a Chaplygin gas.
For open sets of smooth initial data with non-trivial vorticity and entropy,
our main results yield a constructive proof of the formation, structure, and stability of the singular boundary,
which is the set of points where the solution forms a shock singularity, i.e., where some
first-order Cartesian coordinate partial derivatives of the velocity and density blow up.
We prove that in the solution regime under study,
the singular boundary has the structure of a degenerate, acoustically null
$3D$ submanifold-with-boundary. Our approach yields the full structure of a neighborhood
of a connected component of the crease, which is a $2D$
acoustically spacelike
submanifold equal to the past boundary of the singular boundary.
In the study of shocks, the crease plays the role of the “true initial singularity”
from which the singular boundary emerges,
and it is a crucial ingredient for setting up the shock development problem.
These are the first results revealing the totality of these structures
without symmetry, irrotationality, or isentropicity assumptions. Moreover,
even within the sub-class of irrotational and isentropic solutions,
these are the first constructive results revealing these structures
without a strict convexity assumption on the shape of the singular boundary.
Our proof relies on a new method: the construction of rough foliations
of spacetime, dynamically adapted to the exact shape of the singular boundary
and crease, where the latter is provably two degrees less differentiable than the fluid.
Our results also set the stage for our forthcoming paper,
in which we will prove the emergence and stability of a Cauchy horizon,
which emanates from the crease and “evolves” in a direction that is “opposite”
the singular boundary in a sense determined by the intrinsic acoustic geometry of the flow.

Cauchy horizon;
characteristic current;
compressible Euler equations;
eikonal function;
maximal development;
null condition;
null hypersurface;
null structure;
shock development problem;
shock formation;
singular boundary;
stable singularity formation;
wave breaking;
vectorfield method

Mathematics Subject Classification (2010)
Primary: 35L67 - Secondary: 35L05, 35Q31, 74J40, 76N10

§ INTRODUCTION

This is the first of two papers in which we construct a large (though compact) portion of
the maximal classical globally hyperbolic development
(which we refer to as the “maximal development” for short from now on)
– up to the boundary –
of the initial data for open sets of initially smooth, shock-forming solutions
to the $3D$ compressible Euler equations with vorticity and dynamic entropy,
and without symmetry assumptions.
Roughly speaking, the maximal classical development
is the largest possible classical solution
region that is determined by the initial data.
Our main results are the constructions of the (singularity-forming) solution
$+$ localized region depicted in Fig. <ref>
(see Remark <ref> for comments on our notation in the figure)
and a proof of their stability under Sobolev-class perturbations of the initial data on
the Cauchy hypersurface $\Sigma_0 \eqdef \lbrace t = 0 \rbrace$.

Our papers provide the first results that
construct and fully justify Fig. <ref> for open sets of initial data
without symmetry, irrotationality, or isentropicity assumptions.

The present paper concerns the analysis up to the singular boundary,
which we denote by[More precisely,
we follow the solution up to a compact portion of the singular boundary that we denote by
$\mathcal{B}^{[0,\muxmulevelsetvalue_0]}$ in our main theorems. ]
while our companion work <cit.> concerns the analysis up to the Cauchy horizon,
which we denote by “$\underline{\mathcal{C}}$.”
Roughly, $\mathcal{B}$ is the submanifold-with-boundary[It is not obvious that the set of blowup-points
has the structure of a submanifold-with-boundary in Cartesian coordinate space.
Indeed, uncovering this structure is one of the main results of
the paper. While this structure holds for open sets of solutions, including the solutions we handle
in this paper, for other solutions, the set of blowup-points might fail to have the structure of a
of points where the fluid's first derivatives blowup
(though the fluid variables themselves remain bounded),
while $\underline{\mathcal{C}}$ is a future boundary that “feels the influence” of
the singularity, even though the solution remains classical along $\underline{\mathcal{C}} \backslash \partial_- \mathcal{B}$.
Readers can jump to Theorem <ref>
for an abbreviated version of the main results of this paper,
and Theorems <ref> and <ref>
for precise, extended statements.
As we describe in Sect. <ref>, the papers resolve several open problems
and allow one to properly set up the shock development problem,
which is the problem of (locally) describing the transition of the solution
from classical to weak, past the “initial singularity,”
which is the acoustically[The word “acoustic” refers to
the acoustical metric $\gfour$ of definition (<ref>),
which is the solution-dependent Lorentzian metric that dictates the geometry of sound waves.
Throughout the paper, all Lorentzian geometric notions such as spacelike,
timelike, null, etc. are with respect to $\gfour$. ]
spacelike co-dimension $2$
submanifold that we denote by “$\partial_- \mathcal{B}$”
(the past boundary of $\mathcal{B}$) in the figure.
In full generality, the shock development problem is open, though there has been
inspiring progress, which we describe in Sect. <ref>.
It is crucially important, for example in setting up the shock development problem,
that our approach yields a complete description of the
initial singularity. We highlight the following key point:

The “initial singularity” should not be thought
of as a point in spacetime, but rather the set
$\partial_- \mathcal{B}$ depicted in Fig. <ref>, which we refer to
as the crease, and which we prove has the structure of a $2$-dimensional acoustically spacelike
submanifold in the solution regime under study.

[scale=.36, grid = false, tics=5, trim=-.5cm 0cm -1cm -.5cm, clip]MinterestingCartesian.pdf
(50,7) $\Sigma_0$
(4,20) $(x^2,x^3) \in \mathbb{T}^2$
(2,34) $t$
(16,27) $x^1 \in \mathbb{R}$
(39,57) $\mbox{\upshape Crease}$
(46,53) 90$=$
(42,49) $\partial_- \mathcal{B}$
(60,44) $\Lunit$
(59,72.5) $\mathcal{B}$
(23,68.5) $\underline{\mathcal{C}}$

A localized subset of the maximal classical development in Cartesian coordinate space

(50,8) $\Sigma_0$
(1.2,22) $(x^2,x^3) \in \mathbb{T}^2$
(0,36) $t$
(41,32) $\nullhyparg{u}$
(14,28) $x^1 \in \mathbb{R}$
(68,80) $\mathcal{B}$
(23,74) $\underline{\mathcal{C}}$

Infinite density of the characteristics $\nullhyparg{u}$ on $\mathcal{B}$

Cartesian coordinate space illustrations of the main results

While we treat in detail a specific regime in our two papers
– perturbations of simple isentropic plane-symmetric solutions –
the methods we develop are robust and could be applied to other regimes,
such as perturbations of non-vacuum steady state solutions in $\mathbb{R}^{1+3}$.
We already stress that our results apply to situations in
which $\mathcal{B}$ fails to be strictly convex
(that is, whenever it fails to be strictly concave up),
and that substantial new ideas are needed in this case to follow the solution up to $\mathcal{B}$
and to understand its structure;
in Fig. <ref> we have depicted
a $\mathcal{B}$ that, while “convex in the $x^1$-direction,” it
fails to be convex “in the $(x^2,x^3)$-directions.”

§.§ Outline of the remainder of the introduction

In Sect. <ref>, we introduce the compressible Euler equations,
though not in the form we use to prove our main results.
In Sect. <ref>, we highlight some the main
challenges in the proofs of our main results and
discuss a few of the most important new ideas we use to overcome them.
In Sect. <ref>, we provide an overview of the two open problems
that we resolve in the paper, and we describe how they are connected to the shock development problem.
In Sect. <ref>, we state Theorem <ref>,
which is a first, somewhat informal, abbreviated version of our main results; see
Theorems <ref> and <ref>
for the extended, precise statements.
In Sect. <ref>,
having stated Theorem <ref>,
we provide some extended remarks on the ideas and methods we use in the proof.
In Sect. <ref>,
we describe the most relevant precursor results to this paper,
focusing on those works that yielded methods that we use here.
In Sect. <ref>, we discuss the history of the study of shock formation
and highlight some important developments in the subject.
In Sect. <ref>, we provide an overview of
the main ideas of the proofs of our main results, in particular highlighting
various technical issues that were not discussed earlier in the introduction.
In Sect. <ref>, we provide an outline of the remainder of the paper.

§.§ First version of the equations and local well-posedness

Our main results concern the Cauchy problem for the $3D$ compressible Euler equations,
which can be formulated as a quasilinear hyperbolic PDE system
in the velocity $v : \mathbb{R} \times \Sigma \rightarrow \mathbb{R}^3$,
the density $\varrho : \mathbb{R} \times \Sigma \rightarrow [0,\infty)$,
and the entropy $\Ent : \mathbb{R} \times \Sigma \rightarrow \mathbb{R}$.
In this paper, $\Sigma \eqdef \mathbb{R} \times \mathbb{T}^2$ denotes the “space manifold,”
that is, we assume that “space” is diffeomorphic to $\mathbb{R} \times \mathbb{T}^2$.
In our setup, the Cartesian coordinates on spacetime are $(t,x^1,x^2,x^3)$,
where $t$ is the standard Cartesian time function,
$x^1$ denotes the standard Cartesian coordinate on $\mathbb{R}$, and $(x^2,x^3)$ denote
standard Cartesian coordinates on $\mathbb{T}^2 \eqdef [-\pi,\pi]^2$ (with the endpoints identified).
The spatial topology $\mathbb{R} \times \mathbb{T}^2$
allows us to simplify various aspects of our approach, leading to a cleaner
presentation of the analysis. However, our analysis
is local in spacetime and, with modest additional effort,
all our results could be extended to other spatial topologies such
as $\mathbb{R}^3$, without cut-off functions or partitions of unity.
Relative to the Cartesian coordinates $(t,x^1,x^2,x^3)$,
the $3D$ compressible
equations can be expressed as follows[Throughout, if $\mathbf{V}$ is a vectorfield and $f$ is a scalar function,
then $\mathbf{V} f \eqdef \mathbf{V}^{\alpha} \partial_{\alpha}f$ denotes the derivative of $f$
in the direction of $\mathbf{V}$. ]
(see <cit.> for background on the equations):

\begin{align}
\Transport v^i
& =
- \frac{\partial_i p}{\varrho},
&& (i=1,2,3),
\label{E:INTROTRANSPORTVI}
\\
\label{E:INTROTRANSPORTDENSITY}
\Transport \varrho
& = - \varrho \Flatdiv v,
\\
\Transport \Ent
& = 0,
\label{E:INTROBS}
\end{align}

$p$ is the pressure,
$\Transport$ denotes the material vectorfield:

\begin{align}\label{E:MATERIALDERIVATIVEVECOTRFIELD}
\Transport
& \eqdef \p_t + v^a \p_a,
\end{align}

and $\Flatdiv$ is the standard Euclidean divergence operator (see Def. <ref>).
To close the equations, we assume an equation of state $p = p(\varrho,\Ent)$.
Our results apply for any sufficiently smooth equation of state
– except for that of the Chaplygin gas (see Sect. <ref>) –
with positive sound speed $\Speed$ defined by:

\begin{align} \label{E:SOUNDSPEED}
\Speed
& \eqdef \sqrt{p_{;\varrho}},
\end{align}

where $p;_{\varrho}$ is the partial derivative of the equation of state with respect
to the density at fixed entropy.

We assume that smooth initial data
for (<ref>)–(<ref>)
are prescribed along the spacelike hypersurface
$\Sigma_0 \eqdef \lbrace t = 0 \rbrace = \lbrace 0 \rbrace \times \mathbb{R} \times \mathbb{T}^2$.
We consider only initial data such that $\varrho|_{\Sigma_0} > 0$, thereby avoiding the severe
degeneracies that can occur at fluid-vacuum boundaries.
It is well-known that the equations (<ref>)–(<ref>)
are locally well-posed for
non-vacuum initial data on $\Sigma_0 \eqdef \lbrace t= 0 \rbrace$
such that $(v,\varrho,\Ent)|_{\Sigma_0} \in H^3(\Sigma_0)$.
To follow the solution to the singular boundary, we assume that the data belong
to a sufficiently high order Sobolev space, where different solution variables
have distinct, directionally dependent amounts of regularity;
see Sect. <ref> for our detailed assumptions.
All of the known approaches to studying shocks away from symmetry rely on
the assumption that the data belong to a high order
Sobolev space. This is due to possibly singular energy estimates at the high
derivative levels – even in the “good” geometric
coordinate system $(t,u,x^2,x^3)$, described below, which “hides” the singularity
at the low to mid derivative levels; see Sects. <ref>
and (<ref>)
for further discussion of these fundamental technical issues.

We stress up front that our analysis crucially relies on a
geometric reformulation of (<ref>)–(<ref>)
as a system of covariant wave equations coupled to transport-div-curl equations,
where the nonlinear terms exhibit remarkable null and regularity properties.
The reformulation was derived in <cit.> (see also the precursor <cit.>),
and we recall it in Theorem <ref>.

§.§ An overview of the degeneracies and difficulties in the problem

The analysis needed to fully justify Fig. <ref> is fraught with degeneracies and difficulties.
While many prior works on shocks have constructed the solution in strict subsets
of the region in Fig. <ref>, our papers are the first
to fully grapple with the degeneracies and
construct the solution in the entire region.

To handle solutions with non-zero vorticity and dynamic entropy,
we rely on an arsenal of geometric and analytic techniques, developed in
earlier works dC2007,jLjS2018,jS2019c,jLjS2020a,lAjS2020,
combined with key new ideas that we describe below.
Here, we highlight some of the main challenges in the analysis and
mention some of the methods we use to overcome them.

* (Singularities). The solution forms a shock along
the submanifold-with-boundary $\mathcal{B}$, i.e., some of its first-order partial derivatives with respect
to the Cartesian coordinates blow up, though the solution itself remains bounded.
The blowup-dynamics are extremely rich: some quantities exhibit blowup of their derivatives
– but only derivatives in certain directions[Roughly, along $\mathcal{B}$, many quantities'
derivatives in directions transversal to the characteristics blow up,
while their tangential derivatives remain bounded, much like in the simple case of
Burgers' equation in $1D$: $\partial_t \Psi + \Psi \partial_x \Psi = 0$.] – while other quantities
and their derivatives in all directions remain bounded.
We refer to Theorem <ref> for the details.

* (Nonlinear geometric optics and geometric coordinates) As in many prior works on shock formation,
to follow the solution up to $\mathcal{B}$ and to obtain a precise understanding of the singularity formation,
we cannot rely on the Cartesian coordinates, which are not adapted to the singularity.
Instead, we rely on nonlinear geometric optics.
Specifically, nonlinear geometric optics yields a “geometric coordinate system”
(see Def. <ref>)
– one that is globally homeomorphic (in the compact region under study)
to the Cartesian coordinate system but not diffeomorphic to it along $\mathcal{B}$ –
relative to which the solution remains rather smooth; see Sect. <ref>.
The key ingredient in implementing nonlinear geometric optics
is an eikonal function $u$ solving the acoustic eikonal equation
$(\gfour^{-1})^{\alpha \beta} \partial_{\alpha} u \partial_{\beta} u = 0$,
where the acoustical metric $\gfour$ is the solution-dependent Lorentzian metric (see (<ref>))
that captures the intrinsic geometry of sound waves.
The level sets of $u$, which we denote by $\nullhyparg{u}$, are characteristic for the compressible Euler equations;
roughly, the $\nullhyparg{u}$ represent surfaces along which sound waves can propagate.
In the present paper, our geometric coordinate system is $(t,u,x^2,x^3)$, where $u$ is the
eikonal function and $t,x^2,x^3$ are standard Cartesian coordinate functions.

* (Degeneracies in the acoustic geometry)
The fluid singularity along $\mathcal{B}$ is intimately tied to degeneracies in the acoustic geometry;
roughly, along $\mathcal{B}$, the characteristic hypersurfaces $\nullhyparg{u}$
(which we also refer to as “characteristics,”
“null hypersurfaces,” “acoustically null hypersurfaces,” or “$\gfour$-null hypersurfaces”),
develop infinite density along $\mathcal{B}$. In the present paper,
the infinite density of the $\nullhyparg{u}$ is characterized
by the vanishing of a function $\upmu$, the inverse foliation density, which we describe below in detail;
$\upmu$ is positive in the maximal development, except along $\mathcal{B}$.
We depict this infinite density in Fig. <ref>, where we show
three distinct null hypersurfaces “piling up” along $\mathcal{B}$.
It is important to appreciate that in the solution regime we are studying,
distinct characteristic hypersurfaces, viewed as submanifolds in Cartesian coordinate space,
never actually intersect[This fact follows from the homeomorphism property
of the map $\Upsilon$, as described in Theorem <ref>.]
on $\mathcal{B}$, even though their density becomes infinite.
This phenomenon is crucial for properly setting up the shock development problem.

[scale=.6, grid = false, tics=5, trim=-.5cm -1cm -1cm -.5cm, clip]CHcutsoffxmupositive.pdf
(73,68) $\mathcal{B}$
(52,55.5) $\partial_-\mathcal{B}$
(18,66) $\underline{\mathcal{C}}$
(48,10) $\Sigma_0$
(-3,15.5) $(x^2,x^3) \in \mathbb{T}^2$
(6,27) $t$
(-5,22.5) $u \in \mathbb{R}$

The singular boundary and Cauchy horizon in geometric coordinates

The singular boundary $\mathcal{B}$, viewed as a subset of geometric coordinate space
(see Fig. <ref>), is:

\begin{align} \label{E:INTROSINGBOUNDARY}
\mathcal{B}
& = \lbrace \upmu = 0 \rbrace
\cap
\lbrace \muX \upmu \leq 0 \rbrace,
\end{align}

where the vectorfield $\muX$, described later on
(see Fig. <ref> for a depiction of $\muX$ in Cartesian coordinate space),
is transversal to the $\nullhyparg{u}$ and satisfies $\muX u = 1$.
The complementary set
\lbrace \upmu = 0 \rbrace
\cap
\lbrace \muX \upmu > 0 \rbrace
is not part of the singular boundary.
In the context of Fig. <ref> in geometric coordinate space,
the irrelevance of
\lbrace \upmu = 0 \rbrace
\cap
\lbrace \muX \upmu > 0 \rbrace
for $\mathcal{B}$ and the maximal development can be understood as follows:
the singular boundary portion
$\mathcal{B}$ in the figure cannot be extended to the left into the region where
$\muX \upmu$ would be positive
because before that region has a chance to dynamically develop,
it will be cut off by a Cauchy horizon
emanating from the crease $\partial_- \mathcal{B}$,
which is the left boundary of $\mathcal{B}$ in the figure.
The crease $\partial_- \mathcal{B}$ is characterized by:

\begin{align} \label{E:INTROCREASE}
\partial_- \mathcal{B}
& = \lbrace \upmu = 0 \rbrace
\cap
\lbrace \muX \upmu = 0 \rbrace.
\end{align}

From (<ref>)–(<ref>), the status of $\partial_- \mathcal{B}$
as a boundary of $\mathcal{B}$ is clear.
In the regime under study, the two level sets on RHS (<ref>)
intersect transversally, which is what gives the crease the structure of a $2D$ submanifold.
The transversality of the intersection is a consequence of acoustical transversal convexity,
which is mild condition satisfied by open sets of solutions and which we describe below in more detail.
By doing formal Taylor expansions starting from the crease,
one can check that the unattainable portion $\{\upmu = 0\} \cap \{ \muX \upmu >0\}$
would indeed lie to the causal future of the Cauchy horizon $\underline{\mathcal{C}}$;
see the dotted portion in Fig. <ref>,
which formally depicts $\{\upmu = 0\} \cap \{ \muX \upmu >0\}$.

* (Null hypersurfaces and PDE energy degeneracies). In Fig. <ref>,
$\mathcal{B}$ and $\underline{\mathcal{C}}$ are acoustically null hypersurfaces emanating from the crease,
and in particular, $\mathcal{B}$ is ruled by acoustically null curves whose tangent vectors
are denoted by $\Lunit$ in the figure; see Prop. <ref> for a
proof of these properties of $\mathcal{B}$, and see Remarks <ref> and
Remark <ref>
for a discussion of some degeneracies that occur along $\mathcal{B}$.
As is well-known, any $L^2$-type energy that one uses to control solutions
necessarily degenerates along null hypersurfaces,
becoming only positive semi-definite instead of positive definite.
This is a particularly challenging issue in the present context,
where singularities are forming along all of $\mathcal{B}$.

* (Regularity and rough foliations).
Many objects in the construction
(in particular, the crease $\partial_- \mathcal{B}$)
are less regular than the fluid solution, which leads to difficult regularity theory for the problem.
To “detect” these rough objects as they emerge in the course of the evolution,
we rely on a new family of rough foliations
given by the level sets of rough time functions, described below.
The word “rough” refers to the fact that the rough time functions are also less regular than the fluid solution.
There is another way in which the problem of shock formation can be viewed as a low regularity problem:
the piling up of the characteristics is tied to the blowup of
a Euclidean-unit-length derivative of various fluid variables in directions transversal
to the characteristics, even though the solution remains rather smooth in directions tangent
to the characteristics. In particular,
we are forced to close the estimates knowing that with respect
to the Cartesian differential structure, there will be no differentiability
in transversal directions at the end of the classical evolution.
As we already mentioned, the geometric coordinates $(t,u,x^2,x^3)$
partially ameliorate this difficulty in the sense
that the solution remains rather smooth with respect to them.
Nonetheless, as in other works on shock formation,
our high order geometric energies can still become singular;
see Sect. <ref>.
This is one of the main technical challenges in the PDE analysis
since singular high order energy estimates make it difficult for us to prove that
the solution's partial derivatives with respect to the geometric coordinates
remain bounded at the lower derivative levels.

* (One submanifold of $\mathcal{B}$ at a time via a family of rough time functions).
Our approach to constructing $\mathcal{B}$ is to show that it can be foliated by
a family of $\muxmulevelsetvalue$-parameterized submanifolds $\twoargmumuxtorus{0}{-\muxmulevelsetvalue}$
with $\muxmulevelsetvalue \geq 0$ a real parameter,
and to construct each $\twoargmumuxtorus{0}{-\muxmulevelsetvalue}$, “one $\muxmulevelsetvalue$ at a time;”
see Fig. <ref>. Our construction is such that the crease $\partial_- \mathcal{B}$
coincides with $\twoargmumuxtorus{0}{0}$.
One might wonder why we didn't try to derive all of $\mathcal{B}$ “at the same time.”
From the discussion three points above, we see that
that approach would have effectively required us to work with foliations of spacetime that contain
or are asymptotic to level sets of $\upmu$ in regions where $\upmu$ is small
(recall that $\mathcal{B}$ is a portion of $\lbrace \upmu = 0 \rbrace$).
The difficulty is that for real numbers $\mulevelsetvalue$ small and positive, near the crease,
the level sets $\lbrace \upmu = \mulevelsetvalue \rbrace$ have
a $\gfour$-timelike portion,
along which top-order $L^2$ estimates for the solution are not available.
This can formally be understood as the statement that
in the solution regime under study, along the surface $\lbrace \upmu = 0 \rbrace$,
RHS (<ref>) would be positive[In the regime under study,
the term $\Lunit \upmu$ on RHS (<ref>) is strictly negative
and the term $|\angD \upmu|_{\gtorus}^2$ is of negligible size.]
in the regions
where $\muX \upmu > 0$.
See Fig. <ref>,
where these level sets become $\gfour$-timelike in the region near their intersection with
the Cauchy horizon $\underline{\mathcal{C}}$.
In fact, for $\mulevelsetvalue > 0$ small,
any surface that agreed with the level set $\lbrace \upmu = \mulevelsetvalue \rbrace$ up to
second order along the surface $\lbrace \muX \upmu = 0 \rbrace$
would suffer from the same difficulty:
it would necessarily contain a $\gfour$-timelike portion,
along which $L^2$ estimates for the solution are not available.

The upshot is that in the solution regime under study,
it is impossible to detect the entire singular boundary $\mathcal{B}$
by deriving estimates on $C^2$ (relative to the geometric coordinates)
$\gfour$-spacelike foliations of spacetime
such that $\mathcal{B}$ (including its past boundary $\partial_- \mathcal{B}$)
is contained in the interior of one of the leaves of the foliation.
This difficulty is connected to the following issue: in our main theorem,
namely Theorem <ref>,
the region $\MInteresting$
that we study (which contains $\mathcal{B}$) is foliated by the level sets of a $C^{1,1}$
(relative to the geometric coordinates) time function $\newtimefunction$,
and this $C^{1,1}$ regularity is optimal given the shape of region;
see Remark <ref>.
In our approach to the PDE analysis, the $C^{1,1}$ regularity of $\newtimefunction$
would be insufficient for our proofs of some of our estimates
(e.g., the co-dimension two Gauss curvature estimates we derive in Lemma <ref>).
As we explain below, we avoid these difficulties by avoiding working directly with
$\newtimefunction$; we instead work with a one-parameter family of time functions
$\timefunctionarg{\muxmulevelsetvalue}$, which are more regular than $\newtimefunction$.

[scale=.6, grid = false, tics=5, trim=-.5cm -1cm -1cm -.5cm, clip]levelsetsofmutimelikepiece.pdf
(73,74) $\mathcal{B}$
(51.5,60.5) $\partial_-\mathcal{B}$
(18.5,72.5) $\underline{\mathcal{C}}$
(48,10) $\Sigma_0$
(1,22) $(x^2,x^3) \in \mathbb{T}^2$
(6,32) $t$
(87,45) $\{\upmu = \mulevelsetvalue\}$
(80,30.5) $\{\upmu = 2 \mulevelsetvalue\}$
(49,22) $\{\muX \upmu = 0\}$
(62,31) $\{\muX \upmu < 0\}$
(30,31) $\{\muX \upmu > 0\}$
(38,24) $\muX$
(-5,28) $u \in \mathbb{R}$

The singular boundary and the level sets $\{\upmu = \mulevelsetvalue\}$
and $\{ \upmu = 2\mulevelsetvalue\}$ for a small $\mulevelsetvalue > 0$, in geometric coordinates

To construct the submanifolds $\twoargmumuxtorus{0}{-\muxmulevelsetvalue} \subset \mathcal{B}$ and circumvent
the difficulties noted in Remark <ref>,
we proceed as follows.
For each $\muxmulevelsetvalue \in [0,\muxmulevelsetvalue_0]$, where $\muxmulevelsetvalue_0 > 0$ is a small, data-dependent constant,
we construct a rough time function $\timefunctionarg{\muxmulevelsetvalue}$,
which is $C^{2,1}$ relative to the geometric coordinates;
see Sect. <ref>.
Note that $\timefunctionarg{\muxmulevelsetvalue}$ is one degree more differentiable
than the time function $\newtimefunction$ mentioned above, which is sufficient
for avoiding the difficulties we highlighted in Remark <ref>.
Our construction of $\timefunctionarg{\muxmulevelsetvalue}$
depends on the eikonal function $u$, i.e., our construction of
$\timefunctionarg{\muxmulevelsetvalue}$ relies on nonlinear geometric optics.
For each $\muxmulevelsetvalue \in [0,\muxmulevelsetvalue_0]$, the level sets of $\timefunctionarg{\muxmulevelsetvalue}$
yield a foliation of spacetime by $\gfour$-spacelike hypersurfaces, and
$\twoargmumuxtorus{0}{-\muxmulevelsetvalue} \subset \lbrace \timefunctionarg{\muxmulevelsetvalue} = 0 \rbrace$.
Hence, to construct $\twoargmumuxtorus{0}{-\muxmulevelsetvalue}$, it suffices to control the fluid
up to the level set $\lbrace \timefunctionarg{\muxmulevelsetvalue} = 0 \rbrace$; the vast majority of
our efforts in this paper are dedicated towards that task.
The union $\bigcup_{\muxmulevelsetvalue \in [0,\muxmulevelsetvalue_0]} \twoargmumuxtorus{0}{-\muxmulevelsetvalue}$ is the portion of $\mathcal{B}$
that we construct in our main theorem.
The portion of $\mathcal{B}$ that formally corresponds to $\muxmulevelsetvalue < 0$
never has a chance to emerge in the maximal classical development because it is cut off by the
Cauchy horizon $\underline{\mathcal{C}}$. In Fig. <ref>,
the “irrelevant portion” of $\mathcal{B}$, corresponding to $\muxmulevelsetvalue < 0$,
is formally delineated by dotted curves
(we cannot actually construct this portion, and we have displayed it only to illustrate that
$\underline{\mathcal{C}}$ lies below it).
We could have extended our results to handle a larger range of $\muxmulevelsetvalue$ values, i.e., $\muxmulevelsetvalue > \muxmulevelsetvalue_0$;
we avoided this because we would have had to modify our construction of the $\timefunctionarg{\muxmulevelsetvalue}$ for large
$\muxmulevelsetvalue$, which would have lengthened the paper.
We refer Sect. <ref> for a
more detailed overview of our construction of the rough time functions and the $\twoargmumuxtorus{0}{-\muxmulevelsetvalue}$.

The word “rough” in “rough time function”
refers to the fact that the elements of
$\lbrace \timefunctionarg{\muxmulevelsetvalue} \rbrace_{\muxmulevelsetvalue \in [0,\muxmulevelsetvalue_0]}$
are less regular than the fluid; see Sect. <ref>
for further details.
In the present paper, each $\timefunctionarg{\muxmulevelsetvalue}$ is $C^{2,1}$
(with respect to the geometric coordinates)
because the fluid variables are $C^{3,1}$.
The fluid variables are $C^{3,1}$ because
we have only assumed limited differentiability on the fluid data in directions
transversal to the characteristics $\nullhyparg{u}$,
even though we assumed they are much smoother in the $\nullhyparg{u}$-tangential directions;
see Sect. <ref> for our data assumptions.
Despite the terminology “rough time function,” if
we had instead assumed that the fluid data were $C^{\infty}$, then
each $\timefunctionarg{\muxmulevelsetvalue}$ would also have been $C^{\infty}$.
In contrast, even with $C^{\infty}$ fluid data, the time function $\newtimefunction$ from
Remark <ref> would have had only $C^{1,1}$ regularity.

* (Lack of strict convexity).
Observe that in Fig. <ref>,
there are points $q \in \mathcal{B}$ such that the tangent plane
to $\mathcal{B}$ at $q$, which we denote by $T_q \mathcal{B}$,
does not lie below $\mathcal{B}$.
Moreover, there are points $q \in \partial_- \mathcal{B}$ such that
every (three-dimensional) causal plane containing the two-dimensional
subspace $T_q \partial_- \mathcal{B}$
(which appears to be one-dimensional in Fig. <ref>, due to our
suppression of a spatial dimension), fails to lie below $\mathcal{B}$
even locally near $q$.
Let us informally refer to these phenomena as “absence of strict convexity,”
where here, “convexity” informally refers to “upwards bending,”
and “strict convexity”
– though not featured Fig. <ref> – would refer to “upwards bending in all directions.”
We clarify that these notions implicitly refer to the $1+3$-dimensional Cartesian coordinate space
since the tangent planes referred to above are understood to be standard Cartesian-flat hyperplanes
and Fig. <ref> is a picture in Cartesian coordinate space.
The absence of strict convexity

poses serious technical difficulties.

In the absence of symmetry and strict convexity, the entirety of the crease $\partial_- \mathcal{B}$
has never before been fully constructed for any open set of shock-forming solutions to any hyperbolic PDE.
In particular, this aspect of our main results is new even in the case of irrotational and isentropic solutions.

The rough time function $\timefunctionarg{\muxmulevelsetvalue}$ and corresponding rough foliations allow us to derive,
through a fully constructive approach,
the structure of the singular boundary,
even if it is not strictly convex.
Instead of strict convexity, we rely on
acoustical transversal convexity
(which we refer to as transversal convexity for short),
which allows us to handle, for example, perturbations of symmetric solutions,
where strict convexity of the singular boundary
can fail due to the approximate symmetry.
Roughly, transversal convexity is a form of convexity only in
a particular direction, specifically in a direction that is transversal to the level sets of the eikonal function $u$.
Note that transversal convexity refers to the structure of $\mathcal{B}$ viewed as an embedded submanifold
of geometric coordinate space (as opposed to Cartesian coordinate space).
In the solution regime under study, transversal convexity is captured by our
data assumption (<ref>) on the inverse foliation density $\upmu$,
which we are able to propagate throughout the evolution (see (<ref>)).
The singular boundary $\mathcal{B}$ in Fig. <ref>
enjoys transversal convexity
(roughly, it is “parabolic in the $u$ direction” at fixed $(x^2,x^3)$ near $\partial_- \mathcal{B}$).
One could check that in the solution regime under study,
the transversal convexity of $\mathcal{B}$ in geometric coordinates also implies,
in the Cartesian coordinate picture,
the convexity of the $x^1$-parameterized curves in $\mathcal{B}$ along which $(x^2,x^3)$ are fixed;
while we do not directly need this “Cartesian transversal convexity” in our analysis,
we have exhibited the “upwards bending[In $1+1$ dimensions,
under transversal convexity,
the embedding of the singular boundary in Cartesian coordinate space
can be modeled by the $u$-parameterized curve $t = u^2$,
$x^1 = u^2 + u^3$, for $u \leq 0$.
Near the origin (which models the crease) in $(t,x^1)$-space,
this singular boundary-modeling curve is asymptotic to the graph of
$t = (x^1)^2 + (x^1)^{3/2}$ for $x^1$ small and positive, which bends upwards
(c.f. the singular boundary in Fig. <ref>)
and has regularity $C^{1,1/2}$ (c.f. the
regularity of $\Upsilon(\mathcal{B}^{[0,\muxmulevelsetvalue_0]})$ stated in
Theorem <ref>).
of $\mathcal{B}$ in the $x^1$-direction”
in Fig. <ref>.
Our assumption of transversal convexity is close to optimal in the sense that without it,
the qualitative character of the singular boundary can dramatically change,
even for plane-symmetric solutions; see Sect. <ref> for
further discussion.

We also highlight that the assumption of transversal convexity played a crucial role in Christodoulou's resolution
<cit.> of the restricted shock development problem, which we describe below.

* (Degenerate wave and elliptic estimates on curved domains).
To close the high order $L^2$ estimates, we must
adapt a variety of hyperbolic energy estimates and
“top-order” elliptic estimates for the vorticity, entropy, and geometry to the
precise shape of
$\partial_- \mathcal{B}$,
and $\underline{\mathcal{C}}$,
which are not known in advance. The shock singularity introduces degeneracies into these estimates,
and when we control the top-order derivatives of the vorticity and entropy using elliptic estimates,
our handling of these degeneracies
requires our observation of special cancellations
within delicately constructed “elliptic-hyperbolic” identities.
To exhibit the cancellations, we rely on the full nonlinear structure of the
geometric formulation of compressible Euler flow provided by Theorem <ref>.
In constructing the elliptic-hyperbolic identities (see Sect. <ref>),
we rely on the framework we developed in <cit.>. However,
for the purposes of the present paper, we had to substantially upgrade that framework to accommodate the structure of the singularity.
The key new object that we use to derive the elliptic-hyperbolic identities is a well-constructed characteristic current,
defined in Def. <ref>.

§.§ The two open problems that we resolve and connection to the shock development problem

In his breakthrough 2007 monograph <cit.>,
Christodoulou gave a sharp description of the stable formation of shock singularities,
starting from open sets of smooth initial data,
in solutions to the $3D$ irrotational and isentropic relativistic Euler equations.
Together with Miao, he later extended his results to the $3D$ compressible Euler equations <cit.>,
again for irrotational and isentropic solutions.
These results revealed a large subset of the maximal classical development,
including a portion of the boundary.

In the wake of <cit.>, there have been many exciting developments on
the formation of shock singularities and the subsequent evolution of the solution as a weak solution,
after the shock; see Sect. <ref>
for further discussion. However, two fundamental problems have remained open:

* (The full structure of bounded portions of the maximal classical development).
As is explained on <cit.>*Pages 929, 968–969, Christodoulou's approach
yields a union of developments of the initial data, where each of his developments
can be foliated by portions of Cartesian-flat[By a “Cartesian-flat” hypersurface, we mean a plane with
respect to the standard rectangular coordinates on $\mathbb{R}^{1+3}$.]
spacelike hypersurfaces and portions of characteristic hypersurfaces.
By varying the “angle of tilt” of these Cartesian-flat hypersurfaces and varying the initial data of the characteristic
one obtains (see <cit.>*Pages 929, 968–969) “a larger part” of the maximal development.
While this approach yields a sharp description of some portion of the maximal development,
the precise portion that it reveals is not made explicit through the construction.
Moreover, from Fig. <ref>, one can infer that
for some solutions, there are portions of the boundary (in particular, portions of $\partial_- \mathcal{B}$)
that are not accessible
through Cartesian-flat spacelike foliations.
This is connected to the lack of strict convexity of the singular boundary, as we discussed in
Sect. <ref>;
see also Fig. <ref> and, in Sect. <ref>,
our discussion of the points $b_1$ and $b_2$ featured in the figure.
Hence, the following problem is glaring:

Can one derive the full structure of the maximal development,
at least in some region of spacetime that includes a neighborhood of the boundary
that contains the crease?

This problem, while mathematically rich in itself, is important for two other fundamental reasons:
I) The breakthrough result <cit.>
shows that in general, one cannot ensure uniqueness of the maximal classical development
until one constructs it and proves that it enjoys some crucial structural qualitative properties.
For quasilinear hyperbolic PDEs, the question of uniqueness of the maximal classical development is global in nature.
However, the results we derive in the present article and the companion <cit.>
exhibit a localized version of the crucial property that the maximal development “lies on one side of its boundary.”
In <cit.>, the authors showed, roughly speaking, that if this property holds globally,
then the maximal classical development is unique.
II) As is explained in <cit.>,
the full structure of a neighborhood of the boundary of the maximal development
is an essential ingredient for properly setting up the
aforementioned shock development problem.

* (Removing the irrotationality and isentropicity assumptions).
Assuming a positive resolution to the first problem,
a second one of clear physical importance stands out:

Can one extend the result away from the irrotational and isentropic
class of solutions, that is, to handle solutions with non-zero vorticity and dynamic entropy?

It is of fundamental importance to understand such “general solutions” because vorticity and entropy
will form in the weak solution (see, e.g., <cit.>), i.e., “after the first shock,” even if the original
initial data are irrotational, isentropic, and $C^{\infty}$. That is, if one aims towards developing
a global-in-space-and-time theory that accommodates the formation of shocks and their subsequent interactions,
vorticity and entropy are an unavoidable aspect of the dynamics.

In the present paper and its companion <cit.>, we resolve both problems
for the $3D$ compressible Euler equations (<ref>)–(<ref>).
In the notation of Fig. <ref>,
starting from smooth data on $\Sigma_0$, we construct the classical solution
in the region lying in between $\Sigma_0$ and[This union is not disjoint, for the crease $\partial_- \mathcal{B}$
is a past boundary of both of the closed sets $\underline{\mathcal{C}}$ and $\mathcal{B}$.]
$\underline{\mathcal{C}} \cup \partial_- \mathcal{B} \cup \mathcal{B}$. More precisely, our results apply to
open sets of smooth data satisfying assumptions described below
(see, for example, Theorem <ref>),
and we control the solution in bounded regions
of spacetime that contain all[More precisely, we construct an entire connected component of $\partial_- \mathcal{B}$.]
of $\partial_- \mathcal{B}$
and a full neighborhood of it in $\underline{\mathcal{C}}$ and $\mathcal{B}$.
For convenience, in the present paper, we have studied the solution only in
a single, “spatially local” region.
However, our approach could be used as a building block to study the solution
across space, at least in regions where the solution exhibits the property of
acoustical transversal convexity, mentioned in Sect. <ref>.

[scale=.36, grid = false,trim=-.5cm -1cm -1cm -1.9cm, clip=true]singhypcreasecartesian.pdf
(36,11) $\Sigma_0$
(11,19) $(x^2,x^3) \in \mathbb{T}^2$
(10,29) $t$
(20,23.5) $x^1 \in \mathbb{R}$
(31,90) $\twoargmumuxtorus{0}{0}$
(34,85) 90$=$
(31,82) $\partial_- \mathcal{B}$
(54,93) $\mathcal{B}$
(12,89) $\hypthreearg{0}{[-\rightu,\leftu]}{0}$
(43,36.5) $\twoargmumuxtorus{-\timefunction}{0}$
(55,37) $\datahypfortimefunctionarg{0}$
(18,30) $\hypthreearg{\timefunction}{[- \rightu,\leftu]}{0}$
(24.5,65) $\Wtransarg{0}$
(31,52) $\Wtransarg{0}$

Rough foliations adapted to the crease

[scale=.36,grid=false,trim=-.5cm -1cm -1cm -.5cm, clip=true]singhypkappacartesian.pdf
(36,11) $\Sigma_0$
(11.5,19) $(x^2,x^3) \in \mathbb{T}^2$
(10,29) $t$
(20,23.5) $x^1 \in \mathbb{R}$
(39,88) $\twoargmumuxtorus{0}{-\muxmulevelsetvalue}$
(58,98) $\mathcal{B}$
(46,35) $\twoargmumuxtorus{-\timefunction}{-\muxmulevelsetvalue}$
(15,90) $\hypthreearg{0}{[-\rightu,\leftu]}{\muxmulevelsetvalue}$
(59,39) $\datahypfortimefunctionarg{-\muxmulevelsetvalue}$
(30,27.5) $\hypthreearg{\timefunction}{[- \rightu,\leftu]}{\muxmulevelsetvalue}$
(29,66) $\Wtransarg{\muxmulevelsetvalue}$
(34,52.5) $\Wtransarg{\muxmulevelsetvalue}$

Rough foliations adapted to a non-crease torus $\twoargmumuxtorus{0}{-\muxmulevelsetvalue} \subset \mathcal{B}$

Rough foliations adapted to the singular boundary, depicted in Cartesian coordinate space

In the present paper, we follow the solution up
to the singular portion of the boundary, denoted by $\mathcal{B}$
in Figs. <ref> and <ref>,
where some first derivatives of the density and velocity blow up.
We again emphasize that $\mathcal{B}$ contains its past boundary, $\partial_- \mathcal{B}$,
which is the aforementioned crease.
In <cit.>, we construct the Cauchy horizon, denoted by
$\underline{\mathcal{C}}$ in Fig. <ref>,
which is an acoustically null hypersurface such that no fluid singularity forms along it,
but it is nonetheless a boundary of the maximal classical development
because its causal past intersects the singularity.
As we mentioned earlier, our analysis in the present paper fundamentally relies on
a new family of acoustically spacelike rough foliations that are precisely and dynamically adapted
to the shape of $\mathcal{B}$, which is not known in advance.
The foliations are level sets of a one-parameter family of rough time functions
$\lbrace \timefunctionarg{\muxmulevelsetvalue} \rbrace_{\muxmulevelsetvalue \in [0,\muxmulevelsetvalue_0]}$,
where $\muxmulevelsetvalue_0 > 0$ is a constant depending on the initial data.
Each $\timefunctionarg{\muxmulevelsetvalue}$ is defined on a portion of the maximal classical development union its boundary
and has a range $[\timefunction_0,0]$ for some constant $\timefunction_0  < 0$ depending on the initial data.
Moreover, $\timefunctionarg{\muxmulevelsetvalue}$ has the crucial property that its zero level set $\lbrace \timefunctionarg{\muxmulevelsetvalue} = 0 \rbrace$
is tangent to $\mathcal{B}$ and intersects it in a submanifold. More precisely,
$\mathcal{B} \cap \lbrace \timefunctionarg{\muxmulevelsetvalue} = 0 \rbrace = \twoargmumuxtorus{0}{-\muxmulevelsetvalue}$
is a torus with spacetime co-dimension $2$ such that
the tori $\twoargmumuxtorus{0}{-\muxmulevelsetvalue}$ foliate a neighborhood of $\partial_- \mathcal{B}$ in $\mathcal{B}$,
i.e., $\bigcup_{\muxmulevelsetvalue \in [0,\muxmulevelsetvalue_0]} \twoargmumuxtorus{0}{-\muxmulevelsetvalue}$
is a neighborhood of $\partial_- \mathcal{B}$ in $\mathcal{B}$.
In particular, the crease $\partial_- \mathcal{B}$ is equal to
In Fig. <ref>, we exhibit two level sets of $\timefunctionarg{0}$,
where the top one contains the crease.
In Fig. <ref>, for some $\muxmulevelsetvalue > 0$,
we exhibit two level sets of $\timefunctionarg{\muxmulevelsetvalue}$,
where the top one contains the torus $\twoargmumuxtorus{0}{-\muxmulevelsetvalue}$.

We already highlight that the hypersurfaces
in Fig. <ref>
play a crucial role in our construction of
$\timefunctionarg{\muxmulevelsetvalue}$ respectively.
In particular, by construction,
solve transport equations with initial data
given on
$\datahypfortimefunctionarg{-\muxmulevelsetvalue}$ respectively;
see Sect. <ref>
for the details.
The surfaces
are transversal to $\mathcal{B}$
and intersect $\mathcal{B}$ in the tori
While the transversality of $\datahypfortimefunctionarg{0}$
to $\mathcal{B}$ is crucial for our analysis,
the causal structure of
(i.e., whether they are acoustically timelike, spacelike or null)
is not important because
we do not have to derive any energy estimates for the fluid
In particular,
can be acoustically timelike, spacelike, or null.

In both papers, our approach relies on giving a complete description of the dynamics that in particular shows
what blows up and what remains regular and yields a sharp description of the structure and
regularity of $\mathcal{B}$, $\partial_- \mathcal{B}$, and $\underline{\mathcal{C}}$.
To handle the presence of vorticity and entropy,
we develop modified versions of the integral identities that we discovered in <cit.>,
adapted here to the precise structure of $\mathcal{B}$.
More precisely, in the present paper,
we rely on a new family of “elliptic-hyperbolic”
integral identities that are adapted to the rough foliations and the structure of $\mathcal{B}$,
and we develop new analytic techniques to handle the following difficulty, which permeates the paper:
the singular boundary is ruled by curves that have acoustically null tangent vectors,[In the context of non-degenerate
Lorentzian metrics,
hypersurfaces that are ruled by curves with null tangent vectors are usually called “null hypersurfaces,”
and as is well-known, those tangent vectors are also normals to the surface. However, we avoid referring to the singular boundary
as a “null hypersurface” because of the following severe degeneracy, which is tied to the many analytical difficulties
we must handle in the PDE analysis: the singular boundary is
a portion of the zero level set of a function “$\upmu$” whose gradient with respect to the Cartesian coordinates
blows up along its zero level set; below we discuss $\upmu$ in detail.
Hence, in the Cartesian differential structure, the gradient one-form of $\upmu$,
which in non-degenerate contexts
can be viewed as the metric dual to its normal, is ill-defined. ]
which leads to severe degeneracies in the estimates,
especially in the top-order energy estimates for the vorticity and entropy.
See Prop. <ref> and Remark <ref>
for detailed information
on the structure of the singular boundary, including its causal structure and its
differential-topological properties as a subset of Cartesian coordinate space.

[scale=.4, grid = false, tics=5, trim=-.5cm -1cm -1cm -.5cm, clip]weakshockhyp.pdf
(50,15) $\Sigma_0$
(4,26) $(x^2,x^3) \in \mathbb{T}^2$
(3,38) $t$
(32,15) $x^1 \in\mathbb{R}$
(35,31) $\partial_- \mathcal{B}$
(66,75) $\mathcal{B}$
(66,32) $\mathcal{K}$
(21,65) $\underline{\mathcal{C}}$
(56,76) $\Lunit$

A localized subset of the maximal classical development and the shock hypersurface in Cartesian coordinate space

In the solution regime that we treat in our main results,
the crease $\partial_- \mathcal{B}$ has the structure of a co-dimension $2$
spacelike submanifold; see Fig. <ref>.
This co-dimension $2$ spacelike submanifold structure
is essential for properly setting up the shock development problem using
the known techniques <cit.>.
In fact, without this structure, it is not even clear whether the shock development problem is well-posed.
In particular, the shock hypersurface (with boundary), denoted by “$\mathcal{K}$” in
Fig. <ref>, emanates from the crease
and starts off tangent to the singular boundary. Roughly, $\mathcal{K}$ is the hypersurface
of discontinuity for the weak solution that develops to the future of the crease.
$\mathcal{K}$ is not part of the constructions we provide in this paper;
instead, its construction would be part of the resolution of the shock development problem.
Let us again highlight one of the key difficulties in the problem:
$\mathcal{B}$ is ruled by acoustically null curves with tangent vectors
denoted by “$\Lunit$” in Fig. <ref>.
This leads to rather severe degeneracies
in the analysis because, as is well-known, the coercive quantities that can be used to control the solution
become degenerate along null hypersurfaces.
Here, the notion of null is with respect to the acoustical metric
$\gfour$, the solution-dependent Lorentzian metric (see (<ref>))
that captures the intrinsic geometry of sound waves in
the flow. See below for extended discussion of these fundamental issues.

§.§ Abbreviated statement of the main results

The full statement of our main results is quite lengthy,
due to the intricate geometric structures
and the highly tensorial nature of the singularity.
In Theorem <ref>
we provide an abbreviated, slightly informal statement of the main results.
In Theorems <ref> and <ref>,
we provide full statements of the main results.

We first provide two pictures illustrating the region that we study in the theorem.
In Fig. <ref>, we display the region
in geometric coordinates $(t,u,x^2,x^3)$. In Fig. <ref>,
we display the region in Cartesian coordinates $(t,x^1,x^2,x^3)$,
where in the labels, $\Upsilon(t,u,x^2,x^3) = (t,x^1,x^2,x^3)$
is the change of variables map between the two coordinate systems.

In most of the article, we consider sets such as the singular boundary $\mathcal{B}$
to be subsets of geometric coordinate space $\mathbb{R}_t \times \mathbb{R}_u \times \mathbb{T}^2$,
and we denote the image of these sets in Cartesian coordinate space $\mathbb{R}_t \times \mathbb{R}_{x^1} \times \mathbb{T}^2$
by explicitly indicating the change of variables map $\Upsilon$, e.g., by
We use this notation in particular in Fig. <ref>.
However, in many of the other figures that depict regions in Cartesian coordinate space,
such as Figs. <ref> and <ref>,
we have suppressed the map $\Upsilon$ so as to not clutter the figure.

[scale=.36, grid = false, tics=5, trim=-.5cm -1cm -1cm -.5cm, clip]Minterestingmainresults.pdf
(40,10) $\Sigma_0$
(61,18.5) $(x^2,x^3) \in \mathbb{T}^2$
(78,28) $t$
(62.5,25) $u \in\mathbb{R}$
(15,50) $\MLeft$
(45,53) $\MSingular$
(65,70) $\MRight$
(30,25) $\inthyp{\timefunction_0}{[- \rightu,\leftu]}$
(30,100) $\inthyp{0}{[- \rightu,\leftu]}$
(32,35) $\twoargmumuxtorus{-\timefunction}{0}$
(61,45) $\twoargmumuxtorus{-\timefunction}{-\muxmulevelsetvalue_0}$
(26,81) $\twoargmumuxtorus{0}{0}$
(49,89.5) $\twoargmumuxtorus{0}{-\muxmulevelsetvalue_0}$
(34,80) $\mathcal{B}^{[0,\muxmulevelsetvalue_0]}$
(42,35) $\datahypfortimefunctiontwoarg{0}{[\timefunction_0,0]}$
(66,49) $\datahypfortimefunctiontwoarg{-\muxmulevelsetvalue_0}{[\timefunction_0,0]}$
(17,65) $\Wtransarg{0}$
(28,52) $\Wtransarg{0}$
(15,53) $\Wtransarg{0}$
(25,40) $\Wtransarg{0}$
(43,66) $\Wtransarg{\muxmulevelsetvalue_0}$
(53,59) $\Wtransarg{\muxmulevelsetvalue_0}$
(54,41) $\Wtransarg{\muxmulevelsetvalue_0}$
(45,56) $\Wtransarg{\muxmulevelsetvalue_0}$
(29,69) $b_1$
(53,77) $b_2$

The region $\MInteresting = \MLeft \cup \MSingular \cup \MRight$ in
geometric coordinate space

[scale=.6, grid = false, tics=5, trim=-.5cm -1cm -1cm -.5cm, clip]MinterestingCartesiancoord.pdf
(40,10) $\Sigma_0$
(35,24) Cauchy stability region
(-2,20) $(x^2,x^3) \in \mathbb{T}^2$
(2,29) $t$
(9,26) $x^1 \in\mathbb{R}$
(15,50) $\Upsilon(\MLeft)$
(65,69) $\Upsilon(\twoargmumuxtorus{\timefunction_0}{-\muxmulevelsetvalue_0})$
(30,28) $\Upsilon(\inthyp{\timefunction_0}{[- \rightu,\leftu]})$
(45,99) $\Upsilon(\inthyp{0}{[- \rightu,\leftu]})$
(75,75) $\Upsilon(\MRight)$
(7,45) $\Upsilon(\twoargmumuxtorus{-\timefunction_0}{0})$
(39,78) $\Upsilon(\twoargmumuxtorus{0}{0})$
(59,89) $\Upsilon(\twoargmumuxtorus{0}{-\muxmulevelsetvalue_0})$
(18,28) $\Upsilon(\datahypfortimefunctiontwoarg{0}{[\timefunction_0,0]})$
(61,47) $\Upsilon(\datahypfortimefunctiontwoarg{-\muxmulevelsetvalue_0}{[\timefunction_0,0]})$
(38,54) $\Upsilon(\MSingular)$
(45,81) $\Upsilon(\mathcal{B}^{[0,\muxmulevelsetvalue_0]})$

The region $\Upsilon(\MInteresting)$ from Fig. <ref>,
i.e., $\MInteresting$ mapped into Cartesian coordinate space under $\Upsilon$

* Fix any of the “admissible” background (shock-forming) simple isentropic plane-symmetric
solutions that we construct in Appendix <ref>,
where we define “admissible” in Def. <ref>;
for each background solution, only a single Riemann invariant[By our definition
(<ref>) of Riemann invariants,
the fluid state $\RRiemannPS = 0$ corresponds to a fluid with a non-zero constant density
$\bar{\varrho} > 0$.]
is non-vanishing, and we denote it by $\RRiemannPS$.
For each such background solution, there exist numbers
$\timefunction_0 < 0$, $\muxmulevelsetvalue_0 > 0$,
$\rightu > 0$,
$\farrightu > \rightu$,
and $\leftu > 0$
such that the data of $\RRiemannPS$ is
compactly supported[We made the assumption that the initial data of the background Riemann invariant
are compactly supported only for convenience;
this assumption could be eliminated without much additional effort.]
in the Cauchy hypersurface portion
\Sigma_0^{[-\rightu,\leftu]}
\left\lbrace
(0,-x^1,x^2,x^3) \ | \
-x^1 \in [-\rightu,\leftu],
\,
(x^2,x^3) \in \mathbb{T}^2
\right\rbrace
in Cartesian coordinate space $\mathbb{R}_t \times \mathbb{R}_{x^1} \times \mathbb{T}^2$.

* Let $\Ntop \geq 24$ be an integer.

* Assume that the data-norm $\mathring{\Delta}_{\Sigma_0^{[-\farrightu,\leftu]}}^{\Ntop+1}$
defined in (<ref>),
which measures the $H^{\Ntop+1}(\Sigma_0^{[-\rightu,\leftu]})$-closeness of the perturbed
data to the background data, is sufficiently small. Note that the perturbed data do not have to be compactly supported.

Then the following conclusions hold, where the role of $\timefunction_0$ and $\muxmulevelsetvalue_0$ is described below.

Classical existence in the Cauchy stability region.
The perturbed solution exists classically with respect to the Cartesian coordinates
and remains uniformly bounded in the “Cauchy stability region”
trapped in between the flat hypersurface portion $\Sigma_0^{[-\rightu,\leftu]}$ and
the curved hypersurface portion $\inthyp{\timefunction_0}{[- \rightu,\leftu]}$,
where $\inthyp{\timefunction_0}{[- \rightu,\leftu]}$ is a portion of the $\timefunction_0$-level set of the
time function $\newtimefunction$, described below.
We will not further discuss the behavior of the
solution in the Cauchy stability region,
noting only that it is depicted in Fig. <ref>.

Acoustic geometry and regular behavior in the geometric coordinates.
There exists an eikonal function $u$ solving the eikonal equation
$(\gfour^{-1})^{\alpha \beta} \partial_{\alpha} u \partial_{\beta} u = 0$
such that the geometric coordinates $(t,u,x^2,x^3)$ form a global coordinate system on
the compact subset $\MInteresting = \MLeft \cup \MSingular \cup \MRight$
depicted in Fig. <ref>.
The figure is a subset of geometric coordinate space
$\mathbb{R}_t \times \mathbb{R}_u \times \mathbb{T}^2$,
\subset
\lbrace
\in
\mathbb{R} \times \mathbb{R} \times \mathbb{T}^2
\ | \ - \rightu \leq u \leq \leftu
\rbrace$.
The acoustical metric
$\gfour$ is the fluid-solution-dependent Lorentzian metric (see (<ref>))
that captures the intrinsic geometry of sound waves in the flow.
The level sets of $u$, which we denote by $\nullhyparg{u}$ and which are depicted as subsets of Cartesian coordinate space
in Fig. <ref>,
are characteristic surfaces for the system (<ref>)–(<ref>).
The fluid blowup is not visible in the geometric coordinates:
the fluid solution and its up-to-mid-order partial derivatives in the geometric coordinate system
remain bounded on $\MInteresting$.

The singularity is caused by the infinite density of the characteristics.
With respect to the Cartesian coordinates, there is singularity formation,
described below, that coincides with the vanishing of the inverse foliation density
$\upmu \eqdef - \frac{1}{(\gfour^{-1})^{\alpha \beta}\partial_{\alpha} t \partial_{\beta} u}$.
The vanishing of $\upmu$ signifies the infinite density of the level sets of $u$
(viewed as a function of the Cartesian coordinates). Within $\MInteresting$,
this singular behavior occurs
precisely along the singular boundary portion $\mathcal{B}^{[0,\muxmulevelsetvalue_0]}$
depicted in Fig. <ref>.

Behavior of the change of variables map.
The change of variables map $\Upsilon(t,u,x^2,x^3) = (t,x^1,x^2,x^3)$
is a homeomorphism from $\MInteresting$ onto its image in Cartesian coordinate space.
It Jacobian determinant is $\approx - \upmu$, and since $\upmu > 0$ on
$\MInteresting \backslash \mathcal{B}^{[0,\muxmulevelsetvalue_0]}$,
the map is a diffeomorphism on $\MInteresting \backslash \mathcal{B}^{[0,\muxmulevelsetvalue_0]}$.

The shock singularity in the Cartesian differential structure.
There exists a $\Sigma_t$-tangent vectorfield $X$
of Euclidean length $\sqrt{\sum_{a=1}^3 (X^a)^2} \approx 1$ such that the following occurs.
Within the compact subset
$\Upsilon(\MInteresting) = \Upsilon(\MLeft \cup \MSingular \cup \MRight)$
of Cartesian coordinate space
depicted in Fig. <ref>,
with respect to the Cartesian coordinates,
$|X^a \partial_a v^1|$
$|X^a \partial_a \varrho|$
blow up precisely on the singular boundary portion
$\Upsilon(\mathcal{B}^{[0,\muxmulevelsetvalue_0]})$, described below.
The solution is smooth in $\Upsilon\left(\MInteresting \backslash \mathcal{B}^{[0,\muxmulevelsetvalue_0]} \right)$,
and for $\alpha = 0,1,2,3$, $i=1,2,3$, and $A=2,3$,
the fluid quantities
$\Flatcurl v^i$,
$\partial_i \Ent$,
$\gfour_{ab} \Yvf{A}^a \partial_{\alpha} v^b$
remain bounded on all of $\Upsilon(\MInteresting)$,
where $\Yvf{2}$ and $\Yvf{3}$ are the vectorfields from (<ref>),
and they are tangent to $\nullhyparg{u} \cap \Sigma_t$.
Moreover, the up-to-mid-order derivatives of these quantities
with respect to the $\nullhyparg{u}$-tangent vectorfields
$\lbrace \Lunit, \Yvf{2}, \Yvf{3} \rbrace$
from (<ref>)
also remain bounded on all of $\Upsilon(\MInteresting)$.

The structure of $\mathcal{B}^{[0,\muxmulevelsetvalue_0]}$ and $\partial_- \mathcal{B}^{[0,\muxmulevelsetvalue_0]}$.
The singular boundary $\mathcal{B}^{[0,\muxmulevelsetvalue_0]}$,
viewed as a subset of geometric coordinate space
$\mathbb{R}_t \times \mathbb{R}_u \times \mathbb{T}^2$,
is contained in the top boundary of $\MInteresting$
and is a subset of the level set $\lbrace \upmu = 0 \rbrace$, specifically, a portion that
can be realized as the limit as $\mulevelsetvalue \downarrow 0$ of portions of the level sets
$\lbrace (t,u,x^2,x^3) \ | \ \upmu(t,u,x^2,x^3) = \mulevelsetvalue \rbrace$ that are
either null or spacelike with respect to the acoustical metric $\gfour$.
$\mathcal{B}^{[0,\muxmulevelsetvalue_0]}$ is
an embedded $3$-dimensional submanifold-with-boundary of geometric coordinate space
$\mathbb{R}_t \times \mathbb{R}_u \times \mathbb{T}^2$.
Moreover, we can decompose
$\mathcal{B}^{[0,\muxmulevelsetvalue_0]} = \bigcup_{\muxmulevelsetvalue \in [0,\muxmulevelsetvalue_0]} \twoargmumuxtorus{0}{-\muxmulevelsetvalue}$,
where each $\twoargmumuxtorus{0}{-\muxmulevelsetvalue}$
is a $C^{1,1}$, $2$-dimensional, spacelike submanifold of geometric coordinate space
$\mathbb{R}_t \times \mathbb{R}_u \times \mathbb{T}^2$
that is a graph over $\mathbb{T}^2$
such that $\upmu X \upmu \equiv - \muxmulevelsetvalue$ along $\twoargmumuxtorus{0}{-\muxmulevelsetvalue}$.
The torus $\twoargmumuxtorus{0}{0}$,
which we refer to as the crease and also denote by $\partial_- \mathcal{B}^{[0,\muxmulevelsetvalue_0]}$,
is the past boundary of $\mathcal{B}^{[0,\muxmulevelsetvalue_0]}$.

The structure of $\Upsilon(\mathcal{B}^{[0,\muxmulevelsetvalue_0]})$ and
$\Upsilon(\partial_- \mathcal{B}^{[0,\muxmulevelsetvalue_0]})$.
The change of variables map $\Upsilon(t,u,x^2,x^3) = (t,x^1,x^2,x^3)$
is a homeomorphism from the singular boundary portion
onto its image $\Upsilon(\mathcal{B}^{[0,\muxmulevelsetvalue_0]})$
in Cartesian coordinate space. Moreover, $\Upsilon$ is a diffeomorphism from
$\mathcal{B}^{[0,\muxmulevelsetvalue_0]} \backslash \partial_- \mathcal{B}^{[0,\muxmulevelsetvalue_0]}$
onto its image
$\Upsilon(\mathcal{B}^{[0,\muxmulevelsetvalue_0]} \backslash \partial_- \mathcal{B}^{[0,\muxmulevelsetvalue_0]})$.
The image set $\Upsilon(\mathcal{B}^{[0,\muxmulevelsetvalue_0]})$
is an embedded $3$-dimensional submanifold-with-boundary in Cartesian coordinate space
that has regularity[The $C^{1,1/2}$ embedding of $\Upsilon(\mathcal{B}^{[0,\muxmulevelsetvalue_0]})$ is
provided by Prop. <ref>; it is the map
$(z,x^2,x^3) \rightarrow \Upsilon \circ \embeddingofsingularboundaryintogeometriccoordinatespace(\sqrt{z},x^2,x^3)$
from the proposition.
$C^{1,1/2}$ with respect to the Cartesian coordinates.
In addition, $\Upsilon(\mathcal{B}^{[0,\muxmulevelsetvalue_0]} \backslash \partial_- \mathcal{B}^{[0,\muxmulevelsetvalue_0]})$
is a null hypersurface with respect to the acoustical metric $\gfour$ on $\Upsilon(\MInteresting)$,
$\Upsilon(\mathcal{B}^{[0,\muxmulevelsetvalue_0]} \backslash \partial_- \mathcal{B}^{[0,\muxmulevelsetvalue_0]})$ is ruled,
in a degenerate sense explained in Remark <ref>,
by integral curves of the $\gfour$-null vectorfield $\Lunit = \Lunit^{\alpha} \partial_{\alpha}$
on $\Upsilon(\MInteresting)$.
for $\muxmulevelsetvalue \in [0,\muxmulevelsetvalue_0]$,
the images $\Upsilon(\twoargmumuxtorus{0}{-\muxmulevelsetvalue})$ of the tori $\twoargmumuxtorus{0}{-\muxmulevelsetvalue}$,
including the crease $\partial_- \mathcal{B}^{[0,\muxmulevelsetvalue_0]} = \twoargmumuxtorus{0}{0}$,
are $2$-dimensional, $C^{1,1}$, $\gfour$-spacelike submanifolds of Cartesian coordinate space
$\mathbb{R}_t \times \mathbb{R}_{x^1} \times \mathbb{T}^2$
that are graphs over $\mathbb{T}^2$.

Rough time functions reveal the tori foliating $\mathcal{B}^{[0,\muxmulevelsetvalue_0]}$.
There exists a one-parameter family of rough time functions
$\lbrace \timefunctionarg{\muxmulevelsetvalue} \rbrace_{\muxmulevelsetvalue \in [0,\muxmulevelsetvalue_0]}$,
each with range $[\timefunction_0,0]$,
such that the level sets $\lbrace \timefunctionarg{\muxmulevelsetvalue} = \timefunction \rbrace$
with $\timefunction \in [\timefunction_0,0)$
do not intersect $\mathcal{B}^{[0,\muxmulevelsetvalue_0]}$,
while $\lbrace \timefunctionarg{\muxmulevelsetvalue} = 0 \rbrace \cap \mathcal{B}^{[0,\muxmulevelsetvalue_0]} = \twoargmumuxtorus{0}{-\muxmulevelsetvalue}$.
That is, each rough time function reveals the structure of the two-dimensional torus $\twoargmumuxtorus{0}{-\muxmulevelsetvalue}$
in $\mathcal{B}^{[0,\muxmulevelsetvalue_0]}$,
and $\twoargmumuxtorus{0}{-\muxmulevelsetvalue} \subset \lbrace \timefunctionarg{\muxmulevelsetvalue} = 0 \rbrace$.
Moreover, each $\timefunctionarg{\muxmulevelsetvalue}$ is one degree less differentiable with respect
to the geometric coordinates than the fluid solution,
and the tori $\twoargmumuxtorus{0}{-\muxmulevelsetvalue}$ are two degrees
less differentiable than the fluid solution.

Moreover, there exists a rough time function $\newtimefunction$ with range $[\timefunction_0,0]$ that foliates
$\newtimefunction$ is
$C^{1,1}$ with respect to the geometric coordinates - and not more regular;
see Remark <ref>.
Two of its level sets intersected with $\lbrace u \in [- \rightu,\leftu] \rbrace$,
namely $\inthyp{0}{[- \rightu,\leftu]}$ and
$\inthyp{\timefunction_0}{[- \rightu,\leftu]}$, are depicted in Fig. <ref>.
Finally, $\mathcal{B}^{[0,\muxmulevelsetvalue_0]} \subset \lbrace \newtimefunction = 0 \rbrace$.

In our main results,
we have chosen to assume the smallness of the norm $\mathring{\Delta}_{\Sigma_0^{[-\farrightu,\leftu]}}^{\Ntop+1}$
defined in (<ref>)
because this immediately allows us to conclude that
our results hold for open sets of initial data that are $H^{\Ntop+1}$-close
to the data of a background solution.
However, the main assumptions on the initial data that we use in our PDE analysis
are actually stated in Sect. <ref>,
in terms of data-size parameters that satisfy
assumptions stated in Sect. <ref>.
All these assumptions follow as consequences of
the smallness of the norm $\mathring{\Delta}_{\Sigma_0^{[-\farrightu,\leftu]}}^{\Ntop+1}$.
That is, the proof of our main results would go through under only
the assumptions stated in Sect. <ref>
and the parameter-size assumptions stated in Sect. <ref>,
modulo the remarks we make at the beginning of Sect. <ref>.

§.§ Remarks on the main results and methods

Before proceeding, we make a series of remarks about our main results
and our proof framework,
with an emphasis on the new methods we introduce in this paper,
how they connect to methods developed in other papers,
and how they connect to open problems.

* (Use of acoustic geometry to detect the singularity and to “hide it” by unfolding the characteristics).
As in Christodoulou's breakthrough work <cit.> on irrotational and isentropic shock formation,
our analysis fundamentally relies
on nonlinear geometric optics (which we also loosely refer to as “the acoustic geometry”),
implemented via an eikonal function $u$, which solves
the eikonal equation $(\gfour^{-1})^{\alpha \beta} \partial_{\alpha} u \partial_{\beta} u = 0$.
The acoustical metric $\gfour$ is a Lorentzian metric
whose Cartesian components $\gfour_{\alpha \beta}$
are functions of the fluid variables;
see (<ref>) for the precise formula.
The level sets $\nullhyparg{u}$ of $u$ are null hypersurfaces (also known as “characteristics” or “characteristic hypersurfaces”),
and they correspond to the propagation of sound waves.
Infinite density of the level sets of $u$ (viewed as a function of the Cartesian coordinates)
signifies the formation of a shock,
and in the regime under study, it coincides with the blowup of first-order partial derivatives of $v$ and $\varrho$
in directions transversal to the $\nullhyparg{u}$.
As is shown in Fig. <ref>,
these phenomena occur along the entire singular boundary.
Moreover, with the help of $u$, we can construct a “geometric coordinate system”
$(t,u,x^2,x^3)$, relative to which the solution remains rather smooth.
This is crucial for our derivation of PDE estimates up to top-order.

* (The structure of the singular boundary and the crease).
We prove that in the solution regime under study,
relative to a differential structure on spacetime tied to the eikonal function,
the singular boundary $\mathcal{B}$ has the structure of a $3D$ submanifold-with-boundary;
see Fig. <ref>.
Its past boundary $\partial_- \mathcal{B}$ is the crease,
which we prove is a $2D$ acoustically spacelike submanifold,
where our notion of “spacelike” is
relative to $\gfour$.
These structures are stable, and their availability is fundamental for our approach.
These structures also have important implications for the shock development problem,
described below.

* (Prior works).
There are many prior works on shock formation for compressible Euler solutions,
including Riemann's famous work <cit.> in one spatial dimension.
Although much less is known in multi-dimensions, there has been dramatic progress over the last several decades,
starting with Alinhac's foundational works on quasilinear wave equations
sA1999a,sA1999b,sA2002 and
Christodoulou's breakthrough monograph on the irrotational and isentropic
relativistic Euler equations <cit.>;
see Sect. <ref> for a discussion of some additional key developments in the history of the subject.
While prior results have led to a revolution in our understanding of multi-dimensional
shock formation, they all were limited in one or more of the following ways:
i) they treated only irrotational and isentropic solutions;
ii) the methods allowed one only to follow the solution to the constant-Cartesian-time
hypersurface of first blowup;
iii) the methods applied only to fully non-degenerate singularities (see below for their definition)
which, as it turns out, corresponds to understanding the blowup at the unique first (relative to Cartesian time)
point that is contained inside a strictly convex crease;
iv) the methods yielded a description only of some implicit portion of the
boundary of the maximal development of the data
and in particular yielded only the portion of the crease
that is tangent to some flat spacelike hypersurface $\Sigma$
such that a neighborhood of the crease lies in its future.
A key point is that without strict convexity,
there can be points in the crease such that
no open neighborhood of them is accessible
through such an approach.
In particular, for the solutions
we treat in our main results, the full structure of the crease
is not accessible through any of these approaches;
see Fig. <ref>.

* (Relevance for the shock development problem).
The crease and its sub-manifold structure are crucial ingredients
needed to set up the shock development problem,
which, as we mentioned earlier,
is the problem of describing the transition of the solution
from classical to weak, past the initial singularity.
It turns out that the crease is part of the maximal classical development
and, once one constructs the weak solution, it will also be part of the weak solution;
see  <cit.>*Section 1.5 for a detailed discussion of the connection between
the maximal classical development and the weak solution.
Away from symmetry, the shock development problem is an outstanding open problem.
In Christodoulou's approach to studying shock developments
under the simplifying assumption that the vorticity and entropy are vanishing[The irrotational and isentropic weak solutions
constructed by Christodoulou in <cit.> are not solutions to the compressible Euler equations because they do not
respect the jump in entropy and vorticity that must occur across the shock hypersurface. Instead, they are solutions
to a closely related hyperbolic PDE system that is equivalent to the compressible Euler equations for classical
isentropic and irrotational classical solutions. ]
across the shock hypersurface <cit.>,
the rest of the singular boundary is also important because it
plays the role of a mathematical barrier in the construction
of the weak solution. However, aside from the crease, the singular boundary
is not part of the weak solution. More precisely, in general,
there is a portion of the maximal classical development that does not
agree with the weak solution because the classical development does not account for the shock hypersurface.
The shock hypersurface, which we denote by $\mathcal{K}$ in Fig. <ref>,
is not part of the constructions of this paper, and we mention it to highlight that the
weak solution – once it is constructed – will disagree with the classical solution in the region
in between $\mathcal{B}$ and $\mathcal{K}$.
We also highlight that $\mathcal{K}$ is supersonic relative to the acoustical metric in the
maximal classical development,
and that $\mathcal{K}$ is “allowed” to emerge from the crease only if one imposes a weak
formulation of the flow, i.e., there is no such thing as “$\mathcal{K}$” in the classical formulation.

* (The Cauchy horizon).
In our forthcoming paper, we will describe the emergence of a Cauchy horizon
from the crease. The Cauchy horizon, denoted by $\underline{\mathcal{C}}$ in
Fig. <ref>,
is a $\gfour$-null hypersurface with past boundary equal to the crease,
and it is also a crucial ingredient in the setting up the shock development problem.
Like the crease, the Cauchy horizon is part of the maximal classical development
and the weak solution (once one constructs it). However, in that paper, we will show that
the solution remains smooth up to the Cauchy horizon (away from the crease, that is),
which is in stark contrast to what happens along the singular boundary.

* (Solution regimes other than perturbations of simple isentropic plane-waves).
Simple isentropic plane-waves, mentioned already in
Theorem <ref>, are compressible Euler solutions such that
$\RRiemann = \RRiemann(t,x^1)$
$\LRiemann = v^2 = v^3 \equiv 0$,
where the “almost[For isentropic plane-symmetric solutions,
$\almostRiemann_{(\pm)}$ are Riemann invariants.]
Riemann invariants” $\almostRiemann_{(\pm)}$
are defined in Def. <ref>.
Although for definiteness we have focused on general small perturbations of such solutions,
the techniques we develop here are robust and can be applied to other solution regimes of physical interest,
such as the regime corresponding to small, spatially-decaying perturbations[In this regime, the dispersive tendency
of sound waves competes against the transport phenomena associated to vorticity and entropy,
possibly leading to exceptionally complicated long-time dynamics.]
of non-vacuum constant fluid states in $\mathbb{R}^{1+3}$,
where dispersive effects play a fundamental role in the dynamics.

* (Building blocks that can accommodate degeneracies).
We have focused our attention on perturbations of
simple isentropic plane-symmetric waves because such solutions are of physical interest
and, from the point of view of analysis, they are challenging to study because
they exhibit degeneracies tied to lack of strict convexity (as we mentioned above)
of their singular boundaries. This is the first paper to fully grapple with
these degeneracies and as such, our results are new even in the sub-class of irrotational
and isentropic solutions.
The presence of these degeneracies forced us
to develop techniques that we anticipate can be used,
in view of finite speed of propagation,
as building blocks to study
the global structure of much more general shock-forming solutions across space.

* (New, solution-dependent foliations are needed).
To access the entire crease/singular boundary, in general, one cannot exclusively rely on arguments based
on analysis on regions bounded by the characteristic surfaces $\nullhyparg{u}$,
surfaces $\Sigma_t$ of constant Cartesian time,
or, for that matter, any other family of surfaces that are “pre-specified” in the sense that they are
explicitly parameterized with respect to the Cartesian coordinates
(e.g., one can not generally use “tilted” spacelike hypersurfaces that are planes with respect to the Cartesian coordinates).
In particular, the crease is typically not contained in a fixed $\nullhyparg{u}$
or $\Sigma_t$; see Fig. <ref>.
Hence, one of our key new ideas in the paper is to replace the surfaces $\Sigma_t$ with better ones.
That is, we construct foliations by spacelike hypersurfaces
that are precisely adapted to the shape of the singular boundary
(though we also fundamentally rely on foliations by the $\nullhyparg{u}$).
We construct these foliations by “flowing out” (see below), along
the integral curves of a “generating vectorfield” that is transversal to $\nullhyparg{u}$,
from co-dimension $2$ topological tori that are adapted to the anticipated shape of the singularity.
In Fig. <ref>, the generating vectorfield is denoted by $\Wtransarg{0}$, and
crucially, it is also transversal to the hypersurface $\datahypfortimefunctionarg{0}$,
which is foliated by the tori, two of which are denoted by
in Fig. <ref>,
the generating vectorfield is denoted by $\Wtransarg{\muxmulevelsetvalue}$, and
it is also transversal to the hypersurface $\datahypfortimefunctionarg{- \muxmulevelsetvalue}$,
which is foliated by the tori, two of which are denoted by
While there have been many
other works on the formation of shocks without symmetry assumptions, described in
Sect. <ref>,
those works have relied on a “background” time function (typically the Cartesian time function $t$)
and the corresponding foliations.
Note that the PDE analysis of compressible Euler flow,
which inevitably involves some kind of energy estimates with respect to the spacelike foliations,
necessarily halts when the spacelike foliations first intersect the
singular boundary, because singularities form at the intersection points.
This limits the portion of the singular
boundary that can be detected through background foliations.

* (Remarks on the simpler fully non-degenerate sub-regime).
Despite the previous point, we note that there is one sub-regime in which background time functions
can be used to detect the structure of the crease
(though perhaps not a full neighborhood of it in the singular boundary):
the regime in which the crease and singular boundary are acoustically strictly convex
(see below and Sect. <ref>).
In Fig. <ref>, we depict an acoustically strictly convex singular boundary.

(50,10) $\Sigma_0$
(67,24) $(x^2,x^3) \in \mathbb{T}^2$
(92,38) $t$
(70,34) $u \in\mathbb{R}$
(48,35) $\partial_- \mathcal{B}$
(60,76) $\mathcal{B}$
(52,44) $b_*$

An acoustically strictly convex crease and singular boundary in geometric coordinate space

[scale=.36, grid = false, tics=5, trim= 0cm 0cm 0cm -1.5cm, clip]crease_convex_cartesian.pdf
(50,10) $\Sigma_0$
(1,22) $(x^2,x^3) \in \mathbb{T}^2$
(.5,35) $t$
(14,28) $x^1 \in \mathbb{R}$
(39,24) $\partial_- \mathcal{B}$
(54,65) $\mathcal{B}$
(40.5,33) $b_*$

The same crease and singular boundary in Cartesian coordinate space

Acoustically strictly convex crease and singular boundary

We emphasize that symmetric solutions and their general perturbations fall outside of
the acoustically strictly convex regime.
Away from symmetry, the acoustically strictly convex regime was first studied by Alinhac
in the context of shock-forming solutions
to quasilinear wave equations, where he followed the solution to the constant-Cartesian-time hypersurface
of first blowup.
In particular, due to the strict convexity of the crease, the singular set
within the constant-Cartesian-time hypersurface of first blowup is an isolated point,
denoted by $b_*$ in Fig. <ref>,
and the isolated nature of this point was fundamental for his approach.
We refer to such singularities as fully non-degenerate;
see Sect. <ref> for further discussion.

* (Acoustical transversal convexity)
Here and throughout, by an “acoustically strictly convex”
hypersurface $\mathcal{B}$, we mean that relative to a system of geometric
coordinates $(t,u,x^2,x^3)$ tied to an acoustic eikonal function $u$,
the tangent planes to $\mathcal{B}$ lie below $\mathcal{B}$;
see Fig. <ref> for an example of such a surface.
The key point is that acoustical strict convexity
is absent
in the singular boundaries of open sets of physically relevant shock-forming solutions,
including perturbations of isentropic plane-symmetric simple waves.
To prove our main results, we rely on an acoustical transversal convexity
(we will refer to it as “transversal convexity” for short),
which is substantially weaker than strict convexity
and it is close to optimal in the following sense:
without transversal convexity,
the crease can fail to have the structure of a $2D$ sub-manifold,
which would drastically alter the qualitative character of the singular boundary.
We highlight that transversal convexity can be ensured along a portion of the singular boundary
containing the crease by making appropriate
open assumption on the initial data. In the present paper, we
ensure transversal convexity via the data assumptions stated in (<ref>).
Roughly, by an “acoustically transversally convex”
hypersurface $\mathcal{B}$,
we mean that relative to the geometric
coordinates $(t,u,x^2,x^3)$, there exists
a family of lines that are tangent to $\mathcal{B}$  – as opposed to the full tangent plane –
such that
i) the lines lie below[Our approach to studying the singular boundary is local,
and thus, in our main results, we only need the tangent lines to locally lie below it.]
and ii) the lines are transversal to the level sets of $u$.
In the context of Fig. <ref>,
lack of strict convexity is exhibited by the observation that at $b_1$,
the straight line that is tangent to $\twoargmumuxtorus{0}{0}$ (this line is not drawn in the figure)
lies above $\mathcal{B}^{[0,\muxmulevelsetvalue_0]}$,
and similarly, at $b_2$,
the straight line that is tangent to $\twoargmumuxtorus{0}{-\muxmulevelsetvalue_0}$
lies above $\mathcal{B}^{[0,\muxmulevelsetvalue_0]}$.
In Fig. <ref>,
transversal convexity is exhibited, for example,
by the fact that along the curve in $\mathcal{B}^{[0,\muxmulevelsetvalue_0]}$ joining $b_1$ and $b_2$,
the tangent lines (which are not drawn) lie below $\mathcal{B}^{[0,\muxmulevelsetvalue_0]}$.

* (Rough time functions).
Our new foliations are level sets of a $\muxmulevelsetvalue$-indexed family of
“rough time functions,” denoted by $\timefunctionarg{\muxmulevelsetvalue}$,
where $\muxmulevelsetvalue$ is a non-negative real parameter.
The $\timefunctionarg{\muxmulevelsetvalue}$
solve a $\muxmulevelsetvalue$-dependent transport equation (see (<ref>))
with coefficients
and initial data (see (<ref>))
determined by the acoustic geometry
and the up-to-first-order derivatives of the solution.
Moreover, we choose the initial
data surface for $\timefunctionarg{\muxmulevelsetvalue}$
to be a level set of a first derivative of the acoustic geometry,
more precisely, the $-\muxmulevelsetvalue$ level set of a first derivative of the null lapse in a direction transversal
to the characteristics. More precisely, each level set of
$\timefunctionarg{\muxmulevelsetvalue}$ is the solution a transport equation
with (constant) initial data given on a topological torus that itself
is equal to the intersection of a level set of the null lapse with a level set
of one of its transversal derivatives.
In Fig. <ref>, we denote the initial data hypersurfaces for the transport equations
by $\datahypfortimefunctionarg{0}$
and $\datahypfortimefunctionarg{-\muxmulevelsetvalue}$,
we denote some tori in the data hypersurfaces by
and we denote the vectorfields that transport $\timefunctionarg{0}$ and $\timefunctionarg{\muxmulevelsetvalue}$
respectively by $\Wtransarg{0}$ and $\Wtransarg{\muxmulevelsetvalue}$.
These vectorfields are tangent to the level sets of
$\timefunctionarg{0}$ and $\timefunctionarg{\muxmulevelsetvalue}$ respectively,
and they are also tangent to the singular boundary along
the tori $\twoargmumuxtorus{0}{0}$ and $\twoargmumuxtorus{0}{-\muxmulevelsetvalue}$
We emphasize that this construction is quite delicate and that
it is crucial that we consider only the case $\muxmulevelsetvalue \geq 0$ in order to ensure that
$\timefunctionarg{\muxmulevelsetvalue}$ has a $\gfour$-timelike normal or $\gfour$-null normal.
In total, the construction is “dynamic” and “solution-adapted,”
and we achieve it through a bootstrap argument.
The singular boundary is provably one degree less differentiable than the fluid solution,
and the same is true for our rough time functions. Moreover, the tori (such as $\twoargmumuxtorus{-\timefunction}{-\muxmulevelsetvalue}$)
are two degrees less differentiable.
To handle this fundamental difficulty,
we work with multiple coordinate systems, as we describe in the next point.

* (Three distinct coordinate systems).
There is no known approach to studying multi-dimensional shock formation
that relies only on commuting the equations
with the Cartesian coordinate partial derivative vectorfields.
Because of the roughness of the time functions $\timefunctionarg{\muxmulevelsetvalue}$,
we cannot close our top-order energy estimates using only commutation vectorfields
that are adapted[However, for some crucial low-derivative-level estimates,
we do use commutators that are adapted to the level sets of $\timefunctionarg{\muxmulevelsetvalue}$.]
to their level sets;
such an approach would lead to the loss of a derivative.
Hence, as in other works on shocks, we construct appropriate commutators with the help of
the eikonal function $u$. In total, our study of the flow relies on understanding the behavior of the
fluid in three coordinate systems as well as carefully controlling the relationships
– which degenerate near the singular boundary –
between the coordinate systems:

* The standard Cartesian coordinates $(t,x^1,x^2,x^3)$, relative to which the
compressible Euler equations are initially formulated
and relative to which the singularity is visible.

* The geometric coordinates $(t,u,x^2,x^3)$ (where $u$ is the eikonal function),
for constructing suitable commutator
and multiplier vectorfields that are adapted to the singularity and have sufficient regularity,
and relative to which the solution remains rather smooth;
see Sect. (<ref>).

* The rough adapted coordinates $(\timefunctionarg{\muxmulevelsetvalue},u,x^2,x^3)$, where
the level sets of $\timefunctionarg{\muxmulevelsetvalue}$
are precisely adapted to the shape of the singular boundary.

* (Elliptic-hyperbolic estimates for the vorticity and entropy).
It is difficult to control the top-order derivatives of the vorticity and entropy
on regions that are adapted to the shape of the singular boundary.
In particular, the analytic framework we use, which is based on the
formulation of compressible Euler flow provided by Theorem <ref>,
requires that we prove that the vorticity and entropy gradient are exactly as differentiable
as the velocity and density, all the way up to the singular boundary.
To this end, we develop an upgraded version of the “elliptic-hyperbolic” integral identities
from <cit.>. The upgraded version, which we implement with the help of a new characteristic current,
allows us to handle severe degeneracies that arise near the crease.
The degeneracies are tied to the fact that the level sets of our time functions
$\timefunctionarg{\muxmulevelsetvalue}$ become asymptotically null the near crease
for $\muxmulevelsetvalue$ small and positive; see Figs. <ref>
We refer to Sect. <ref> for our derivation of the elliptic-hyperbolic identities.

* (Transversal convexity).
As we have mentioned, our analysis relies on an open data-assumption that we call
acoustical transversal convexity
(which we refer to as transversal convexity for short);
we refer to (<ref>) for the precise, technical
data estimates that capture transversal convexity.
In studying the flow,
we are able to propagate the transversal convexity all the way up to the singular boundary,
which in particular ensures that the level sets $\lbrace \upmu = 0 \rbrace$
and $\lbrace \muX \upmu = 0 \rbrace$
on RHS (<ref>) (which defines the crease) intersect transversally,
Our assumption of transversal convexity is weaker than the assumption that the graph
of $\lbrace \upmu = 0 \rbrace$, viewed as a subset of geometric coordinate space,
is strictly convex, and our weaker assumption is what allows us, for example,
to handle perturbations of symmetric solutions.
Equivalently, under our approach, transversal convexity allows
us to understand the structure of the crease and the singular boundary for singularities
that are more general than the fully non-degenerate ones featured in Fig. <ref>,
including the perturbations of simple isentropic plane-waves that we treat in detail in our main theorems,
as is depicted in Fig. <ref>.

Our transversal convexity assumption is close to optimal in the sense that without it,
the qualitative character of the singular boundary can dramatically change,
e.g., without transversal convexity,
the crease could fail to have the structure of a $2D$ submanifold,
and then one could not even properly set up the shock development problem
using the known approach <cit.>; see the next point.
The qualitative change of the structure of the crease in the absence of transversal convexity
can readily be seen in the context of the simple isentropic plane-symmetric solutions
that we study in Appendix <ref>. More precisely,
using the explicit formulas provided by Cor. <ref>,
one can construct simple isentropic plane-symmetric fluid initial data
such that transversal convexity fails and such that
the crease $\lbrace \upmu = 0 \rbrace \cap \lbrace \muX \upmu = 0 \rbrace$,
viewed as subset of $1+3$-dimensional geometric coordinate space,
fails to have the structure of a $2D$ submanifold. For example,
one can any consider smooth initial data such that the data-term
\frac{d}{du}
\antiderivativePSLmusourcetermfunction[\dataRRiemannPS(u)]
on RHS (<ref>)
achieves its negative minimum along a closed interval $I$ of $u$-values.
It is then easy to see,
using (<ref>),
and the fact that $\muX = \geop{u}$ for simple isentropic plane-symmetric solutions,
that $\muX \upmu = 0$ at all points in $\Sigma_{T_{\textnormal{Shock}}}$
where $\upmu$ vanishes, where $T_{\textnormal{Shock}}$ is the Cartesian time of first blowup.
In particular, the level sets $\lbrace \upmu = 0 \rbrace$
and $\lbrace \muX \upmu = 0 \rbrace$ coincide on the interval of $u$ values $I$ within
$\Sigma_{T_{\textnormal{Shock}}}$, signifying a dramatic failure of
transversal convexity. Moreover, using the characterization (<ref>)
of the crease, we see that the crease,
viewed as a subset of $1+3$-dimensional geometric coordinate space,
is the following set:
\left\lbrace (t,u,x^2,x^3) \ | \ t = T_{\textnormal{Shock}}, u \in I, (x^2,x^3) \in \mathbb{T}^2 \right\rbrace$,
which is a $3D$ submanifold-with-boundary (in particular, it is not a $2D$ submanifold).

* (Homeomorphism property of the change of variables map).
Finally, we highlight that  – thanks to the transversal convexity – we are able to prove
that the change of variables map
$\Upsilon(t,u,x^2,x^3) \eqdef (t,x^1,x^2,x^3)$
from geometric to Cartesian coordinates
is a homeomorphism all the way up to the singular boundary,
and it is a diffeomorphism away from the singular boundary;
see Prop. <ref>.
In particular, $\Upsilon$ is a bijection
on the region $\MInteresting$ depicted in Fig. <ref>,
and the level sets of $u$ never actually intersect
$\MInteresting$ (even along the singular boundary portion $\mathcal{B}^{[0,\muxmulevelsetvalue_0]}$!),
despite the fact that the density of the level sets of $u$ in Cartesian coordinate space becomes infinite along the singular boundary.
Under transversal convexity,
the map $u \rightarrow - u^3$ provides a crude scalar caricature of
the degenerate behavior of the $x^1$ Cartesian coordinate
along a constant Cartesian time slice $\Sigma_t$ that happens to be tangent to the crease.
This not just a mathematical curiosity; the bijective property of $\Upsilon$ up to and including the singular boundary
is crucial for formulating the shock development problem,
and this property can dramatically fail without transversal convexity.
As in the previous point, this difficulty
can readily be seen in the context of the simple isentropic plane-symmetric solutions
that we study in Appendix <ref>. For example, in the example discussed in the previous point,
using that $\geop{u} x^1 = - \upmu$ in plane symmetry, one finds that
$x^1 = x^1(t,u,x^2,x^3)$ is constant along the entire crease
\left\lbrace (t,u,x^2,x^3) \ | \ t = T_{\textnormal{Shock}}, u \in I, (x^2,x^3) \in \mathbb{T}^2 \right\rbrace
thereby exhibiting a dramatic failure of the injectivity of $\Upsilon$.
In the Cartesian coordinate space picture in $1+1$ dimensions (i.e., ignoring the $x^2$ and $x^3$ directions),
the shock formation in this example corresponds to a continuum of characteristic curves
all intersecting in a single point in Cartesian coordinate space.
This is in stark contrast to the behavior of the characteristic hypersurfaces
on the singular boundary in the transversally convex regime,
as we show in Prop. <ref>
and Fig. <ref>.

§.§ Works building up toward Theorem <ref>

In this section, we describe some key prior works that developed some of the technology
we use here. While we have already mentioned some of them,
here we provide additional details.
We refer to Sect. <ref> for a more extensive (though far from comprehensive)
discussion of the history of works on shock formation.

§.§.§ Christodoulou's framework from <cit.>

Our proof of Theorem <ref> relies on the framework of nonlinear geometric optics
developed by Christodoulou <cit.> in his groundbreaking proof of stable shock formation
for the irrotational and isentropic relativistic Euler equations in three spatial dimensions.
In this setting, the equations reduce to a scalar quasilinear wave equation for a potential function

\begin{align} \label{E:QLW}
(\gfour^{-1})^{\alpha \beta}(\pmb{\partial} \Phi) \partial_{\alpha} \partial_{\beta} \Phi
& = 0.
\end{align}

In (<ref>), $\gfour_{\alpha \beta}$ is the acoustical metric,
a Lorentzian metric whose rectangular component functions $\gfour_{\alpha \beta}$
are nonlinear functions of $\pmb{\partial} \Phi \eqdef (\partial_t \Phi, \partial_1 \Phi, \partial_2 \Phi, \partial_3 \Phi)$
depending on the equation of state:
$\gfour_{\alpha \beta} = \gfour_{\alpha \beta}(\pmb{\partial} \Phi)$.
To implement nonlinear geometric optics, Christodoulou constructed an
eikonal function $u$, that is a solution to the eikonal equation

\begin{align} \label{E:INTROEIKONAL}
(\gfour^{-1})^{\alpha \beta}(\pmb{\partial} \Phi) \partial_{\alpha} u \partial_{\beta} u
& = 0.
\end{align}

The level sets of $u$, denoted by $\nullhyparg{u}$, are characteristic for equation (<ref>).
Notice that the principal coefficients on LHS (<ref>) depend on $\pmb{\partial} \Phi$
and thus the evolution of $u$ is coupled to that of $\Phi$, signifying dynamic nature of the geometry
and the quasilinear nature of the flow.
In <cit.>, Christodoulou used $u$ to construct a system of geometric coordinates
$(t,u,\vartheta^1,\vartheta^2)$ on spacetime such that the solution remains rather smooth relative to these
coordinates. The formation of the singularity can be recovered as a degeneracy between the
geometric coordinates and the standard ones, $(t,x^1,x^2,x^3)$.
The degeneracy is signified by the vanishing of the inverse foliation density $\upmu$, defined by:

\begin{align} \label{E:INTROINVERSEFOLIATION}
\upmu
\eqdef - \frac{1}{(\gfour^{-1})^{\alpha \beta}(\pmb{\partial} \Phi) \partial_{\alpha} t \partial_{\beta} u}.
\end{align}

In the region of classical existence, one has $\upmu > 0$,
and when $\upmu \to 0$, the density of the $\nullhyparg{u}$ in Cartesian coordinate space
becomes infinite, signifying the “piling up” of the characteristics.
In Figs. <ref> and <ref>,
we show three distinct characteristic hypersurfaces piling up along the singular boundary,
along which $\upmu$ vanishes.
The hard part of the proof is to derive energy estimates in regions where $\upmu$ is small.
We refer readers to Sect. <ref> for further discussion on Christodoulou's framework.

§.§.§ Wave equations beyond fluid mechanics

The wave equations studied by Christodoulou in <cit.>
were Euler–Lagrange equations that were invariant with respect to the Poincaré group.
In <cit.>, we extended the results of <cit.>
to apply to all wave equations of type (<ref>)
and of type:

\begin{align} \label{E:INTROCOVARIANTWAVE}
\square_{\gfour(\Psi)} \Psi
& = \nullform(\pmb{\partial} \Psi, \pmb{\partial} \Psi)
\end{align}

that fail to satisfy Klainerman's null condition <cit.>,
where $\square_{\gfour(\Psi)}$ is the covariant wave operator of
$\gfour(\Psi)$ (see Def. <ref>)
and $\nullform$ is a null form relative to $\gfour$ (see Def. <ref>).

§.§.§ Nearly simple plane-symmetric waves

In <cit.>, the second author and his collaborators extended the methods of <cit.> and <cit.>
to prove stable shock formation
for a large class of quasilinear wave equations on the spacetime $\mathbb{R} \times \Sigma$,
where $\Sigma = \mathbb{R} \times \mathbb{T}$ was the two-dimensional spatial manifold.
The initial data we treated were analogs of the data from Theorem <ref>.
More precisely, the data were (asymmetric) perturbations of simple plane-symmetric waves,
which are solutions that depend only on $(t,x^1) \in \mathbb{R} \times \mathbb{R}$
and which feature a wave moving only in one direction[In plane symmetry, one can study
the flow by constructing Riemann invariants, in which case simple plane-symmetric waves would feature only a single non-zero Riemann invariant. See Appendix <ref>, in which we use Riemann invariants to construct the background
solutions whose perturbations we study in our main results.]
(say to the right).
As in the relativistic case, for irrotational and isentropic solutions to the $2D$ compressible Euler equations,
the dynamics reduces to a quasilinear wave equation of type (<ref>).
Hence, as a special case, the results of
yielded stable shock formation
for nearly simple and isentropic plane-symmetric solutions
to the $2D$ irrotational and isentropic compressible Euler equations.

§.§.§ A new formulation of the flow and stable shock formation in the presence of vorticity

In <cit.>, the second author and Luk
developed a new formulation of barotropic[Barotropic equations of state are such that the pressure
$p$ can be expressed as a function of the density $\varrho$, i.e., with no dependence on $\Ent$.]
compressible Euler flow with remarkable
regularity and null structures, which in many regimes, allows one to study the flow as if it was
a perturbation of the quasilinear wave equation (<ref>).
In <cit.>, the second author
derived a similar new formulation for all equations of state
in which the pressure is a function of the density and entropy,
thus allowing one to incorporate thermodynamic effects into the framework of <cit.>.
In particular, the equations of <cit.> include a system of
transport-div-curl equations for the vorticity and entropy, which allows one to propagate a gain of one
derivative for these quantities relative to standard estimates.
In this article, we use the equations of <cit.> to prove our main results,
and the gain in regularity is crucial for our approach.
In Theorem <ref>, we recall the new formulation of the flow derived in <cit.>.
In <cit.>, we derived a similar new formulation for the relativistic Euler equations.

In <cit.>, the second author and Luk used the equations of
<cit.> and the technology of dC2007,jS2016b
to prove the first stable shock formation result for the $2D$ compressible Euler equations.
The authors treated open sets of initial data with vorticity that are close to the
data of an irrotational simple plane-wave solution.
The main theorem yielded the full structure of the set of blowup-points
within the constant-time hypersurface $\Sigma_{T_{\textnormal{Shock}}}$ of first blowup;
in the context of Fig. <ref>,
the authors understood the structure of
$\Sigma_{T_{\textnormal{Shock}}} \cap \mathcal{B}$.

In <cit.>, the second author and Luk used the equations of <cit.>
to extend the results of <cit.>
to the $3D$ case in the presence of vorticity and entropy.
The proof was much more difficult than the $2D$ case because
our regularity theory for the vorticity and entropy relied
on elliptic estimates across space,
which are difficult to derive near the shock.
In particular, the authors used the elliptic estimates to handle the vorticity stretching term in the equations,
which, as is well-known, vanishes for $2D$ barotropic solutions.
These elliptic estimates relied on the
transport-div-curl equations for the vorticity and entropy that we derived in <cit.>,
and our approach yielded elliptic estimates only along complete, flat hypersurfaces of constant time.
As in the $2D$ case, this approach yielded the full structure of the set of blowup-points
within the constant-time hypersurface $\Sigma_{T_{\textnormal{Shock}}}$ of first blowup,
in the context of Fig. <ref>,
the singular set portion $\Sigma_{T_{\textnormal{Shock}}} \cap \mathcal{B}$.
As we have already mentioned, in order to study the flow beyond the
Cartesian-flat hypersurface $\Sigma_{T_{\textnormal{Shock}}}$,
one needs additional ingredients,
including the identities that we describe
in Sect. <ref>.

§.§.§ Remarkable localized integral identities

In <cit.>, we derived new localized integral identities for $3D$ compressible Euler flow
that allow us to extend the elliptic estimates for the vorticity and entropy to yield
estimates on arbitrary spacetime regions that are bounded by
spacelike or null hypersurfaces.
We use this crucial ingredient in the present paper to derive
top-order vorticity and entropy estimates up to the singular boundary.

§.§ Additional history and results tied to shocks and singularities

A fundamental issue, discovered by Riemann <cit.> in the context of one spatial dimension, is that initially smooth
compressible Euler solutions typically develop shock singularities in finite time.
Roughly, shocks are singularities such that $\rho$, $v$, $\Ent$ remain bounded
but some first derivative of $\rho$ and $v$ blows up. Riemann's analysis relied on his discovery of Riemann invariants for isentropic solutions in $1D$. Relative to Riemann invariants, the equations reduce to a quasilinear system for two transport equations with distinct characteristic directions; see Appendix <ref>. The proof of the blowup of the solution's first derivatives
then follows from differentiating the equations to obtain a Riccati-type structure, much like in the model case of Burgers' equation $\partial_t \Psi + \Psi \partial_x \Psi = 0$.

§.§.§ The advanced state of the $1D$ theory and the key difference with multi-dimensions

In $1D$, for large class of hyperbolic quasilinear systems (such as strictly hyperbolic systems),
there is an advanced theory capable of describing the global behavior of solutions,
including the formation of shocks and the subsequent interactions of the shock waves.
We refer readers to the compendium <cit.> for a history of the subject
and a comprehensive introduction to the main techniques.

A principal reason for the advanced state of the $1D$ theory is that the
equations are well-posed for initial data in appropriate bounded variation (BV) spaces.
The state of affairs is dramatically different in multi-dimensions;
Rauch's fundamental work <cit.> showed that in multi-dimensions,
quasilinear hyperbolic systems are typically ill-posed for data in BV spaces.
In fact, the only known well-posedness results in multi-dimensions are for initial
data in $L^2$-type Sobolev spaces. For this reason, in multi-dimensions,
one is forced to derive energy estimates, which can be incredibly difficult in regions
containing singularities; this is the main reason why the theory
of multi-dimensional shock waves is so much less developed compared to the $1D$ case.

§.§.§ The first blowup-result in multi-dimensions without symmetry: proof by contradiction

Providing a constructive proof of shock formation in higher dimensions turns out to be a very hard problem.
In a nutshell, the reason is that away from $1D$, it seems necessary to carefully track the evolution
of characteristic hypersurfaces, which is much more difficult compared to the $1D$ case.
The characteristic geometry (e.g., the eikonal function $u$ in the context of the present article)
is much more difficult to construct and control, and, crucially, all proofs of even local well-posedness
rely on energy estimates in $L^2$-based Sobolev spaces, which are difficult to derive near singularities.
In <cit.>, Sideris proved an influential, non-constructive
stable blowup-result in $3D$, the first one for multi-dimensional
compressible Euler flow. He assumed that the equation of state is
barotropic[Barotropic equations of state are such that $p = p(\varrho)$.]
and that it satisfies a convexity assumption.
His proof applied to a large set of data, but it did not reveal the nature of the blowup;
his arguments relied on virial-type identities,
and he showed blowup through a contradiction-argument.

§.§.§ Alinhac's constructive proof of the formation of isolated singularities

The first constructive results on shock formation in multidimensions were by Alinhac
sA1999a,sA1999b,sA2002, who used nonlinear geometric optics (i.e., he constructed characteristic surfaces)
and Nash–Moser estimates to prove stable shock formation in $2D$
for a class of quasilinear wave equations
of the form:

\begin{align} \label{E:ALINHACSWAVEEQUATIONS}
(\gfour^{-1})^{\alpha \beta}(\pmb{\partial} \Phi) \partial_{\alpha} \partial_{\beta} \Phi
& = 0,
\end{align}

whenever they fail to satisfy Klainerman's null condition <cit.>.
He used the Nash–Moser estimates to avoid derivative loss in his control of the regularity of the characteristic surfaces.
His proof applied to open sets of non-degenerate initial data,
which he showed lead to the formation of fully non-degenerate singularities
in which $\pmb{\partial}^2 \Phi$ blows up while $\Phi$ and $\pmb{\partial} \Phi$ remain bounded.
Roughly, Alinhac's framework
applied to initial data such the singular boundary has the strictly convex
structure depicted in Fig. <ref>,
and it allowed him to follow the solution up to the lowest point
on the singular boundary, but not further.
In Fig. <ref>,
we denote this lowest point by “$b_*$.”
Alinhac's framework has also been applied to other wave equations;
see, for example, bDiWhY2012,bDiWhY2015,dBiWyH2016.

§.§.§ Christodoulou's breakthrough on irrotational, isentropic shock formation

In <cit.>, Christodoulou proved a stable shock formation result
for open sets of irrotational and isentropic initial solutions to the relativistic
Euler equations in three spatial dimensions.
In this context, one can study the flow with the help of a potential function $\Phi$,
and relativistic Euler flow reduces to a quasilinear wave equation of type
(<ref>). Aside from a single exceptional equation of state corresponding
to the graph of a timelike minimal surface in Minkowski space, all wave equations
of irroational and isentropic relativistic fluid mechanics fail to satisfy the null condition,
i.e., the basic mechanism driving shock formation is present all equations but one.
The data that Christodoulou treated were compact perturbations of non-vacuum constant state data,
and his shock formation results revealed the instability of these states under irrotational and isentropic perturbations.
To study the solution, Christodoulou
used a refined version of nonlinear geometric optics,
based on techniques that he co-developed with Klainerman in their proof of the stability of the Minkowski space <cit.>.
Specifically, Christodoulou constructed an eikonal function $u$, i.e., a solution to the eikonal equation
$(\gfour^{-1})^{\alpha \beta}(\pmb{\partial} \Phi) \partial_{\alpha} u \partial_{\beta} u = 0$,
and he showed that the wave equation solution's first Cartesian derivatives,
$\pmb{\partial} \Phi$, remain quite smooth relative to the geometric
coordinates. The level sets $\nullhyparg{u}$ of $u$ are characteristic hypersurfaces for the wave equation.
The surfaces depend on the solution itself, reflecting the quasilinear nature of the flow.

As we mentioned already in Sect. <ref>,
as in the present paper, in Christodoulou's framework,
the formation of the shock corresponds to the vanishing of the inverse foliation density $\upmu$.
In the region of classical existence, one has $\upmu > 0$,
and when $\upmu \to 0$, the density of the $\nullhyparg{u}$ in Cartesian coordinate space
becomes infinite, signifying the piling up of the characteristics.
The blowup of $\pmb{\partial}^2 \Phi$ is a consequence of the degeneracy
between the Cartesian and a system[In <cit.>, the “angular” coordinate functions $\vartheta^1$
and $\vartheta^2$ were constructed so as to be constant along the integral curves of the null generator $\Lunit$.
In contrast, in the present paper, in the role of the “angular” coordinate functions, we use the standard
Cartesian coordinates $x^2$ and $x^3$, which are not typically constant along the integral curves of $\Lunit$.
This minor difference turns out to have no substantial effect on the analysis.
of “geometric coordinates” $(t,u,\vartheta^1,\vartheta^2)$,
where we schematically represent the degeneracy caused by the vanishing of $\upmu$
as follows: $\pmb{\partial} \sim \frac{1}{\upmu} \frac{\partial}{\partial u}$;
the blowup of $\pmb{\partial}^2 \Phi$
then follows from this relation,
from proving that $\upmu \to 0$ in finite time, and
from proving a lower of the form
$|\frac{\partial}{\partial u} \partial \Phi| \gtrsim 1$.

Importantly, Christodoulou's proof relied on geometric energy estimates relative to a foliation of spacetime
by the $\nullhyparg{u}$,
rather than Nash–Moser estimates.
To control the geometry without derivative loss, he relied on techniques and renormalized quantities
that have their roots in <cit.>.
Crucially, Christodoulou's geometric approach allowed him to prove shock formation for
a larger class of singularities than the fully non-degenerate ones treated by Alinhac.
He was able to show that blowup occurs at one or more points
even for solutions whose singular boundaries do not have to enjoy the strict convexity displayed in
Fig. <ref>, though in the absence of strict convexity, the full structure of the maximal
development was not revealed. Importantly, the work <cit.>
also yielded a sharp conditional global existence result,
which showed that irrotational, isentropic near-constant-state solutions are global
unless a shock forms.

§.§.§ Shock-formation works that built upon Christodoulou's framework

Many authors have used Christodoulou's approach to prove stable shock formation
for various quasilinear hyperbolic PDEs.
For example, there are stable shock formation results for:

* A larger class of wave equations sMpY2017,jS2016b.

* Solution regimes that are different than the small, compactly supported data regime:

* Solutions that exist classically precisely on a past-infinite half-slab <cit.>.

* Various systems involving multiple speeds of propagation,
some with symmetry dCdRP2016,xAhCsY2020,
and some without jS2018b,jS2019a.

§.§.§ Shock formation in the presence of vorticity and entropy

The aforementioned paper <cit.> was to first to prove stable shock formation for open
sets of $2D$ compressible Euler solutions with vorticity under an arbitrary[As in our main results,
the Chaplygin gas equation of state is exceptional and is not known to lead to shock formation.]
barotropic equation of state. The initial data were not required to be fully non-degenerate,
i.e., transversal convexity did not play a role
in the proof of blowup.
The main theorem followed the solution to the constant-time hypersurface of first blowup,
and it gave a complete description of what blows up and what does not, as well as a precise description of
the set of singular points at the time of first blowup,
i.e., in the context of Fig. <ref>,
a complete description of $\Sigma_{T_{\textnormal{Shock}}} \cap \mathcal{B}$.
Since the geometric setup relied on the methods of <cit.>, the
singularities were allowed to be more general than the
fully non-degenerate ones described in Sect. <ref>.

The recent work <cit.> extended <cit.>
to the case of the $3D$ compressible Euler equations under an arbitrary (non-Chaplygin gas)
equation of state with vorticity and entropy.
As in <cit.>, stable shock formation was proved in <cit.> without transversal convexity,
though in a sub-regime of solutions with transversal convexity, additional information
on the Hölder regularity of the solution with respect to the Cartesian coordinates was derived.[With
transversal convexity, the solutions was shown to enjoy $C^{1/3}$-Hölder regularity with
respect to the Cartesian coordinates up to the singularity, while not enjoying
$C^{(1/3})^+$-Hölder regularity. Providing a more detailed treatment of this regime
was inspired by Hölder regularity results derived
in tBsSvV2019a,tBsSvV2019b,tBsSvV2020 for some fully non-degenerate solutions.]

§.§.§ A new approach to proving the formation of fully-non-degenerate shock singularities via self-similarity

In tBsSvV2019b,tBsSvV2020, the authors provided a philosophically interesting new approach
for proving the formation of shock singularities in $3D$ compressible Euler solutions
under an adiabatic equation of state without symmetry,
and with vorticity (and also entropy in <cit.>).
The approach allows one to follow the solution to the time of first blowup,
and the singularities produced are isolated within the constant-time hypersurfaces
of first blowup. That is, in the context of
Fig. <ref>, the approach allows one
to follow the solution up to the point $b_*$.
Such singularities are analogs of the non-degenerate singularities that
Alinhac studied sA1999a,sA1999b,sA2002
in the case of quasilinear wave equations.
The framework of tBsSvV2019b,tBsSvV2020
relied on modulation parameters to show that
for open sets of smooth data, a singularity develops
in finite time, and it is a perturbation of a self-similar Burgers'-type shock.

The aforementioned $3D$ works were preceded by the work <cit.> in $2D$ azimuthal symmetry with vorticity.
In the recent work <cit.>, in the same symmetry class,
the authors constructed shock-forming solutions whose cusp-like spatial behavior (with respect to the standard coordinates)
is non-generic; such solutions are unstable.
In the language of the present paper,
such solutions do not exhibit the quantitative transversal convexity (<ref>)
that we use in proving our main results.

§.§.§ Self-similar blowup for non-hyperbolic PDEs

Notably, there are non-hyperbolic PDES with solutions that exhibit self-similar blowup
modeled on a self-similar Burgers shock. Examples include the Euler–Poison equations
Burgers' equation with transverse viscosity <cit.>,
the Burgers–Hilbert equations <cit.>, the fractal Burgers equation <cit.>,
and various dispersive or dissipative perturbations of the Burgers equation <cit.>;
see also cCteGsInM2018,cCteGnM2019.
It would be interesting to investigate the extent to which the set of singular points
can be understood, i.e., to try to derive results for these equations that are analogs of the results of the present paper.

§.§.§ Implosion singularities

The recent breakthrough work <cit.> in spherical symmetry showed that for
adiabatic equations of state $p = \varrho^{\upgamma}$ with $\upgamma > 1$,
there exist $C^{\infty}$ initial data such that the corresponding
solution's density and velocity blow up at the center of symmetry in finite time.
In fact, there are infinitely many such singularity-forming solutions,
collectively exhibiting a discrete sequence of blowup-rates.
These “implosion singularities” are much more severe than shocks.
The methods of <cit.> suggest that the implosion might enjoy co-dimension stability
under perturbations of the initial data without symmetry,
though possibly not full stability corresponding to
open sets of data.

§.§.§ Inviscid limits

A physically important and mathematically interesting problem is to study
the relationship between the formation of shock singularities in classical solutions to hyperbolic PDEs
and the behavior of classical solutions when a small amount of viscosity is added to the equations.
Of particular interest is to understand the zero viscosity limiting behavior of classical solutions.
The recent paper <cit.>, which concerns Burgers' equation in $1D$ and its viscous analog,
is the only work to date that yields information about the zero viscosity limit
of classical solutions all the way up until the time of first singularity formation for the inviscid solution.
The main results show that the $1D$ viscous Burgers' equation
solution can be decomposed into a singular piece and a smoother piece,
and that the viscous solution converges to the singular piece in $\| \cdot \|_{L^{\infty}}$ as the viscosity vanishes
(the $L^{\infty}$ norm is shown to be bounded from above by the viscosity parameter to a positive power).
Here, the $\| \cdot \|_{L^{\infty}}$ norm is over the entire region of classical existence of the inviscid solution,
i.e., the $L^{\infty}$-convergence holds all the way until the time of first blowup for the solution to the inviscid Burgers' equation. The results of <cit.>
apply for non-degenerate initial data, where the notion of non-degeneracy in <cit.>
is essentially equivalent to the transversal convexity assumption satisfied by the solutions featured in our main theorems.
An important open problem is to extend the results of <cit.> to the compressible Euler equations
(where the corresponding viscous model is the compressible Navier–Stokes equations)
and to multiple spatial dimensions.

§.§.§ Progress on the shock development problem

In Majda's celebrated works aM1981,aM1983a, he proved linear stability and local well-posedness
results for a class of weak solutions to the $3D$ compressible Euler equations arising from a set of
discontinuous initial data.
More precisely, he assumed that the data on $\mathbb{R}^3$ were piecewise smooth and jumped across
a smooth two-dimensional hypersurface such that the jumps are consistent with
well-known Rankine–Hugoniot jump conditions. He also assumed the data satisfy a certain “stability condition.”
His main result was the construction of a local, unique weak solution in a subset of $\mathbb{R}^{1+3}$,
and a corresponding three-dimensional shock hypersurface,
across which the solution jumps in accordance with the Rankine–Hugoniot jump conditions.
This is known as the shock front problem.

An important problem, distinct from the shock front problem,
is to determine whether/how Majda's discontinuous initial conditions
can develop from an initially smooth solution. That is, one would like to describe the transition of compressible Euler solutions
from being smooth, to developing a “first singularity” – which in the language of the present paper is the crease –
and finally to becoming a weak solution that develops a shock hypersurface (emanating from the crease),
across which the solution jumps. This is known as the shock development problem.
In Fig. <ref>, we denote the shock hypersurface by “$\mathcal{K}$,”
and we show its emergence from the crease $\partial_- \mathcal{B}$. We stress
that Majda's works did not study the flow in a neighborhood of $\partial_- \mathcal{B}$, but
rather started from initial conditions on a flat hypersurface $\Sigma_{t_0}$
such that the co-dimension surface $\Sigma_{t_0} \cap \mathcal{K}$, across which the data jump,
is already assumed to exist.
A crucial issue is that the jump conditions
imply that initially irrotational and isentropic smooth solutions will develop dynamic entropy and non-trivial vorticity as the solution
jumps across the shock hypersurface. This means that it is not possible to solve the true shock development problem in the class of irrotational and isentropic solutions. The full problem in $3D$ without symmetry assumptions remains open, but there has been
inspiring progress in recent years:

* In <cit.>, Christodoulou–Lisibach used Riemann invariants
and a pair of eikonal functions (one “ingoing” and the other “outgoing”)
to solve the shock development problem
for spherically symmetric, barotropic solutions to the relativistic Euler equations.
Mathematically, the problem comprised a system of quasilinear transport equations
coupled to the free boundary problem of tracking the location of the
shock hypersurface. The problem was solved through an iteration scheme that relied
on the assumption that the crease exhibits transversal convexity.
Relative to geometric coordinates analogous to the ones that we use in our main results,
the data are assumed to be smooth along the crease and Cauchy horizon, and this allows one to
use Taylor expansions to approximate the expected location
of the shock hypersurface and the behavior of the solution along it;
these Taylor expansions form a crucial ingredient in controlling the iterates.

* In his breakthrough monograph <cit.>, Christodoulou extended the methods from <cit.>
and solved the “restricted” shock development problem
without symmetry assumptions for the compressible Euler equations and relativistic Euler equations
in an arbitrary number of spatial dimensions. The word “restricted” means that he studied only irrotational and isentropic solutions,
and that 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 real one.
Christodoulou assumed that along the crease, the solution satisfies
the same kind of transversal convexity that the solutions produced by our main results enjoy.
In particular, this class of singularities is more general than the fully non-degenerate ones
treated by Alinhac in Sect. <ref>.
The main new difficulty in <cit.> is deriving energy estimates for the solution
and the acoustic geometry. In particular, relative to a system of geometric coordinates,
the high order energies exhibit degenerate behavior, much like the high order energies in the present paper.

* In <cit.>,
Huicheng–Lu solved the shock development problem
for a class of first-order, scalar, divergence-form hyperbolic equations in $2D$,
starting from initial singularities that are fully non-degenerate,
i.e., such that the crease is strictly convex,
as in Fig. <ref>.

* In <cit.>,
Buckmaster–Drivas–Shkoller–Vicol extended the methods from tBsSvV2019a,tBsSvV2019b,tBsSvV2020
to solve the shock development problem
for the $2D$ compressible Euler equations
in azimuthal symmetry with vorticity and entropy under the
adiabatic equations of state $p = \frac{1}{\upgamma} \varrho^{\upgamma} \exp(\Ent)$
for initial singularities that are fully non-degenerate
with respect to variations in the radial and angular variable.
In particular, <cit.> provides the first solution to a compressible Euler
shock development problem with vorticity.
As in <cit.>, the problem comprised transport equations coupled to the issue of tracking the
location of the free boundary, and it was solved via an iteration scheme based on Taylor expansions
relative to a coordinate system in which the solution is rather smooth.
A new feature compared to <cit.> is that there is a third characteristic
direction in <cit.>, corresponding to the transporting of vorticity and entropy by
the material derivative vectorfield.

We also highlight that recently, there have been other interesting works on weak solutions to the compressible
Euler equations in $1D$ without entropy.
In <cit.>, Lisibach proved local existence for the shock reflection (off of a wall)
problem in plane-symmetry.
In <cit.>, he studied the interaction of two shocks in plane-symmetry
and proved local existence of a weak solution near the interaction point.

§.§ Ideas behind the proof of the main results

In this section, we provide a more detailed overview of
some key ideas behind the proof of our main results.

§.§.§ Almost Riemann invariants

To study perturbations of simple isentropic plane-symmetric waves,
we find it convenient to replace the scalar functions $\LogDensity$
(where $\LogDensity$ is the logarithmic density from Def. <ref>)
and $v^1$ with the “almost[For plane-symmetric isentropic solutions,
$\RRiemann$, $\LRiemann$ are Riemann invariants, but away from symmetry, they are not.]
Riemann invariants” (see Def. <ref>)
$\RRiemann$, $\LRiemann$ which carry the same information.

Here and in the rest of the paper,
$\wavearray \eqdef (\RRiemann,\LRiemann,v^2,v^3,\Ent)$
are the fluid “wave variables”
and $\gfour_{\alpha \beta} = \gfour_{\alpha \beta}(\wavearray)$ denotes
the acoustical metric (see definition (<ref>)),
which drives the propagation of sound waves.

§.§.§ “Late-time” assumptions on the data for perturbations of simple isentropic plane-waves

We consider the “bona fide” initial data to be the data specified along
the flat spacelike hypersurface $\Sigma_0$ and a portion of a null hypersurface,
denoted by $\nullhyparg{- \rightu}^{[0,\frac{4}{\mathring{\updelta}_*}]}$ in Fig. <ref>
(i.e., we are studying the Cauchy problem for spacelike-characteristic data),
where the parameter $\mathring{\updelta}_* > 0$ (see (<ref>)) depends on the data and
is such that the Cartesian time of first blowup is approximately
However, the structures we need to detect the singular boundary appear only late
in the evolution, and we therefore find it convenient to
state our data-assumptions on a “late” hypersurface portion of constant rough time,
denoted by $\hypthreearg{\timefunction_0}{[- \rightu,\leftu]}{\muxmulevelsetvalue}$
in Fig. <ref>,
that will end up being close to the singular boundary $\mathcal{B}^{[0,\muxmulevelsetvalue_0]}$.
That is, we find it convenient
to state the assumptions on the spacelike level sets
\hypthreearg{\timefunction_0}{[- \rightu,\leftu]}{\muxmulevelsetvalue}
\eqdef
\lbrace \timefunctionarg{\muxmulevelsetvalue} = \timefunction_0 \rbrace \cap \lbrace - \rightu \leq u \leq \leftu \rbrace$
of the rough time functions,
which we describe in detail in Sect. <ref>.
In Sect. <ref>,
we state all of these data-assumptions.
In Appendix <ref>, we use
Cauchy stability arguments to show that these assumptions
are satisfied by perturbations of the bona fide data on $\Sigma_0$
corresponding to a class of shock-forming simple isentropic plane-wave solutions.
In Appendix <ref>, we use standard arguments to construct these
simple isentropic plane-wave solutions.

[scale=.4, grid = false, tics=5]Datahypersurfacesgeometric.pdf
(42,76) $\twoargmumuxtorus{0}{-\muxmulevelsetvalue}$
(32,71) $\mathcal{B}^{[0,\muxmulevelsetvalue_0]}$
(38,15) $\Sigma_0$
(4,5) $\ell_{0,\leftu}$
(24,5) $\ell_{0,u}$
(61,5) $\ell_{0,-\rightu}$
(26,25) $\nullhyparg{u}$ (truncated)
(21,45) $\twoargroughtori{\timefunction_0,u}{\muxmulevelsetvalue}$
(32,-2.5) $(x^2,x^3) \in \mathbb{T}^2$
(41,6) $t$
(31,3) $u \in\mathbb{R}$
(10,33) $\hypthreearg{\timefunction_0}{[- \rightu,\leftu]}{\muxmulevelsetvalue}$
(10,50) $\hypthreearg{0}{[- \rightu,\leftu]}{\muxmulevelsetvalue}$
(64,80) $\nullhyparg{- \rightu}^{[0,\frac{4}{\mathring{\updelta}_*}]}$

Data hypersurfaces in geometric coordinate space

§.§.§ The geometric div-curl-transport formulation of the flow

To derive estimates for the fluid variables, we fundamentally rely on
a geometric reformulation of the flow, which was derived in
<cit.> and which we restate below in Theorem <ref>.
To aid our discussion in the remainder of the introduction,
here we recall the formulation in schematic form:

\begin{align} \label{E:WAVEINTRO}
\square_{\gfour(\wavearray)} \Psi
\nullform(\pmb{\partial} \wavearray, \pmb{\partial} \wavearray)
(\vortrenormalized,\GradEnt) \cdot \pmb{\partial} \wavearray
(\VortVort, \DivGradEnt),
\\
\Transport (\vortrenormalized,\GradEnt)
& =  (\vortrenormalized,\GradEnt) \cdot \pmb{\partial} \wavearray,
\label{E:TRANSPORTINTRO} \\
\Transport (\VortVort,\DivGradEnt)
& = \nullform(\pmb{\partial} \wavearray,\partial \vortrenormalized)
\nullform(\pmb{\partial} \wavearray,\partial \GradEnt)
\GradEnt \cdot \nullform(\pmb{\partial} \wavearray, \pmb{\partial} \wavearray)
\GradEnt \cdot \GradEnt \cdot \pmb{\partial} \wavearray,
\label{E:TRANSPORTFORMODIFIEDINTRO} \\
(\Flatdiv \vortrenormalized, \Flatcurl \GradEnt)
& = \partial \wavearray.
\label{E:DIVCURLFORMODIFIEDINTRO}
\end{align}

In (<ref>)–(<ref>),
\vortrenormalized = \frac{(\Flatcurl v)^i}{\exp (\LogDensity)}