The relativistic Euler equations with a physical vacuum boundary: Hadamard local well-posedness, rough solutions, and continuation criterion

In this paper we provide a complete local well-posedness theory for the free boundary relativistic Euler equations with a physical vacuum boundary on a Minkowski background. Specifically, we establish the following results: (i) local well-posedness in the Hadamard sense, i.e., local existence, uniqueness, and continuous dependence on the data; (ii) low regularity solutions: our uniqueness result holds at the level of Lipschitz velocity and density, while our rough solutions, obtained as unique limits of smooth solutions, have regularity only a half derivative above scaling; (iii) stability: our uniqueness in fact follows from a more general result, namely, we show that a certain nonlinear functional that tracks the distance between two solutions (in part by measuring the distance between their respective boundaries) is propagated by the flow; (iv) we establish sharp, essentially scale invariant energy estimates for solutions; (v) a sharp continuation criterion, at the level of scaling, showing that solutions can be continued as long as the velocity is in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^1_t Lip$$\end{document}Lt1Lip and a suitable weighted version of the density is at the same regularity level. Our entire approach is in Eulerian coordinates and relies on the functional framework developed in the companion work of the second and third authors on corresponding non relativistic problem. All our results are valid for a general equation of state \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p(\varrho )= \varrho ^\gamma $$\end{document}p(ϱ)=ϱγ, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma > 1$$\end{document}γ>1.


Introduction
In this article, we consider the relativistic Euler equations, which describe the motion of a relativistic fluid in a Minkowski background M d+1 , d ≥ 1. The fluid state is represented by the (energy) density ≥ 0, and the relativistic velocity u. The velocity is assumed to be a forward time-like vector field, normalized by u α u α = −1. (1.1) The equations of motion consist of where T is the energy-momentum tensor for a perfect fluid, defined by T αβ :=( p + )u α u β + p m αβ . (1.3) Here m is the Minkowski metric, and p is the pressure, which is subject to the equation of state p = p( ).
Projecting (1.2) onto the directions parallel and perpendicular to u, using definition (1.3), and the identity (1.1), yields the system u μ ∂ μ + ( p + )∂ μ u μ = 0 ( p + )u μ ∂ μ u α + μ α ∂ μ p = 0, (1.4) with u satisfying the constraint (1.1), which is in turn preserved by the time evolution. Here is is the projection on the space orthogonal to u and is given by Throughout this paper, we adopt standard rectangular coordinates in Minkowski space, denoted by {x 0 , x 1 , . . . , x d }, and we identify x 0 with a time coordinate, t:=x 0 . Greek indices vary from 0 to d and Latin indices from 1 to d.
The system (1.4) can be seen as a nonlinear hyperbolic system, which in the reference frame of the moving fluid has the propagation speed which is subject to 0 ≤ c s < 1, implying that the speed of propagation of sound waves is always non-negative and below the speed of light (which equals to one in the units we adopted).
In this article we consider the physical situation where vacuum states are allowed, i.e. the density is allowed to vanish. The gas is located in the moving domain t :={x ∈ R d | (t, x) > 0}, whose boundary t is the vacuum boundary, which is advected by the fluid velocity u.
The distinguishing characteristic of a gas, versus the case of a liquid, is that the density, and implicitly the pressure and the sound speed, vanish on the free boundary t , = 0, p = 0, c s = 0 on t .
Thus, the equations studied here provide a basic model of relativistic gaseous stars (see Section 1. 6). An appropriate equation of state to describe this situation is 1 (see,e.g., [[37], Section 2.4] or [35]): The decay rate of the sound speed at the free boundary plays a critical role. Precisely, there is a unique, natural decay rate which is consistent with the time 1 Observe that the requirement 0 ≤ c 2 s < 1 imposes a bound on . This occurs because power-law equations of state such as (1.5) are no longer valid if the density is very large [17]. evolution of the free boundary problem for the relativistic Euler gas, which is commonly referred to as physical vacuum, and has the form (1.6) where dist(·, ·) is the distance function. Exactly the same requirement is present in the non-relativistic compressible Euler equations. As in the non-relativistic setting, (1.6) should be considered as a condition on the initial data that is propagated by the time-evolution. There are two classical approaches in fluid dynamics, using either Eulerian coordinates, where the reference frame is fixed and the fluid particles are moving, or using Lagrangian coordinates, where the particles are stationary but the frame is moving. Both of these approaches have been extensively developed in the context of the Euler equations, where the local well-posedness problem is very well understood.
By contrast, the free boundary problem corresponding to the physical vacuum has been far less studied and understood. Because of the difficulties related to the need to track the evolution of the free boundary, all the prior work is in the Lagrangian setting and in high regularity spaces which are only indirectly defined.
Our goal in this paper is to provide the first local well-posedness result for this problem. Unlike previous approaches, which were limited to proving energy-type estimates at high regularity and in a Lagrangian setting [12,16], here we consider this problem fully within the Eulerian framework, where we provide a complete local well-posedness theory, in the Hadamard sense, in a low regularity setting. We summarize here the main features of our result, which mirror the results in the last two authors' prior paper devoted to the non-relativistic problem [14], referring to Section 1.5 for precise statements: a) We prove the uniqueness of solutions with very limited regularity v ∈ Li p, ∈ Li p 2 . More generally, at the same regularity level we prove stability, by showing that bounds for a certain nonlinear distance between different solutions can be propagated in time. b) Inspired by [14], we set up the Eulerian Sobolev function space structure where this problem should be considered, providing the correct, natural scale of spaces for this evolution. c) We prove sharp, scale invariant 3 energy estimates within the above mentioned scale of spaces, which guarantee that the appropriate Sobolev regularity of solutions can be continued for as long as we have uniform bounds at the same scale v ∈ Li p. d) We give a constructive proof of existence for regular solutions, fully within the Eulerian setting, based on the above energy estimates. e) We employ a nonlinear Littlewood-Paley type method, developed prior work [14], in order to obtain rough solutions as unique limits of smooth solutions. This also yields the continuous dependence of the solutions on the initial data.

Space-time foliations and the material derivative
The relativistic character of our problem implies that there is no preferred choice of coordinates. On the other hand, in order to derive estimates and make quantitative assertions about the evolution, we have to choose a foliation of spacetime by space-like hypersurfaces. Here, we take advantage of the natural set-up provided by Minkowski space and foliate the spacetime by {t = constant} slices. We then define the material derivative, which is adapted to this specific foliation, as The vectorfield D t is better adapted to the study of the free-boundary evolution than working directly with u μ ∂ μ . Indeed, in order to track the motion of fluid particles on the boundary, we need to understand their velocity relative to the aforementioned spacetime foliation. The velocity that is measured by an observer in a reference frame characterized by the coordinates (t, x 1 , . . . , x d ) is u i /u 0 . This is a consequence of the fact that in relativity observers are defined by their world-lines, which can be reparametrized. This ambiguity is fixed by imposing the constraint u μ u μ = −1. As a consequence, the d-dimensional vectorfield (u 1 , . . . , u d ) can have norm arbitrarily large, while the physical velocity has to have norm at most one (the speed of light). It follows, in particular, that fluid particles on the boundary move with velocity u i /u 0 . These considerations also imply that the standard differentiation formula for moving domains holds with D t , i.e., (1.8) This formula remains valid with the good variable v we introduce below since v i /v 0 = u i /u 0 .

The good variables
The starting point of our analysis is a good choice of dynamical variables. We seek variables that are tailored to the characteristics of the Euler flow all the way to moving boundary, where the sound characteristics degenerate due to the vanishing of the sound speed. Our choice of good variables will (i) better diagonalize the system with respect to the material derivative, (ii) be associated with truly relativistic properties of the vorticity, and (iii) lead to good weights that allow us to control the behavior of the fluid variables when one approaches the boundary.
Property (i) will be intrinsically tied with both the wave and transport character of the flow in that (a) the diagonalized equations lead to good second order equations that capture the propagation of sound in the fluid, see Section 3.2, and (b) it provides a good transport structure that will allow us to implement a time discretization for the construction of regular solution, see Section 6. Property (ii) will ensure a good coupling between the wave-part and the transport-part of the system. Finally, property (iii) will lead to the correct functional framework needed to close the estimates. 4 Our good variables, denoted by (r, v), are defined in (1.9) and (1.15). The corresponding equations of motion are (1.16), which we now derive.
Our first choice of good variables is a rescaled version of the velocity given by where f is given by Although we are interested in the case p( ) = κ+1 , it is instructive to consider first a general barotropic equation of state; see the discussion related to the vorticity further below. In order to understand our choice for f , compute Solving for ∂ μ u α and plugging the resulting expression into the second equation of (1.4) we find We see that the term in parenthesis vanishes if f is given by (1.10), resulting in an equation which is diagonal with respect to the material derivative, and which we write as (1.11) We notice that in terms of v, the material derivative (1.7) reads as In view of the constraint (1.1), we have that v 0 satisfies 12) and in solving for v 0 we chose the positive square root because u, and thus v, is a future-pointing vectorfield. We now show that our choice (1.9) also diagonalizes the first equation in (1.4). First, we use (1.11) with α = 0 and solve to ∂ t v 0 , obtaining where in the second equality we used (1.12) to compute ∂ i v 0 . Using the above identity for ∂ t v 0 , we find the following expression for ∂ μ v μ : where δ is the Euclidean metric. Expressing ∂ μ u μ in terms of ∂ μ v μ (and derivatives of ) and using the above expression for ∂ μ v μ , we see that the first equation in (1.4) can be written as (1.13) Here we are using the notation (1.14) Observe that Eqs. (1.11) and (1.13) are valid for a general barotropic equation of state. We now assume the equation of state (1.5). Then the sound speed is given by c 2 s = (κ + 1) κ and f becomes f ( ) = (1 + κ ) 1+ 1 κ (we choose the constant of integration by setting f (0) = 1, so that v = u when = 0). It turns out that it is better to adopt the sound speed squared as a primary variable instead of because it plays the role of the correct weight in our energy functionals. We thus define 5 the second component of our good variables by Therefore, using (r, v) as our good variables, and p( ) given by (1.5) we find that the Eqs. (1.11) and (1.13) become 5 The factor 1 κ in the definition of r is a matter of convenience. Although r and c 2 s differ by this factor, we slightly abuse the terminology and also call r the sound speed squared.
where we have defined and the coefficients a 0 , a 1 and a 2 are given by Equations (1.16) are the desired diagonal with respect to D t equations, and the rest of the article will be based on them. In writing these equations we consider only the spatial components v i as variables, with v 0 always given by The specific form of the coefficients a 0 , a 1 , and a 2 is not very important for our argument. We essentially only use that they are smooth functions of r and v, and that a 0 , a 2 > 0. The operator G i j ∂ i (·) j can be viewed as a divergence type operator. This divergence structure is related to the fact that Equations (1.16) express the wave-like behavior of r and of the divergence part of v. The symmetric and positive-definite matrix c 2 s G i j is closely related to the inverse of the acoustical metric; precisely, they agree at the leading order near the boundary.
As we will see, Equations (1.16) also have the correct balance of powers of r to allow estimates all the way to the free boundary. The r factor in the divergence of v is related to the propagation of sound in the fluid (see Section 3.2) whereas the r factor in the last term of (1.16a) will allow us to treat ra 1 v i ∂ i r essentially as a perturbation at least in elliptic estimates (see Section 5).
One can always diagonalize Eq. (1.4) by simply algebraically solving for ∂ t ( , u). But it is not difficult to see that this procedure will not lead to equations with good structures for the study of the vacuum boundary problem. In this regard, observe that the choice (1.9) is a nonlinear change of variables, whereas algebraically solving for ∂ t ( , u) is a linear procedure.
We now comment on the relation between v and the vorticity of the fluid ω. It is well-known (see, e.g., [5] Section IX.10.1) that in relativity the correct notion of vorticity is given by the following two-form in spacetim where d st is the exterior derivative in spacetime. This is true not only for the power law equation of state (1.5), but also for an arbitrary barotropic equation of state. A computation using (1.18) (see, e.g., [5]) and the equations of motion implies that v α ω αβ = 0, (1.19) and that ω satisfies the following evolution equation Observe that (1.20) implies that ω = 0 if it vanishes initially. Since we will consider only the spatial components of v as independent, we use (1.19) to eliminate the 0 j components of ω from (1.20) as follows: from (1.19) we can write (1.22) Equation (1.22) will be used to derive estimates for ω i j that will complement the estimates for r and the divergence of v obtained from (1.16).
We remark that in the literature, the use of v, given by (1.9), seems to be restricted mostly to definition and evolution of the vorticity. To the best of our knowledge, this is the first time when it was observed that the same change of variables needed to define the relativistic vorticity also diagonalizes the equations of motion with respect to D t .

Scaling and bookkeeping scheme
Although Eq. (1.16) do not obey a scaling law, it is still possible to identify a scaling law for the leading order dynamics near the boundary. This will motivate the control norms we introduce in the next section, as well as provide a bookkeeping scheme that will allow us to streamline the analysis of many complex multilinear expressions we will encounter.
As we will see, the contribution of last term in (1.16a) to our energies is negligible, due to the multiplicative r factor. Thus, we ignore this term for our scaling analysis. 6 Replacing all coefficients that are functions of (r, v) by 1, while keeping the transport and divergence structure present in the equations, we obtain the following simplified version of (1.16): This system is expected to capture the leading order dynamics near the boundary, and also mirrors the nonrelativistic version of the compressible Euler equations, considered in the predecessor to this paper, see [14]. Equations (1.23) admit the scaling law Based on this leading order scaling analysis, we assign the following order to the variables and operators in Eq. (1.16): (i) r and v have order −1 and −1/2, respectively. More precisely, we only count v as having order −1/2 when it is differentiated. Undifferentiated v's have order zero. (ii) D t and ∂ i have order 1/2 and 1, respectively. (iii) G, a 0 , a 1 , and a 2 , and more generally, any smooth function of (r, v) not vanishing at r = 0, have order 0.
Expanding on (iii) above, the order of a function of r is defined by the order of its leading term in the Taylor expansion about r = 0, being of order zero if this leading term is a constant. The order of a multilinear expression is defined as the sum of the orders of each factor. Here we remark that all expressions arising in this paper are multilinear expressions, with the possible exception of nonlinear factors as in (iii) above. According to this convention, all terms in equation (1.16b) have order zero, and all terms in (1.16a) have order −1/2, except for the last term in (1.16a) which has order −1. Upon successive differentiation of any multilinear expression with respect to D t or ∂, all terms produced are the same (highest) order, unless some of these derivatives apply to nonlinear factors as in (iii); then lower order terms are produced.

Energies, function spaces, and control norms
Here we introduce the function spaces and control norms that we need in order to state our main results. A more detailed discussion is given in Section 2. With some obvious adjustments, here we follow the lat two authors' prior work in [14]. We assume throughout that r is a positive function on t , vanishing simply on the boundary, and so that r is comparable to the distance to the boundary t .
In order to identify the correct functional framework for our problem, we start with the linearization of the Eq. (1.16). In Section 3 we show that the linearized equations admit the following energy which defines the (time dependent) weighted L 2 space H. The motivation for the definition of higher order norms and spaces comes from the good second order equations mentioned in Section 1.2. From Eq. (1.16), we find that the second order evolution is governed at leading order by a wave-like operator which is essentially a variable coefficient version of D 2 t − r . This points toward higher order spaces built on powers of r . Taking into account also the form of the linearized energy above, we are led to the following. We define H 2k as the space of pairs of functions (s, w) in t for which the norm below is finite The definition of H 2k for non-integer k is given in Section 2, via interpolation.
In view of the scaling analysis of Section 1.3, we introduce the critical space H 2k 0 where (1.24) which has the property that its leading order homogeneous component is invariant with respect to the scaling discussed in Section 1.3. Associated with the exponent 2k 0 we define the following scale invariant time dependent control norm here N is a given non-zero vectorfield with the following property. In each sufficiently small neighborhood of the boundary, there exists a x 0 ∈ t such that N (x 0 ) = ∇r (x 0 ). The fact that we can choose such a N follows from the properties of r . The motivation for introducing N is that we can make A small by working in small neighborhood of each reference point x 0 , whereas ∇r L ∞ is a scale invariant quantity that cannot be made small by localization arguments. We further introduce a second time dependent control norm that is associated with H 2k 0 +1 , given by 7 It follows that ∇r C 1 2 scales like theĊ 3 2 norm of r , but it is weaker in that it only uses one derivative of r away from the boundary. The norm B will control the growth of our energies, allowing for a secondary dependence on A.
When the density is bounded away from zero, the relativistic Euler equations can be written as a first-order symmetric hyperbolic system (see, e.g., [1]) and standard techniques can be applied to derive local estimates. The difficulties in our case come from the vanishing of r on the boundary. Using the finite speed of propagation of the Euler flow, we can use a partition of unity to separate the nearboundary behavior, where r approaches zero, from the bulk dynamics, where r is bounded away from zero. Furthermore, we can also localize to a small set where A is small. Such a localization will be implicitly assumed in all our analysis, in order to avoid cumbersome localization weights through the proofs.

The main results
Here we state our main results. Combined, these results establish the sharp local well-posedness and continuation criterion discussed earlier. We will make all our statements for the system written in terms of the good variables (r, v), i.e., Eq. (1. 16). Readers interested in the evolution of (1.4) should have no difficulty translating our statements to the original variables and u.
We recall that Eq. (1.16) are always considered in the moving domain given by for some T > 0, where the moving domain at time t, t , is given by We also recall that we are interested in solutions satisfying the physical vacuum boundary condition x .
For the next Theorem, we introduce the phase space We refer to Section 2 for a more precise definition of H 2k , including its topology.
Since the H 2k norms depend on r , it is appropriate to think of H 2k in a nonlinear fashion, as an infinite dimensional manifold. We also stress that, while k was an integer in our preliminary discussion in Section 1.4, in Section 2 we extend their definition for any k ≥ 0. Consequently, H 2k is also defined for any k ≥ 0, and our Theorems 1.2 and 1.4 below include non-integer values of k. where k 0 is given by (1.24).
• Uniqueness of solutions in a larger class, see Theorem 1.1.
• Continous dependence of solutions on the initial data in the H 2k topology.
Furthermore, in our proof of uniqueness in Section 4 we establish something stronger, namely, that a suitable nonlinear distance between two solutions is propagated under the flow. This distance functional, in particular, tracks the distance between the boundaries of the moving domains associated with different solutions. Thus, our local well-posedness also includes: • Weak Lipschitz dependence on the initial data relative to a suitable nonlinear functional introduced in Section 4.
An important threshold for our results corresponds to the uniform control parameters A and B. Of these A is at scaling, while B is one half of a derivative above scaling. Thus, by Lemma 2.5 of Section 2, we will have the bounds Next, we turn our attention to the continuation of solutions.
where C(A) is a constant depending on A. The energies E 2k will be constructed explicitly only for integer k. Nevertheless, our analysis will show that (1.28) will also hold for any k > 0. This will be done using a mechanism akin to a paradifferential expansion, without explicitly constructing energy functionals for non-integer k. As a consequence, we will obtain Theorem 1.4. Let k be as in (1.27). Then, the H 2k solution given by Theorem 1.2 can be continued as long as A remains bounded and B ∈ L 1 t ( ).

Historical comments
The study of the relativistic Euler equations goes back to the early days of relativity theory, with the works of Einstein [7] and Schwarzschild [38]. The relativistic free-boundary Euler equations were introduced in the '30s in the classical works of Tolman, Oppenheimer, and Volkoff [34,39,40], where they derived the now-called TOV equations. 8 With the goal of modeling a star in the framework of relativity, Tolman, Oppenheimer, and Volkoff studied spherically symmetric static solutions to the Einstein-Euler system for a fluid body in vacuum and identified the vanishing of the pressure as the correct physical condition on the boundary. Observe that such a condition covers both the cases of a liquid, where > 0 on the boundary, as well as a gas, which we study here, where = 0 on the boundary. This distinction is related to the choice of equation of state.
Although the TOV equations have a long history and the study of relativistic stars is an active and important field of research (see, e.g., [37], Part III and [30], Part V), the mathematical theory of the relativistic free-boundary Euler equations lagged behind.
If we restrict ourselves to spherically-symmetric solutions, possibly also considering coupling to Einstein's equations, a few precise and satisfactory mathematical statements can be obtained. Lindblom [20] proved that a static, asymptotically flat spacetime, that contains only a uniform-density perfect fluid confined to a spatially compact region ought to be spherically symmetric, thus generalizing to relativity a classical result of Carleman [4] and Lichtenstein [19] for Newtonian fluids. The proof of existence of spherically symmetric static solutions to the Einstein-Euler system consisting of a fluid region and possibly a vacuum region was obtained by Rendall and Schmidt [36]. Their solutions allow for the vanishing of the density along the interface of the fluid-vacuum region, although it is also possible that the fluid occupies the entire space and the density merely approaches zero at infinity. Makino [21] refined this result by providing a general criterion for the equation of state which ensures that the model has finite radius. Makino has also obtained solutions to the Einstein-Euler equations in spherical symmetry with a vacuum boundary and near equilibrium in [22,23], where equilibrium here corresponds to the states given by the TOV equations. In [24], Makino extended these results to axisymmetric solutions that are slowly rotating, i.e., when the speed of light is sufficiently large or when the gravitational field is sufficiently weak (see also the follow-up works [25,26] and the preceding work in [13]). Another result within symmetry class related to the existence of vacuum regions and relevant for the mathematical study of star evolution is Hadžić and Lin's recent proof of the "turning point principle" for relativistic stars [11].
The discussion of the last paragraph was not intended to be an exhaustive account of the study of the relativistic free-boundary Euler equations under symmetry assumptions, and we refer the reader to the above references for further discussion. Rather, the goal was to highlight that a fair amount of results can be obtained in symmetry classes. This is essentially because some of the most challenging aspects of the problem are absent or significantly simplified when symmetry is assumed. This should be contrasted with what is currently known in the general case, which we now discuss.
Local existence and uniqueness of solutions the relativistic Euler equations in Minkowski background with a compactly supported density have been obtained by Makino and Ukai [27,28] and LeFloch and Ukai [18]. These solutions, however, require some strong regularity of the fluid variables near the free boundary and, in particular, do not allow for the existence of physical vacuum states. Similarly, Rendall [35] established a local existence and uniqueness 9 result for the Einstein-Euler system where the density is allowed to vanish. Nevertheless, as the author himself pointed out, the solutions obtained are not allowed to accelerate on the free boundary and, in particular, do not include the physical vacuum case. Rendall's result has been improved by Brauer and Karp [2,3], but still without allowing for a physical vacuum boundary. Oliynyk [31] was able to construct solutions that can accelerate on the boundary, but his result is valid only in one spatial dimension. A new approach to investigate the free-boundary Euler equations, based on a frame formalism, has been proposed by Friedrich in [8] (see also [9]) and further investigated by the first author in [6], but it has not led to a local well-posedness theory.
In the case of a liquid, i.e., where the fluid has a free-boundary where the pressure vanishes but the density remains strictly positive, a-priori estimates have been obtained by Ginsberg [10] and Oliynyk [32]. Local existence of solutions was recently established by Oliynyk [33] whereas Miao, Shahshahani, and Wu [29] proved local existence and uniqueness for the case when the fluid is in the so-called hard phase, i.e., when the speed of sound equals to one. See also [41], where the author, after providing a proof of local existence for the non-isentropic compressible free-boundary Euler equations in the case of a liquid, discusses ideas to adapt his proof for the relativistic case.
Finally, for the case treated in this paper, i.e., the relativistic Euler equations with a physical vacuum boundary, the only results we are aware of are the a-priori estimates by Hadžić, Shkoller, and Speck [12] and Jang, LeFloch, and Masmoudi [15]. In particular, no local existence and uniqueness (let along a complete local well-posedness theory as we present here) had been previously established.

Outline of the paper
Our approach carefully considers the dual role of r , on the one hand, as a dynamical variable in the evolution and, on the other hand, as a defining function of the domain that, in particular, plays the role of a weight in our energies. An important aspect of our approach is to decouple these two roles. Such decoupling is what allows us to work entirely in Eulerian coordinates. When comparing different solutions (which in general will be defined in different domains), we can think of the role of r as a defining function as leading to a measure of the distance between the two domains (i.e., a distance between the two boundaries), whereas the role of r as a dynamical variable leads to a comparison in the common region defined by the intersection of the two domains. For instance, in our regularization procedure for the construction of regular solutions, the defining functions of the domains are regularized at a different scale than the main dynamical variables.
Although the relativistic and non-relativistic Euler equations, and their corresponding physical vacuum dynamics, are very different, some of our arguments here will closely follow those in the last two authors' prior work [14], where results similar to those of Section 1.5 were established for the non-relativistic Euler equations in physical vacuum. Thus, when it is appropriate, we will provide a brief proof, or quote directly from [14]. This is particularly the case for Sects. 6 and 7.
The paper is organized as follows: is key to our analysis. Similar scales of spaces have been introduced in [14] treating the non-relativistic case and also in [16] where the non-relativistic problem had been considered in Lagrangian coordinates and in high regularity spaces.
Our function spaces H 2k are Sobolev-type spaces with weights r . Since the fluid domain is determined by t := {r > 0}, the state space H 2k is nonlinear, having a structure akin to an infinite dimensional manifold.
Interpolation plays two key roles in our work. Firstly, it allows us to define H 2k for non-integer k without requiring us to establish direct energy estimates with fractional derivatives. This is in particular important for our low regularity setting since the critical exponent (1.24) will in general not be an integer. Secondly, we interpolate between H 2k and the control norms A and B. For this we use some sharp interpolation inequalities presented in Section 2.3. These inequalities are proven in the last two authors' prior work [14] and, to the best of our knowledge, have not appeared in the literature before. In fact, it is the use of these inequalities that allows us to work at low regularity, to obtain sharp energy estimates, and a continuation criterion at the level of scaling.

The linearized equation and the corresponding transition operators, Sect. 3
The linearized equation and its analysis form the foundation of our work, rather than direct nonlinear energy estimates. Besides allowing us to prove nonlinear energy estimates for single solutions, basing our analysis on the linearized equation will also allow us to get good quantitative estimates for the difference of two solutions. The latter is important for our uniqueness result and for the construction of rough solutions as limits of smooth solutions. We observe that there are no boundary conditions that need to be imposed on the linearized variables. This is related to the aforementioned decoupling of the roles of r and signals a good choice of functional framework.
Using the linearized equation we obtain transition operators L 1 and L 2 that act at the level of the linearized variables s and w. These transition operators are roughly the leading elliptic part of the wave equations for s and the divergence part of w. Note that since the wave evolution for the fluid degenerates on the boundary due to the vanishing of the sound speed, so do the transition operators L 1 and L 2 . We refer to L 1 and L 2 as transition operators because they relate the spaces H 2k+2 and H 2k in a coercive, invertible manner. Because of that, these operators play an important role in our regularization scheme used to construct high-regularity solutions.

Difference estimates and uniqueness, Sect. 4
In this section we construct a nonlinear functional that allows us to measure the distance between two solutions. We show that bounds for this functional are propagated by the flow, which in particular implies uniqueness. A fundamental difficulty is that, since we are working in Eulerian coordinates, different solutions are defined in different domains. This difficulty is reflected in the nonlinear character of our functional, which could be thought of as measuring the distance between the boundaries of two different solutions. The low regularity at which we aim to establish uniqueness leads to some technical complications that are dealt with by a careful analysis of the problem.

Energy estimates and coercivity, Sect. 5
The energies that we use contain two components, a wave component and a transport component, in accordance with the wave-transport character of the system. The energy is constructed after identifying Alinhac-type "good variables" that can be traced back to the structure of the linearized problem. This connection with the linearized problem is also key to establish the coercivity of the energy in that it relies on the transition operators L 1 and L 2 mentioned above.

Existence of regular solutions, Sect. 6
This section establishes the existence of regular solutions. It heavily relies on the last two authors' prior work [14], to which the reader is referred for several technical points.
Our construction is based essentially on an Euler scheme to produce good approximating solutions. Nevertheless, a direct implementation of Euler's method loses derivatives. We overcome this by preceding each iteration with a regularization at an appropriate scale and a separate transport step. The main difficulty is to control the growth of the energies at each step.

Rough solutions as limits of regular solutions, Sect. 7
In this section we construct rough solutions as limits of smooth solutions, in particular establishing the existence part of Theorem 1.2. We construct a family of dyadic regularizations of the data, and control the corresponding solutions in higher H 2k norms with our energy estimates, and the difference of solutions in H with our nonlinear stability bounds. The latter allow us to establish the convergence of the smooth solutions to the desired rough solution in weaker topologies. Convergence in H 2k is obtained with more accurate control using frequency envelopes. A similar argument then also gives continuous dependence on the data.

Notation for v, ω and the use of Latin indices
In view of Eq. (1.16) and the corresponding vorticity evolution (1.22), we have now written the dynamics solely in terms of r and the spatial components of v, i.e., v i . We henceforth consider v as a d-dimensional vector field, so that whenever referring to v we always mean (v 1 , . . . , v d ). v 0 is always understood as a shorthand for the RHS of (1.17). Similarly, by ω will stand for ω i j .
Recalling that indices are raised and lowered with the Minkowski metric and that m 0i = 0 = m 0i , m i j = δ i j , we see that tensors containing only Latin indices have indices equivalently raised and lowered with the Euclidean metric.

Function spaces
Here we define the function spaces that will play a role in our analysis. They are weighted spaces with weights given by the sound speed squared r which, in view of (1.25), is comparable to the distance to the boundary. More precisely, since a solution to (1.16) is not a-priori given, in the definitions below we take r to be a fixed non-degenerate defining function for the domain t , i.e., proportional to the distance to the boundary t . In turn, the boundary t is assumed to be Lipschitz.
We denote the L 2 -weighted spaces with weights h by L 2 (h) and we equip them with the norm With these notations the base L 2 space of pairs of functions in for our system, denoted by H, is defined as This space depends only on the choice of r . However, we will often use an equivalent norm that also depends on v, which corresponds to the energy space for the linearized problem and will also be important in the construction of our energies: This uses G to measure the pointwise norm of the one form w. The H norm is equivalent to the H 0 norm (see the definition of H 2k below) since G is equivalent to the the Euclidean inner product with constants depending on the L ∞ norm of (r, v). We continue with higher Sobolev norms. We define H j,σ , where j ≥ 0 is an integer and σ > − 1 2 , to be the space of all distributions in t whose norm is finite. Using interpolation, we extend this definition, thus defining H s,σ for all real s ≥ 0.
To measure higher regularity we will also need higher Sobolev spaces where the weights depend on the number of derivatives. More precisely, we define H 2k as the space of pairs of functions (s, w) defined inside t , and for which the norm below is finite : We extend the definition of H 2k to non-integer k using interpolation. An explicit characterization of H 2k for non-integer k, based on interpolation, was given in the last two authors' prior work [14]. Using the embedding theorems given below, we can show that the H 2k norm is equivalent to the H 2k, 1−κ 2κ +k × H 2k, 1 2κ +k norm.

The state space H 2k .
As already mentioned in the introduction, the state space H 2k is defined for k > k 0 (i.e. above scaling) as the set of pairs of functions (r, v) defined in a domain t in R d with boundary t with the following properties: a) Boundary regularity: t is a Lipschitz surface. b) Nondegeneracy: r is a Lipschitz function in¯ t , positive inside t and vanishing simply on the boundary t . c) Regularity: The functions (r, v) belong to H 2k .
Since the domain t itself depends on the function r , one cannot think of H 2k as a linear space, but rather as an infinite dimensional manifold. However, describing a manifold structure for H 2k is beyond the purposes of our present paper, particularly since the trajectories associated with our flow are merely expected to be continuous with values in H 2k . For this reason, here we will limit ourselves to defining a topology on H 2k .
ii) Domain convergence, r n − r Li p → 0. Here, we consider the functions r n and and r as extended to zero outside their domains, giving rise to Liptschitz functions in R d . iii) Norm convergence: for each > 0 there exist a smooth function (r ,ṽ) in a neighbourhood of so that The last condition in particular provides both a uniform bound for the sequence (r n , v n ) in H 2k ( n ) as well as an equicontinuity type property, insuring that a nontrivial portion of their H 2k norms cannot concentrate on thinner layers near the boundary. This is akin to the the conditions in the Kolmogorov-Riesz theorem for compact sets in L p spaces.
This definition will enable us to achieve two key properties of our flow:

Regularization and good kernels
In what follows we outline the main steps developed in Section 2 of [14], and in which, for a given state (r, v) in H 2k , we construct regularized states, denoted by (r h , v h ), to our free boundary evolution, associated to a dyadic frequency scale 2 h , h ≥ 0. This relies on having good regularization operators associated to each dyadic frequency scale 2 h , h ≥ 0. We denote these regularization operators by h , with kernels K h . These are the same as in [14], and their exact definition can be found in there as well. A brief description on how one should envision these regularization operators is in order.
It is convenient to think of the domain t as partitioned in dyadic boundary layers, denoting by [ j] the layer at distance 2 −2 j away from the boundary. Within each boundary layer we need to understand which is the correct spatial regularization scale. The principal part of the second order elliptic differential operator associated to our system is the starting point. Given a dyadic frequency scale h, our regularizations will need to select frequencies ξ with the property that r ξ 2 2 2h , which would require kernels on the dual scale However, if we are too close to the boundary, i.e. r 2 −2h , then we run into trouble with the uncertainty principle, as we would have δx r . To remedy this issue we select the spatial scale r 2 −2h and the associated frequency scale 2 2h as cutoffs in this analysis. Then the way the regularization works is as follows: The regularized state is obtained by restricting the full regularization to the domain h := r h > 0 .
For completeness we state the result in [14], and refer the reader there for the proof:

3)
and iii) Higher regularity iv) Low frequency difference bound:

Embedding and interpolation theorems
In this section we state some embedding and interpolation results that will be used throughout.They have been proved in the last two authors' prior paper [14], to which the reader is referred to for the proofs. As a corollary of the above lemma we have embeddings into standard Sobolev spaces.

Lemma 2.4.
Assume that σ > 0 and σ ≤ j. Then we have In particular, by standard Sobolev embeddings, we also have Morrey type embeddings into C s spaces: where the equality can hold only if s is not an integer.
Next, we state the interpolation bounds.
Proposition 2.6. Let σ 0 , σ m ∈ R and 1 ≤ p 0 , p m ≤ ∞. Define and assume that Then for 0 < j < m we have Remark 2.7. One particular case of the above proposition which will be used later is when p 0 = p 1 = p 2 = 2, with the corresponding relation in between the exponents of the r σ j weights.
As the objective here is to interpolate between the L 2 type H m,σ norm and L ∞ bounds, we will need the following straightforward consequence of Proposition 2.6: Then for 0 < j < m we have We will also need the following two variations of Proposition 2.8: Then for 0 < j < m we have Then for 0 < j < m we have

The linearized equation
Consider a one-parameter family of solutions (r τ , v τ ) for the main system (1.16) such that (r τ , v τ )| τ =0 = (r, v), Then formally the functions d dτ (r τ , v τ ) τ =0 = (s, w), defined in the moving domain t , will solve the corresponding linearized equation. Precisely, a direct computation shows that, for (s, w) in t , the linearized equation can be written in the form where f and g i represent perturbative terms of the form with potentials V 1,2 and W 1,2 which are linear in ∂(r, v), with coefficients which are smooth functions of r and v. Importantly, we remark that for the above system we do not obtain or require any boundary conditions on the free boundary t . This is related to the fact that our one parameter family of solutions are not required to have the same domain, as it would be the case if one were working in Lagrangian coordinates.
For completeness, we also provide the explicit expressions for the potentials V 1,2 and W 1,2 , though this will not play any role in the sequel. We have where a 3 is a smooth function of (r, v), given by As for the other coefficients, the particular form of a 3 is not relevant, but we wrote it here for completeness.

Energy estimates and well-posedness
We now consider the well-posedness of the linearized problem (3.1) in the time dependent space H. For the purpose of this analysis, we will view H as a Hilbert space whose squared norm plays the role of the energy functional for the linearized equation, We will use this space for both the linearized equation and its adjoint. Our main result here is as follows: This in turn follows from a trivial pointwise bound on the corresponding potentials, We multiply (3.1a) by r Next, we add the two equations above, noting that the second and third terms on the LHS of the first equation combine with the second term on the LHS of the second equation to produce This yields We now integrate the above identity over t , using the formula (1.8) to produce a time derivative of the energy. For this, we need to write the terms on the left as perfect derivatives or material derivatives. When we do so the zero order coefficients do not cause any harm. We only need to be careful with the terms where a derivative falls on r 1−κ κ because this could potentially produce a term with the wrong weight (i.e., one less power of r ). However, this does not occur because we can solve for D t r in (1.16a): r, v)).
Using the above observations, we obtain We now compute the adjoint equation to (3.1) with respect to the duality relation defined by the H inner product determined by the norm (3.3). The terms f and g on the RHS of (3.1) are linear expressions in s and r w and and in s and w, respectively, with ∂(r, v) coefficients. Thus, the source terms in the adjoint equation have the same structure as the original equation. Let us write the LHS of (3.1) as With respect to the H inner product, the adjoint term corresponding to A i ∂ i is modulo terms that are linear expressions ins and rw and ins andw (with ∂(r, v) coefficients) in the first and second components, respectively, wheres andw are elements of the dual. Similarly, the adjoint term corresponding to B is Combining these expressions, we see that the bad term on the lower left corner of the second matrix inÃ i ∂ i cancels with the corresponding terms inB. Therefore, the adjoint problem is the same as the original one, modulo perturbative terms, and it therefore admits an energy estimate similar to the energy estimate (3.4) we have for the linearized Eq. (3.1).
In a standard fashion, the forward energy estimate for the linearized equation and the backward in time energy estimate for the adjoint linearized equation yield uniqueness, respectively existence of solutions for the linearized equation, as needed. This guarantees the well-posedness of the linearized problem.

Second order transition operators
An alternative approach the linearized equations is to rewrite the linearized Eq. (3.1) as a second order system which captures the wave-like part of the fluid associated with the propagation of sound. Applying D t to (3.1a) and using (3.1b) and vice-versa, and ignoring perturbative terms, we find Equations (3.6a) and (3.6b) are akin to wave equations in that the operatorsL 1 and L 2 satisfy elliptic estimates as proved in Section 5.2. More precisely, the operator L 2 is associated with the divergence part of w, and it satisfies elliptic estimates once it is combined with a matching curl operatorL 3 . Even though in this paper we do not use the operatorsL 1 andL 2 directly in connection to the corresponding wave equation, they nevertheless play an important role at two points in our proof. Because slightly different properties ofL 1 andL 2 are needed at these two points, we will take advantage of the fact that only their principal part is uniquely determined in order to make slightly different choices for L 1 andL 2 . Precisely, these operators will be needed as follows: I. In the proof of our energy estimates in Section 5.2, in order to establish the coercivity of our energy functionals. There we will need the coercivity ofL 1 andL 2 +L 3 , but we also want their coefficients to involve only r, ∇r and undifferentiated v.
II. In the constructive proof of existence of regular solutions, in our paradifferential style regularization procedure. There we use functional calculus for bothL 1 andL 2 +L 3 , so we need them to be both coercive and self-adjoint, but we no longer need to impose the previous restrictions on the coefficients.
The two sets of requirements are not exactly compatible, which is why two choices are needed. 10 We begin with the case (I), where we modifyL 1 andL 2 as follows: (3.8a) (3.8b) ToL 2 we associate an operatorL 3 of the form For case (II), we keep the first of the operators, setting L 1 =L 1 but make some lower order changes toL 2 andL 3 as follows: where F is the inverse of the matrix G, i.e., F = G −1 . It is not difficult to show that L 1 is a self-adjoint operator in L 2 (r 1−κ κ ) with respect to the inner product defined by the first component of the H norm in (2.1), and Similarly, both L 2 and L 3 are self-adjoint operators in L 2 (r 1 κ ) with respect to the inner product defined by the second component of the H norm in (2.1) and and similarly for L 3 . We further note that L 2 L 3 = L 3 L 2 = 0 and that the output of L 2 is a gradient, whereas L 3 w depends only on the curl of w.
As seen above, it is the operatorL 2 that naturally come out of the equations of motion rather than L 2 (recall that L 1 =L 1 ). Thus, we need to compare these operators; we also compare L 1 andL 1 for later reference. We have (3.12) We will establish coercive estimates for L 1 , L 2 , and L 3 (see Sects. 5 and 6), from which follows that the above domains can be characterized as

The uniqueness theorem
In this Section we establish Theorem 1.1. It will be a direct consequence of the more general Theorem 4.2 below, which establishes that a suitable functional that measures the difference between two solutions is propagated by the flow.
We consider two solutions (r 1 , v 1 ) and (r 2 , v 2 ) defined in the respective domains 1 t and 2 t . Put t := 1 t ∩ 2 t , t :=∂ t . If the boundaries of the domains 1 t and 2 t are sufficiently close, which will be the case of interest here, then t will have a Lipschitz boundary. Let D 1 t and D 2 t be the material derivatives associated with the domains 1 t and 2 t , respectively. In t define the averaged material derivative We note that the above averaged material derivative is not exactly advecting the domain t . Fortunately exact advection is not at all needed for what follows. See also Remark 4.3 below. To measure the difference between two solutions on the common domain t , we introduce the following distance functional: which is the same as in [14]. We could have used G to measure |v 1 − v 2 |, but the Euclidean metric suffices. This is not only because both metrics are comparable but also because we will not use (4.1) directly in conjunction with the equations.
We will, however, use G further below when we introduce another functional for which the structure of the equations will be relevant.
We observe the following Lemma, which has been proved in [14]: One can view the integral in (4.2) as a measure of the distance between the two boundaries, with the same scaling as D H . We now state our main estimate for differences of solutions: We remark that which implies our uniqueness result. The remaining of this Section is dedicated to proving Theorem 4.2.

A degenerate energy functional
We will not work directly with the functional D H because it is non-degenerate, so we cannot take full advantage of integration by parts when we compute its time derivative. We thus consider the modified difference functional where a 21   This is why we refer to this difference functional as degenerate, and is also the reason why we are able to use the averaged material derivative D t to propagate bounds forD H in time even though it does not exactly advect t .
From [14] we also borrow the equivalence property of the two distance functionals defined above: (4.5)

The energy estimate
To prove the Theorem 4.2 it remains to track the time evolution of the degenerate distance functionalD H . This is the content of the next result, which immediately implies Theorem 4.2 after an application of Gronwall's inequality.

Proposition 4.5. We have
Proof. The difference of the two solutions to (1.16) in the common domain t satisfies and Above, G i and a 1i correspond to G and a 1 for the solutions (r i , v i ), i = 1, 2. We observe that the difference of material derivatives can be written as To compute the time derivative of the degenerate distance we use the standard formula for differentiation in a moving domain t , (4.9) where v is in our case the average velocity. Classically this holds under the assumption that the domain t is advected by D t . But in our case we replace this assumption with the alternative condition that f vanishes on the boundary of t . Using this formula we obtain d dtD where the integrals I i , i = 1, 6, represent contributions as follows: (a) I 1 represents the contribution where the averaged material derivative falls on a or b, Here the derivatives of a and b are homogeneous of order −1, respectively 0. We get Gronwall terms when they get coupled with factors of r 1 + r 2 or r 1 − r 2 from the material derivatives. We discuss I 1 later.
(b) I 2 gathers the contributions where the averaged material derivative is applied to (a 21 + a 22 ) −1 and to G i j mid . These expressions are smooth functions of (r 1 , v 1 ), (r 2 , v 2 ), and thus their derivatives are bounded by B 1 + B 2 , (c) I 3 represents the main contribution of the averaged material derivative that falls onto (r 1 − r 2 ) respectively on v 1 − v 2 which consists of the first and second terms on the RHS in (4.7), and the second term in (4.8): This term will need further discussion. (d) In I 4 we place the contribution of the forth term on the RHS of (4.7): This term will be discussed later.
(e) I 5 is given by the third and fifth terms on the RHS in (4.7) and the third term on the RHS from (4.8) All of the terms in I 5 are straightforward Gronwall terms.
(f) I 6 contains the terms where D t falls on (r 1 + r 2 ) 1−κ κ : We will analyze I 6 later. It remains to take a closer look at the integrals I 1 , I 3 , I 4 , and I 6 . We consider them in succession.
The bound for I 1 . Here we write The first two terms have size O(B(r 1 +r 2 )) so their contribution is a Gronwall term. We are left with the contribution of the last terms, which yields the expressions For the second integral in both expressions, we bound |ṽ 1 −ṽ 2 | |r 1 −r 2 |+|v 1 −v 2 | and estimate their part by and which is discussed later.
We are left with the first integrals in I a 1 and I b 1 , which we record as These integrals are also discussed later. The bound for I 3 . For I 3 , we seek to capture the same cancellation that it is seen in the analysis of the linearized equation. We look at the last term in I 3 , use b = a(r 1 + r 2 ), and integrate by parts; if the derivatives falls on G then this is a straightforward Gronwall term. We are left with three contributions, two of which we pair with the first two terms in I 3 . We obtain For the first integral, I 1 3 we expand the differenceṽ i 1 −ṽ i 2 , seen as a function of v 1 and v 2 , in a Taylor series around the center (v 1 + v 2 )/2. We have where we recognize the matrix on the right as being the main part in G. Thus, we can write where the quadratic v 1 − v 2 terms cancel because we are expanding around the middle, and we used (3.2) to get the second term on the right. The contributions of all of the terms in the last expansion are Gronwall terms. For the second integral I 2 3 we have a simpler expansion where all contributions qualify again as Gronwall terms. Finally, the last integral, J c 1 , is estimated below.

The bound for I 4 . After an integration by parts we have
Writing The bound for I 6 . We use where the first two terms are bounded by (B 1 + B 2 )(r 1 + r 2 ) and yield Gronwall contributions. Then we write To summarize the outcome of our analysis so far, we have proved that It remains to estimate J 2 , J a 1 , J b 1 , J c 1 , and J d 1 , which we write here again for convenience: The integral J 2 is the same as in [14] and can be estimated accordingly, using interpolation inequalities; see Lemma 4.4 in [14].
The bound for the integrals J a 1 , J b 1 , J c 1 and J d 1 matches estimates for similar integrals in [14]. Precisely, the integrals J a 1 and J b 2 are estimated as the integrals called J b 1 and J c 1 in [14], respectively, see Lemmas 4.6-4.13 in [14]. The integral J c 1 is estimated as the integral J d 1 in [14], see Lemmas 4.6-4.13 in [14]. The integral J d 1 is estimated as the integral J a 1 in [14], see Lemmas 4.6-4.13 in [14].
We caution the reader that the arguments in [14] are not straightforward, and involve peeling off a carefully chosen boundary layer, with separate estimates inside the boundary layer and outside it. The only difference in the present paper is the presence of additional weights in the integrals, which are smooth functions of r 1 , r 2 , v 1 , v 2 . For instance, the differenceṽ 1 −ṽ 2 can be expanded as where r, v stand for r 1 , r 2 , v 1 , v 2 and f, g are smooth. The contribution of the first term admits a straightforward Gronwall bound, and the contribution of the second term is akin to the corresponding integral in [14] but with the added smooth weight. The point is that every time we integrate by parts and the derivative falls on the smooth weight, the corresponding contribution is a straightforward Gronwall term. Hence such smooth weights make no difference if added in the arguments in [14].

Energy estimates for solutions
Our goal in this section is to establish uniform control over the H 2k norm of the solutions (r, v) to (1.16) with growth given by the norms A and B. For this, we will use appropriate energy functionals E 2k = E 2k (r, v) constructed out of vector fields naturally associated with the evolution. Our functionals will be associated with the wave and transport parts of the system, which will be considered at matched regularity.
The vector fields we will consider are: • The material derivative D t , which has order 1/2.
• All regular derivatives ∂, which have order 1.
• Multiplication by r , which has order −1.
The wave part of the energy will be associated mainly with D t , whereas the transport part will be associated with all of the above vector fields.

Good variables and the energy functional
Heuristically, higher order energy functionals should be obtained by applying an appropriate number of vector fields to the equation, and then verifying that the output solves the linearized equation modulo perturbative terms. In the absence of the free boundary, there are two equally good choices, (i) to spatially differentiate the equation, using the ∂ j vector fields, or (ii) to differentiate the equation in time, using the ∂ t vector field.
However, in the free boundary setting, both of the above choices have issues, as neither ∂ j nor ∂ t are adapted to the boundary. For ∂ t we do have a seemingly better choice, namely to replace it by the material derivative D t . However, this has the downside that it does not arise from a symmetry of the equations, and consequently the expressions (D 2k t r, D 2k t v) are not good approximate solutions to the linearized equations. To address this matter, we add suitable corrections to these expressions, We also define which we think of as the vorticity counterpart to (s 2k , w 2k ). These we will think of as solving approximate transport equations; using (1.22) we find where h 2k is given by We introduce the wave energy , and the total energy

Energy coercivity
Our goal in this section is to show that the energy (5.7) measures the H 2k size of (r, v). To do so, we would like to consider the energy as a functional of (r, v) defined at a fixed time. This can be done by using Eq. (1.16) to algebraically solve for spatial derivatives of (r, v). Let (r, v) be smooth functions in . Assume that r is positive in and uniformly non-degenerate on the . Then Proof. We begin with the part. We consider the wave part of the energy and the corresponding expressions for (s 2k , w 2k ). Using use Eq. (1.16a) and (1.16b) to successively solve for D t (r, v), we obtain that each (s 2k , w 2k ) is a linear combination of multilinear expressions in r and ∂v (with order zero coefficients). We will use our bookkeeping scheme of Section 1.3 to understand the expressions for (s 2k , w 2k ). It is useful to record here the order and structure of the linear-in-derivatives top order terms obtained by using the equations to inductively compute D 2k t (r, v) and D 2k−1 t (r, v), which involve 2k and 2k − 1 derivatives, respectively: We begin with the expressions of highest order (see Section 1.3), thus we first focus on the multilinear expressions that come from ignoring the last term on LHS (1.16a) and also where no derivative lands on G, a 1 , and a 2 . We also consider first the case when every time we commute D t with ∂, the derivative lands on v i and not on r via v 0 .
In this case, the corresponding multilinear expressions for (s 2k , w 2k ) have the following properties: • They have order k − 1 and k − 1 2 , respectively. • They have exactly 2k derivatives.
• They contain at most k + 1 and k factors of r , respectively.
Thus, a multilinear expression for s 2k in this case has the form where n j , m l ≥ 1, and subject to n j + m l = 2k, a + J + L/2 = k + 1, (5.10) and when J = 0 or L = 0 the corresponding product is omitted. We claim that it is possible to choose b j and c l such that This follows from observing that since a ≤ k. Equality holds only if a = k, J = 1 and L = 0 (i.e., for the leading linear case). This shows that it is possible to make such a choice of b j and c l , which allows us to use our interpolation theorems Observe that the numerators in 1/ p j and 1/q l correspond to the orders of the expressions being estimated and they add up to k − 1 (as needed).
In addition to the principal part discussed above, we also obtain lower order terms in our expression for s 2k . There are three sources of such terms: i) The terms from the commutator [D t , ∂] where derivatives apply to r via v 0 .
This corresponds to the second term in the formal expansion whose order is easily seen to be 1/2 lower. ii) If any derivatives are applied to either r or v via a 0 , a 1 , a 2 or G, this increases the order of the resulting expression by 0, respectively 1/2, compared to the full order of the derivative which is 1. iii) Contributions arising from the last term in (1.15), whose order is, to start with, 1/2 lower than the rest of the terms in the (1.15).
The contributions of all such terms to s 2k have lower order. More precisely, they contain expressions of the form (5.9) but with (5.10) replaced by n j + m l = 2k, where j > 0, and which have lower order k − 1 − j 2 . All such lower order terms can be estimated in a similar fashion, but using lower Sobolev norms for (r, v).
We continue with the part. Applying D t to (5.2a) and (5.2b) and using definitions (5.1) we find the following recurrence relations whereL 1 andL 2 have been defined in Section 3. The next Lemma characterizes the error terms on the RHS of (5.12), and the lemma that follows gives a quantitative relation between the 2 j and 2 j − 2 quantities.

Lemma 5.2.
For j ≥ 2, the terms F 2 j and G 2 j in (5.12) are linear combinations of multilinear expressions in r and ∂v with 2 j derivatives and of order at most j − 1 and j − 1 2 , respectively. Moreover, they are either (i) non-endpoint, by which we mean multilinear expressions of order j − 1 and j − 1 2 , respectively, containing at most j + 1 and j factors of r , respectively, and whose products contain at least two factors of ∂ ≥2 r or ∂ ≥1 v, or (ii) they have order strictly less than j − 1 and j − 1 2 , respectively, and contain at most j + 2 and j + 1 factors of r , respectively.
Proof. We begin with j ≥ 3. We will analyze In order to keep track of terms according to the statement of the Lemma, we observe that s 2 j has order j − 1. We will make successive use of the commutator We begin with the first term on RHS (5.13). From (1.16), we compute.
and we will consider each term on RHS (5.14) separately. The terms ∂ i r G ml and ∂ i (ra 1 v l ∂ l )r have order at most zero, thus has order at most j − 3/2 and belongs to F 2 j . Next, For the first term on the RHS, we have The second and third terms have order ≤ −1 and −1/2, respectively, thus they belong to F 2 j after differentiation by D 2 j−2 t . For the first term, we have The first term satisfies the non-endpoint property while the second has order −1/2, thus both terms belong to F 2 j after differentiation by D 2 j−2 t . Next, The second term has order −1/2 so it belongs to F 2 j upon differentiation by D 2 j−2 t . The first term has order zero, thus producing a top order (i.e., j − 1) term when differentiated by D 2 j−2 t . Nevertheless, it has two ∂ ≥1 v terms so it satisfies the non-endpoint property and hence it also belongs to F 2 j .
We now turn to the other commutator in (5.14): The terms on the RHS have orders ≤ −1/2, −3/2, −1, −1, −1/2, −1, respectively, so they all belong to F 2 j upon differentiation by D 2 j−2 t . For the last two terms on RHS (5.14), we see that they have orders −1 and −1/2, thus also belong to F 2 j after differentiation by D 2 j−2 t . Therefore, writing ≈ for equality modulo terms that belong to F 2 j , (5.14) becomes In the second sum, we can further write The term D 2 j−2− t r D t (G ml ∂ l a 2 ∂ m r ) has order at most j −3/2 and can be absorbed into F 2 j . For the first term, if D t hits G ml a 2 we again obtain a term of order strictly less than j − 1 that is part of F 2 j . Finally, for the term we use that ≤ 2 j − 2 − 1 and (1.16a) to write The first term contains a ∂ ≥1 v and a ∂ ≥2 r so it belongs to F 2 j , whereas the second term has order at most j − 3/2 so it belongs to F 2 j as well. Hence, we have that We now compute the commutator on the second term on the RHS: The first and second terms on the RHS of the second equality have orders ≤ 1/2 and 1, respectively, so they produce terms of order at most j − 3/2 when hit by r D 2 j−3 t and thus can be discarded. Continuing All the terms inside the second parenthesis have orders at most 1 (thus giving order at most j − 3/2 when hit by r D The first term on RHS (5.16) is treated with similar ideas. We notice that the top order term in that expression is The first term belongs toL 1 s s j−2 and the second one to F 2 j . The case j = 2 is done separately (since the definition of s 2 is different, recall (5.1)), but it follows essentially the same steps as above. Finally, the proof for G 2 j is done with the same type of calculations employed above and we omit it for the sake of brevity.
To continue our analysis, we need some coercivity estimates for theL 1 , respec-tivelyL 2 +L 3 .

Lemma 5.3. Assume that A is small. Then
Here we remark that the lower order terms on the right play no role in the proof, and can be omitted if (s, w) are assumed to have small support (by Poincare's inequality), or if we use homogeneous norms on the left.
As a consequence of the second estimate above, we have .
In Section 6 will also need the following straightforward alternative form of the above result: Corollary 5.5. Assume that B is small. Then the same result as in Lemma 5.3 holds for the operators L 1 , respectively L 2 + L 3 .
Here the smallness condition on B allows us to treat the differencesL 1 − L 1 , Proof. We start with two simple observations. First of all, using a partition of unity one can localize the estimates to a small ball. We will assume this is done, and further we will consider the interesting case where this ball is around a boundary point x 0 ; the analysis is standard elliptic otherwise. We can assume that at x 0 on the boundary we have ∇r (x 0 ) = e n so that in our small ball we have |∇r − e n | A 1. (5.18) Secondly, the smallness condition on A guarantees that the coefficients G and a 2 have a small variation in a small ball, and we can freeze these coefficients modulo perturbative errors. Hence, we will simply freeze them, and assume that a 2 and G are constant. Then a 2 only plays a multiplicative role, and will be set to 1 for the rest of the argument.
A preliminary step in the proof is to observe that we have the weaker bounds These bounds can be proved in a standard elliptic fashion by integration by parts, e.g. in the case of the first bound one simply starts with the integral representing L 1 s 2 H 0, 1 2κ and exchange derivatives between the two factors. The details are left for the reader.
In view of the above bounds, it suffices to show that which suffices, by the Cauchy-Schwarz inequality. Now we consider (5.17b). As discussed above, we set a 2 = 1 and assume G is a constant matrix. We recall thatL 2 has the form whileL 3 is given by Then a direct computation shows that We will take advantage of the covariant nature of this operator in order to simplify it. Interpreting G as a dual metric and w as a one form, we see that the above operator viewed as a map from one forms to one forms is invariant with respect to linear changes of coordinates. Here we are interested in changes of coordinates which preserve the surfaces x n = const. But even with this limitation, it is possible to choose a linear change of coordinates, namely the semigeodesic coordinates relative to the surface x n = 0, y = Ax + bx n , y n = x n so that the metric G becomes a multiple of the identity. Then the estimate (5.19b) reduces to its euclidean counterpart, which is discussed in detail in [14] in the corresponding nonrelativistic context.
To finish the proof of Theorem 5.1, we will establish 22) where ε > 0 is sufficiently small. We are using ε here to include two types of small error terms: (a) the terms that we estimate using O(A) as well as (b) the terms that have an extra factor of r and for which we can use smallness of r near the boundary; the latter type arise from the last term of (1.16a). Concatenating these estimates we then obtain the conclusion of the theorem.
To prove (5.22), we first consider (F 2 j , G 2 j ) H 2k−2 j . Using our interpolation inequalities, the non-endpoint property, and the structure of (F 2 j , G 2 j ) described in in Lemma 5.2, we obtain It remains to handle the term (s 2 j , w 2 j ) H 2k−2 j . For j = k the desired estimate is a direct consequence of Lemma 5.3.
We move to treat the case 2 ≤ j < k. The idea is to apply Lemma 5.3 with s 2 j−2 and w 2 j−2 replaced by suitable weighted derivatives of themselves. More precisely, we set s: Applying L to (5.12a), we obtain The term L F 2 j can again be dealt with using Lemma 5.2, as above. Thus we focus on the commutator. To analyze it, we consider induction on a, starting at a = 0, and observe the following: • All terms where at least one r factor gets differentiated twice are non-endpoint terms and can be estimated by interpolation. • The terms where two r factors are differentiated are handled by the induction on a. • Terms where only one r gets differentiated are also handled by induction on a unless a = 0.
Therefore, all terms in the commutator where a > 0 are perturbative terms. We now focus on the case a = 0. Consider a frame (x , x n ) in Minkowski space that is adapted to a point near the boundary in the sense that Then, all terms in the commutator with tangential derivatives only are error terms. For terms involving ∂ n , we find where primed indices run from 1 to n − 1. The last two terms on the RHS can be treated by yet another induction, this time over b. The first term on the RHS can be combined back withL 1 The operatorL b 1 has a similar structure toL 1 , and an inspection in the proof of Lemma 5.3 shows that the corresponding coercive estimate for s holds withL b 1 in place ofL 1 .
The above argument works for j ≥ 2 in that (5.12a) is valid only for j ≥ 2. However, a minor change in the above using the definition s 2 yields the result also for j = 1. This takes care of the s terms in (5.22); the proof for the w terms is similar.

Energy estimates
In this Section we establish Theorem 5.6. The energy functional E 2k defined in (5.7) satisfies the following estimate: Proof. In view of Eqs. (5.2a)-(5.2b) and (5.5), the the energy estimates for the linearized equation in Section 3, and estimates for transport equations, it suffices to show that the terms f 2k , g 2k and h 2k , given by (5.3a), (5.3b), and (5.6), respectively, are perturbative, i.e., they satisfy the estimate To prove this bound we need to understand the structure of ( f 2k , g 2k ), respectively h 2k .
Lemma 5.7. Let k ≥ 1. Then source terms f 2k and g 2k in the linearized Eq. (5.2) for (s 2k , w 2k ), given by (5.3a)-(5.3b) are multilinear expressions in (r, ∇v), with coefficients which are smooth functions of (r, v), which have order ≤ k − 1 2 , respectively ≤ k, with exactly 2k + 1 derivatives, and which are not endpoint, in the sense that there is no single factor in f 2k , respectively g 2k which has order larger that k − 1, respectively k − 1 2 . Similarly, the source term h 2k in the vorticity transport Eq. (5.5), given by (5.6), has the same properties as g 2k above.
Once the lemma is proved, arguing similarly to Section 5.2, we see that this suffices to apply our interpolation results in Propositions 2.6, 2.9 and 2.10 and obtain the desired bound. Here we remark that a scaling analysis shows that in the interpolation estimates we need to use at most one B control norm, with equality exactly in the case of terms of highest order. One should also compare with the situation in the similar computation in [14], where no lower order terms appear. Hence, the poof of the theorem is concluded once we prove the above lemma.
Proof of Lemma 5.7. Consider first f 2k . The fact that all terms in f k have order at most k − 1 2 is obvious. The non-endpoint property can be understood as asking that there are no derivatives of order 2k + 1, and that, in addition, for the terms of maximum order, they have at least two factors of the form ∂ 2+ r or ∂ 1+ v. Notably, this excludes any terms of the form A similar reasoning applies for g 2k and h 2k , where the forbidden terms are those with a factor with 2k + 1 derivatives, as well as those of maximum order of the form We start with a simple observation, which is that, if in (5.3a) or (5.3b), any derivative falls on a coefficient such as G, a 0 , a 1 , or a 2 , then we obtain lower order terms which automatically satisfy the above criteria. Thus, for the purpose of this Lemma we can treat these coefficients as constants.
A second observation is that there are no factors with 2k +1 derivatives in either s 2k or w 2k , due to the commutator structures present in both (5.3a) or (5.3b). This directly allows us to discard all lower order terms, and in particular those containing a 1 and a 3 . By the same token we can set a 0 = 1 and r = 1.
Given the above observations, it suffices to consider the reduced expressions first. When commuting ∂ and D 2k t , this produces at least one ∂v, so [r ∂ i , D 2k t ]v j is not an endpoint term. Similarly, D t (∂ i r ) has order 1/2 so the second expression is also not an endpoint term.
We now investigate g reduced 2k .
Neither of the first two terms is perturbative, but we have a leading order cancellation between them, based on the relations The contribution of the second term is lower order and thus perturbative. The contribution of the first term is combined with the second term in (5.24) to obtain a commutator structure which yields only balanced terms. The third term in (5.24) is also balanced due to the commutator structure, while the last term has a direct good factorization.
We next move to h 2k . From (5.6) we see that we are commuting r a ∂ b with either D t or ∂v, so we always obtain ∂v factors that give non-endpoint terms. The only possible exception is when all derivatives in the commutator with D t are applied to the r term in v 0 . But this yields a lower order term.

Construction of regular solutions
In this section we provide the first step in our proof of local well-posedness, namely, here we present a constructive proof of regular solutions. The rough solutions are obtained in the last section as unique limits of regular solutions.
Given an initial data (r ,v) with regularity where k is assumed to be sufficiently large, we will construct a local in time solution, bounded in H 2k , with a lifespan depending on the H 2k size of the data.
The first property will ensure a uniform energy bound for our sequence. The second property will guarantee that in the limit we obtain an exact solution. There we use a weaker topology, where the exact choice of norms is not so important (e.g. C 2 ).
Having such a sequence of approximate solutions, it is straightforward to produce, as the limit on a subsequence, an exact solution (r, v) on a short time interval which stays bounded in the above topology. The key point is the construction of the above sequence. It suffices to carry out a single step: whereG,å 1 , andå 2 are G, a 1 , and a 2 evaluated at (r ,v).
The strategy for the proof of the theorem is the same as in the last two authors' previous paper [14], by splitting the time step into three: • Regularization, • Transport, • Euler's method, where the role of the first two steps is to improve the error estimate in the third step. The regularization step is summarized in the next Proposition: Proof. We repeat the construction in [14]. There are only a few minor differences, namely • The self-adjoint operators L 1 , L 2 and L 3 there are replaced by their counterparts in this paper, i.e., (3.7a), (3.10), and (3.9) (recall that L 1 =L 1 and L 3 =L 3 ). • Using (3.12), relations similar to (5.12) continue to hold for the self-adjoint operators. Thus, the approximate relations between (s 2k , w 2k ) and (s − 2k , w − 2k ) in Section 6 of [14] also hold here.
• The elliptic estimates of Lemma 5.3 hold for L 1 and L 2 , L 3 , with essentially the same proof.
Aside from the above minor differences, the most important observation in invoking the proof given in [14] is that the counterpart of Lemma 6.3 in [14] still holds with a minor change. For convenience we state here its counterpart (below, Ds 2k and Dw 2k are the differentials of s 2k and w 2k as functions of r and v): We have the algebraic relations where the error terms (F 2k ,G 2k ) are linear in (r −ř ,v −v), Their coefficients are multilinear differential expressions in (ř ,v), have order at most k − 1, respectively k − 1 2 , and whose monomials fall into one of the following two classes: i) Have maximal order but contain at least one factor with order > 0, i.e. ∂ 2+ř or ∂ 1+v , or ii) Have order strictly below maximum.
By comparison, the similar relations in Lemma 6.3 in [14] are homogeneous, so only terms of type (i) arise in the error terms. Here our equations are no longer homogeneous, and lower order terms do appear. In particular, we note that all the contributions coming from the last term in the first equation (1.16a) belong to the class (ii) above. This is correlated with and motivates the fact that this term was neglected in our definition of the operator L 1 .
With these observations in mind, the proof given in [14] applies directly here. We now use Proposition 6.2 in order to prove Theorem 6.1. Proof of Theorem 6.1. For the transport step, we defině where, in agreement with the our definition of the material derivative, we iterate the coordinates by flowing with v i /v 0 , and not simply v i .
Then we carry out the Euler step, and define (r 0 , v 0 ) by To show that (ř ,v) have the properties in the Theorem, the argument is completely identical to the one in [14].

Rough solutions and continuous dependence
The last task of the current work is to construct rough solutions as limits of smooth solutions, and conclude the proof of Theorem 1.2. Fortunately, the arguments in the preceding paper [14] by the last two authors for the similar part of the results apply word for word. This is despite the fact there are several differences between the two problems that play a role on how the energy estimates are obtained, as well as on how uniqueness is proved. However, the functional framework developed in [14] and also implemented here does not see these differences. Furthermore, the proof of the similar result in [14] only uses (i) the regularization procedure in Section 2, (ii) the difference bounds of Theorem 1.1, and (iii) the energy estimates of Theorem 1.3, without any reference to their proof.
Thus, in our current result we rely on the same succession of steps as in the nonrelativistic companion work of the last two authors [14], which we briefly outline here for the reader. These steps are 1. Regularization of the initial data. We regularize the initial data; this is achieved by considering a family of dyadic regularizations of the initial data as described in Section 2. These data generate corresponding smooth solutions by Theorem 1.2. For these smooth solutions we control on the one hand higher Sobolev norms H 2k+2 j using our energy estimates in Theorem 1.3, and on the other hand the L 2 -type distance between consecutive solutions, which is at the level of the H norms, by Theorem 1.1.

Uniform bounds for the regularized solutions.
To prove these bounds we use a bootstrap argument on our control norm B, where B is time dependent. The need for an argument of this kind is obvious. Once we have the regularized data sets (r h ,v h ), we also have the corresponding smooth solutions (r h , v h ) generated by the smooth data (r h ,v h ). A-priori these solutions exist on a time interval that depends on h. Instead, we would like to have a lifespan bound which is independent of h. This step requires closing the bootstrap argument via the energy estimates already obtained in Section 5. 3. Convergence of the regularized solutions. We obtain the convergence of the regular solutions (r h , v h ) to the rough solution (r, v) by combining the high and the low regularity bounds directly. This yields rapid convergence in all H 2k spaces below the desired threshold, i.e. for k < k. Here we rely primarily on results in Section 4, namely Theorem 1.1. 4. Strong convergence. Here we prove the convergence of the smooth solutions to the rough limit in the strong topology H 2k . To gain strong convergence in H 2k we use frequency envelopes to more accurately control both the low and the high Sobolev norms above. This allows us to bound differences in the strong H 2k topology. A similar argument yields continuous dependence of the solutions in terms of the initial data, also in the strong topology. For more details we refer the reader to [14].
Acknowledgements. The first author was partially supported by NSF grant DMS-2107701, a Sloan Research Fellowship provided by the Alfred P. Sloan foundation, a Discovery grant administered by Vanderbilt University, and a Dean's Faculty Fellowship. The second author was supported by a Luce Assistant Professorship, by the Sloan Foundation, and by an NSF CAREER grant DMS-1845037. The last author was supported by the NSF grant DMS-1800294 as well as by a Simons Investigator grant from the Simons Foundation.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/ licenses/by/4.0/.
Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.