Causality of the Einstein-Israel-Stewart Theory with Bulk Viscosity

TL;DR

This paper proves that Einstein's equations coupled with Israel-Stewart-type fluid dynamics with bulk viscosity are causal and well-posed in the full nonlinear regime, enabling more accurate neutron star merger simulations.

Contribution

It demonstrates the causality and well-posedness of Einstein-Israel-Stewart equations with bulk viscosity without symmetry assumptions.

Findings

Equations are causal in the full nonlinear regime.

The system can be written as a symmetric hyperbolic first-order system.

Results apply to arbitrary equations of state.

Abstract

We prove that Einstein's equations coupled to equations of Israel-Stewart-type, describing the dynamics of a relativistic fluid with bulk viscosity and nonzero baryon charge (without shear viscosity or baryon diffusion) dynamically coupled to gravity, are causal in the full nonlinear regime. We also show that these equations can be written as a first-order symmetric hyperbolic system, implying local existence and uniqueness of solutions to the equations of motion. We use an arbitrary equation of state and do not make any simplifying symmetry or near-equilibrium assumption, requiring only physically natural conditions on the fields. These results pave the way for the inclusion of bulk viscosity effects in simulations of gravitational-wave signals coming from neutron star mergers.

Equations

$$\left[\frac{\zeta }{{\tau }_{\Pi }}+n{\left(\frac{\partial P}{\partial n}\right)}_{\epsilon }\right]\frac{1}{\epsilon +P+\Pi }\le 1-{\left(\frac{\partial P}{\partial \epsilon }\right)}_n,$$

$$u^{\alpha }{\nabla }_{\alpha }\epsilon +\left(\epsilon +P+\Pi \right){\nabla }_{\alpha }{u}^{\alpha }$$

$$\left(\epsilon +P+\Pi \right)u^{\beta }{\nabla }_{\beta }{u}_{\alpha }+{\alpha }_1{\Delta }_{\alpha }^{\beta }{\nabla }_{\beta }\epsilon$$

$$+{\alpha }_2{\Delta }_{\alpha }^{\beta }{\nabla }_{\beta }n+{\Delta }_{\alpha }^{\beta }{\nabla }_{\beta }\Pi$$

$${\tau }_{\Pi }{u}^{\alpha }{\nabla }_{\alpha }\Pi +\Pi +\lambda {\Pi }^2+\zeta {\nabla }_{\alpha }{u}^{\alpha }=0,$$

$$u^{\mu }{\nabla }_{\mu }n+n{\nabla }_{\mu }{u}^{\mu }=0.$$

$$R_{\mu \nu }-\frac{1}{2}R{g}_{\mu \nu }+\Lambda g_{\mu \nu }=8\pi G{T}_{\mu \nu },$$

$$\left(\begin{array}{ccccc}\frac{u^0}{(\epsilon +P+\Pi )}& -\frac{u_i}{u_0}& 0& 0& \\ -\frac{u_j}{u_0}& A_{ij}{u}^0\frac{\left(\epsilon +P+\Pi \right)}{{\left(\frac{\partial P}{\partial \epsilon }\right)}_n}& -\frac{u_j}{u_0}\frac{{\left(\frac{\partial P}{\partial n}\right)}_{\epsilon }}{{\left(\frac{\partial P}{\partial \epsilon }\right)}_n}& -\frac{u_j}{u_0}\frac{1}{{\left(\frac{\partial P}{\partial \epsilon }\right)}_n}& \\ 0& -\frac{u_i}{u_0}\frac{{\left(\frac{\partial P}{\partial n}\right)}_{\epsilon }}{{\left(\frac{\partial P}{\partial \epsilon }\right)}_n}& \frac{u^0}{n}\frac{{\left(\frac{\partial P}{\partial n}\right)}_{\epsilon }}{{\left(\frac{\partial P}{\partial \epsilon }\right)}_n}& 0& \\ 0& -\frac{u_i}{u_0}\frac{1}{{\left(\frac{\partial P}{\partial \epsilon }\right)}_n}& 0& \frac{{\tau }_{\Pi }}{\zeta {\left(\frac{\partial P}{\partial \epsilon }\right)}_n}{u}^0& \end{array}\right)$$

$$\left(\begin{array}{cccc}\frac{u^k}{\left(\epsilon +P+\Pi \right)}& {\delta }_i^k& 0& 0\\ {\delta }_j^k& A_{ij}{u}^k\frac{\left(\epsilon +P+\Pi \right)}{{\left(\frac{\partial P}{\partial \epsilon }\right)}_n}& {\delta }_j^k\frac{{\left(\frac{\partial P}{\partial n}\right)}_{\epsilon }}{{\left(\frac{\partial P}{\partial \epsilon }\right)}_n}& \frac{{\delta }_j^k}{{\left(\frac{\partial P}{\partial \epsilon }\right)}_n}\\ 0& {\delta }_i^k\frac{{\left(\frac{\partial P}{\partial n}\right)}_{\epsilon }}{{\left(\frac{\partial P}{\partial \epsilon }\right)}_n}& \frac{u^k}{n}\frac{{\left(\frac{\partial P}{\partial n}\right)}_{\epsilon }}{{\left(\frac{\partial P}{\partial \epsilon }\right)}_n}& 0\\ 0& \frac{{\delta }_i^k}{{\left(\frac{\partial P}{\partial \epsilon }\right)}_n}& 0& \frac{{\tau }_{\Pi }}{\zeta {\left(\frac{\partial P}{\partial \epsilon }\right)}_n}{u}^k\end{array}\right),$$

$${\tau }_{\Pi }{u}^{\alpha }{\nabla }_{\alpha }\Pi +\zeta {\nabla }_{\alpha }{u}^{\alpha }+f(\epsilon ,n,\Pi ,{\tau }_{\Pi },\zeta )=0,$$

$$\left(\begin{array}{cccc}{u}^{\mu }& \left(\epsilon +P+\Pi \right){\delta }_{\nu }^{\mu }& 0& 0\\ {\alpha }_1{\Delta }^{\alpha \mu }& 1_{4\times 4}\left(\epsilon +P+\Pi \right)u^{\mu }& {\alpha }_2{\Delta }^{\alpha \mu }& {\Delta }^{\alpha \mu }\\ 0& n{\delta }_{\nu }^{\mu }& u^{\mu }& 0\\ 0& \zeta {\delta }_{\nu }^{\mu }& 0& {\tau }_{\Pi }{u}^{\mu }\end{array}\right),$$

$$\left(\begin{array}{cccc}0& 0_{1\times 4}& 0& 0\\ 0_{4\times 1}& 0_{4\times 4}& 0_{4\times 1}& 0_{4\times 1}\\ 0& 0_{1\times 4}& 0& 0\\ 0& 0_{1\times 4}& 0& 1+\lambda \Pi \end{array}\right).$$

$$\mathcal{P}(\xi )$$

$$\displaystyle\quad-\Big{(}\frac{\zeta}{\tau_{\Pi}(\varepsilon+P+\Pi)}+\alpha_{1}$$

$$\displaystyle\quad+\frac{\alpha_{2}\,n}{\varepsilon+P+\Pi}\Big{)}\Delta^{\mu\lambda}\xi_{\mu}\xi_{\lambda}\Big{]}.$$

$$-\frac{1}{2}{g}^{\mu \nu }{\partial }_{\mu }{\partial }_{\nu }{g}_{\alpha \beta }$$

$$+\Lambda g_{\alpha \beta }+{\mathcal{Z}}_{\alpha \beta },$$

$$\left(\begin{array}{ccccc}{1}_{10}& 0_{10}& 0_{10}& 0_{10}& 0_{10}\\ 0_{10}& g^{11}{1}_{10}& g^{12}{1}_{10}& g^{13}{1}_{10}& 0_{10}\\ 0_{10}& g^{12}{1}_{10}& g^{22}{1}_{10}& g^{23}{1}_{10}& 0_{10}\\ 0_{10}& g^{13}{1}_{10}& g^{23}{1}_{10}& g^{33}{1}_{10}& 0_{10}\\ 0_{10}& 0_{10}& 0_{10}& 0_{10}& -g^{00}{1}_{10}\end{array}\right),$$

$$\left(\begin{array}{ccccc}{0}_{10}& 0_{10}& 0_{10}& 0_{10}& 0_{10}\\ 0_{10}& 0_{10}& 0_{10}& 0_{10}& g^{k1}{1}_{10}\\ 0_{10}& 0_{10}& 0_{10}& 0_{10}& g^{k2}{1}_{10}\\ 0_{10}& 0_{10}& 0_{10}& 0_{10}& g^{k3}{1}_{10}\\ 0_{10}& g^{k1}{1}_{10\times 10}& g^{k2}{1}_{10}& g^{k3}{1}_{10}& 2g^{k0}{1}_{10}\end{array}\right),$$

$${\mathfrak{A}}^{\mu }=\left(\begin{array}{cc}{\mathcal{A}}^{\mu }& 0\\ 0& {\tilde{\mathcal{A}}}^{\mu }\end{array}\right).$$

$$\det {\mathcal{A}}^0$$

$$\displaystyle\quad-Kv^{2})-v^{2}\Big{]},$$

$$(\epsilon +P+\Pi )\frac{{\tau }_{\Pi }}{\zeta }\left(1-K{v}^2\right)-v^2$$

$$>(\epsilon +P+\Pi )\frac{{\tau }_{\Pi }}{\zeta }\left(1-K\right)-1$$

$$(\epsilon +P+\Pi )\frac{{\tau }_{\Pi }}{\zeta }\left(1-K\right)-1\ge 0,$$

Full text

Israel-Stewart relativistic fluid dynamics, causality, general relativity, symmetric hyperbolic.

1. Introduction. The recent detection of a binary neutron star merger using gravitational waves Abbott et al. (2017a) and electromagnetic signals Abbott et al. (2017b) marked the dawn of the multi-messenger astronomy era Abbott et al. (2017c). Such events are expected to provide key information about the properties of matter at extreme densities and temperatures Abbott et al. (2018), as the density of the inner region of the object left over after the merger can be several times larger than the nuclear saturation density ($\sim 0.16$ fm*-3*) while still subject to temperatures of the order of tens of MeV.

Even though the properties of the equation of state of the highly dense matter formed after the merger are still uncertain Most et al. (2019), for many years it was assumed that this system could be reasonably described as an ideal fluid (coupled to Einstein’s equations) since the time scales for viscous transport to set in were previously determined Bildsten and Cutler (1992) to be outside the ten millisecond range, which is the typical time scale associated with damping due to gravitational wave emission. These estimates were recently revisited in Alford et al. (2018) using state-of-the-art merger simulations where it was concluded that, while neutrino-driven thermal transport and shear dissipation remain unlikely to affect the post-merger gravitational wave signal (unless turbulent motion occurs), damping of high-amplitude oscillations due to bulk viscosity is likely to be relevant if direct Urca processes remain suppressed. Therefore, bulk viscosity is expected to play an important role in gravitational wave emission and, as such, it should be taken into account and thoroughly investigated in merger simulations.

However, as stressed in Alford et al. (2018), the effects of bulk viscosity have not yet been included in merger simulations because this requires a formulation of relativistic fluid dynamics including bulk viscosity that is compatible with the underlying causality property of relativity theory in the strong nonlinear regime probed by the mergers (see also Refs. Duez et al. (2004); Shibata and Kiuchi (2017); Radice et al. (2018) for further discussions regarding viscous effects).

In this letter we solve this issue by proving that the equations of motion of Israel-Stewart (IS) type Mueller (1967); Israel (1976); Israel and Stewart (1979), describing a relativistic fluid with bulk viscosity and baryon charge (without shear viscosity or baryon diffusion) dynamically coupled to gravity, are causal in the full nonlinear regime when

[TABLE]

where $\zeta=\zeta(\varepsilon,n)$ is the bulk viscosity, $\tau_{\Pi}=\tau_{\Pi}(\varepsilon,n)$ is the bulk relaxation time, $\Pi$ is the bulk scalar, $\varepsilon$ is the energy density, $P=P(\varepsilon,n)$ is the equilibrium pressure defined by the equation of state, and $n$ is the baryon density. Requirement (1) is a simple nonlinear generalization of the known condition ensuring the linear stability of the IS equations around equilibrium. In fact, we note that near equilibrium, where $|\Pi/(\varepsilon+P)|\ll 1$, and at zero baryon density where the speed of sound squared is $c_{s}^{2}=dP/d\varepsilon$, (1) reduces to $\frac{\zeta}{\tau_{\Pi}(\varepsilon+P)}+\mathcal{O}\left(\frac{\Pi}{\varepsilon+P}\right)\leq 1-c_{s}^{2}$, which is the standard condition for causality and stability in the linearized regime around equilibrium Denicol et al. (2008); Pu et al. (2010); Romatschke (2010).

We also show how to express the equations of motion of this theory that describes a bulk viscous relativistic fluid coupled to gravity as a first-order symmetric hyperbolic (FOSH) system. This immediately implies local existence and uniqueness of its solutions, and sets the equations in a form for which known numerical algorithms can be applied Gustafsson et al. (2013); Reula and Nagy (1997). We stress that our results remain valid if one considers solely the fluid dynamic equations in a fixed (e.g., Minkowski) background. To the best of our knowledge, this is the first time that such statements (causality, local existence, uniqueness) have been proven for IS-like theories in the nonlinear regime.

2. The equations of motion. IS theories of relativistic fluid dynamics Mueller (1967); Israel (1976); Israel and Stewart (1979) were proposed decades ago as a way to solve the long-standing acausality and instability problem of the relativistic Navier-Stokes (NS) equations Hiscock and Lindblom (1985); Pichon (1965). The basic idea is that a dissipative current such as the bulk scalar $\Pi$ does not instantaneously take its NS form111We use units $c=\hbar=k_{B}=1$. The spacetime metric signature is $(-+++)$. Greek indices run from 0 to 3, Latin indices from 1 to 3. $\Pi_{NS}=-\zeta\nabla_{\mu}u^{\mu}$ (where $u_{\mu}$ is the fluid 4-velocity which obeys $u_{\mu}u^{\mu}=-1$) during the fluid evolution; rather, $\Pi$ obeys a relaxation-type equation that describes how it may relax to $\Pi_{NS}$ within the relaxation time scale $\tau_{\Pi}$. While such theories were originally Israel (1976) derived by imposing the second law of thermodynamics for a judicious choice of the out-of-equilibrium entropy density, their modern versions in Refs. Baier et al. (2008) and Denicol et al. (2012) have focused on different aspects of relativistic hydrodynamics. Ref. Baier et al. (2008) emphasized the effective field theory (EFT) character of relativistic hydrodynamics and its applicability in the strong coupling regime (see Finazzo et al. (2015) for a detailed discussion), elucidating the role played by conformal invariance and how that requires the presence of additional terms in the equations of motion that were not usually taken into account in the original IS theory. In Ref. Denicol et al. (2012) a new moment expansion in relativistic kinetic theory Cercignani and Kremer (2002) was employed, together with a power counting scheme involving Knudsen and inverse Reynolds numbers, to derive the equations of motion of hydrodynamics and obtain their corresponding transport coefficients. These new values for the transport coefficients led to an overall improvement with respect to the original IS theory when comparing hydrodynamic calculations to exact solutions of the Boltzmann equation Denicol et al. (2014a, b, c).

The aforementioned versions of the IS theory differ from their original counterpart while retaining its basic physical insights. As described in the previous paragraph, there are different theories representing an improved, albeit different, formulation of the IS equations. These theories will be referred to here as generalized IS theories.

With applications to neutron star mergers in mind, here we only consider the dissipative effects coming from bulk viscosity. The fluid energy-momentum tensor is given by $T^{\mu\nu}=\varepsilon u^{\mu}u^{\nu}+(P+\Pi)\Delta^{\mu\nu}$, where $\Delta_{\mu\nu}=g_{\mu\nu}+u_{\mu}u_{\nu}$ is the projector orthogonal to $u_{\mu}$ and $g_{\mu\nu}$ is the spacetime metric. Using $u^{\mu}u_{\mu}=-1$, energy-momentum conservation, $\nabla_{\mu}T^{\mu\nu}=0$, is equivalent to the equations

[TABLE]

where $\alpha_{1}=\left(\partial P/\partial\varepsilon\right)_{n}$ and $\alpha_{2}=\left(\partial P/\partial n\right)_{\varepsilon}$. Eqs. (2)–(3) are supplemented by the following bulk scalar relaxation equation

[TABLE]

where $\lambda=\lambda(\varepsilon,n)$ is a transport coefficient. Eq. (4) corresponds to Eq. (63) in Denicol et al. (2012) without shear viscosity or baryon diffusion and with $\delta_{\Pi\Pi}=0$. It is common practice to omit the term proportional to $(\nabla_{\alpha}u^{\alpha})^{2}$ in Denicol et al. (2012) and, therefore, we have also done so here Ryu et al. (2015). We choose to work with this specific generalized IS equation because it contains all the relevant physics while making the interpretation of our results clear. In fact, although it is possible to include other terms in the bulk channel Denicol et al. (2012); Romatschke and Romatschke (2017); Philipsen and Schäfer (2014), the above equation contains the essential terms for our discussion: $\tau_{\Pi}$ is associated with causality, $\lambda$ parametrizes the presence of nonlinear terms that do not contain derivatives of the fields, and the term $\zeta\nabla_{\alpha}u^{\alpha}$ ensures that the NS limit can be recovered. In Sec. 2.3 we explain how our methods apply with no change to other generalized IS theories.

We also include a non-zero baryon current $J_{\mu}=n\,u_{\mu}$, whose conservation gives (our results can be easily adapted for $n=0$)

[TABLE]

The fluid is dynamically coupled to gravity via Einstein’s equations

[TABLE]

where $R_{\mu\nu}$ is the Ricci tensor, $R=g_{\mu\nu}R^{\mu\nu}$, $\Lambda$ is the cosmological constant, and $G$ is Newton’s gravitational constant. Eqs. (2)–(6), with $u^{\mu}u_{\mu}=-1$, define the generalized Einstein-Israel-Stewart (EIS) theory considered in this work.

2.1. Causality. Causality is the idea that no information propagates faster than the speed of light and that no closed timelike curves exist in spacetime, i.e., the future cannot influence the past (see the Appendix for the precise definition of causality). This concept lies at the core of relativity theory and, therefore, the dynamics of the fluid sector must be compatible with it. Despite its importance, causality has not been established in the nonlinear regime for the IS theory or its variations (see Bemfica et al. (2018); Disconzi et al. (2017); Rezzolla and Zanotti (2013) for a discussion). While causality is known to hold in ideal relativistic hydrodynamics at the nonlinear level Fourès-Bruhat (1958); Choquet-Bruhat (2009); Disconzi (2015); Lichnerowicz (1967), only statements valid in the linearized regime around equilibrium are known for the IS theory Hiscock and Lindblom (1983, 1985); Olson (1990); Denicol et al. (2008); Pu et al. (2010). Attempts to go beyond the linear regime, although restricted to $1+1$ dimensions or assuming very strong symmetry conditions, and in flat spacetime, appeared, respectively, in Denicol et al. (2008) and Floerchinger and Grossi (2018) (compare (1) with the causality condition found in Denicol et al. (2008)). Thus, a general proof of causality of IS-like systems (coupled to gravity) is so far lacking.

We show that if $\frac{\zeta}{\tau_{\Pi}(\varepsilon+P+\Pi)}+\alpha_{1}+\frac{\alpha_{2}\,n}{\varepsilon+P+\Pi}\geq 0$ and condition (1) is satisfied, then the generalized EIS system is causal. Causality of Eqs. (2)–(5) in a fixed background also holds under the same assumptions.

We refer to Theorem 1 in the Appendix for a formal statement of the causality of the EIS system as well as its proof. In a nutshell, in the Appendix we establish that the values of the fields $\varepsilon,u_{\alpha},\Pi$, and $g_{\alpha\beta}$ at a point $p$ are influenced only by the dynamics of such fields in the causal past of $p$. In Minkowski space, the causal past of $p$ is simply the bottom half of the lighcone with vertex at $p$. This generalizes to curved spacetimes, where the half-bottom of the light cone is replaced by a curved cone-like region with vertex at $p$.

Condition $\frac{\zeta}{\tau_{\Pi}(\varepsilon+P+\Pi)}+\alpha_{1}+\frac{\alpha_{2}\,n}{\varepsilon+P+\Pi}\geq 0$ is physically very natural as $\alpha_{1}+\frac{\alpha_{2}\,n}{\varepsilon+P}$ is the speed of sound squared in equilibrium and $\zeta$ is non-negative.

2.2. The EIS equations as a FOSH system, existence, and uniqueness. A system of first order partial differential equations for an unknown vector $\Phi$ is said to be a FOSH system if it can be written as $\mathcal{A}^{\mu}(\Phi)\partial_{\mu}\Phi+\mathcal{B}(\Phi)=F$, where $\mathcal{A}^{\mu}$ are matrices and $\mathcal{B}$ is a vector, all possibly depending on $\Phi$ but not on its derivatives, $\mathcal{A}^{\mu}$ are symmetric, $\mathcal{A}^{0}$ is positive definite, and $F$ is a given source term.

Several important properties are readily available for FOSH systems Courant and Hilbert (1991); R. P. Geroch (1996); Frittelli and Reula (1996). One such property is that for these systems the initial-value problem admits existence and uniqueness of solutions. Besides assuring a firm theoretical basis for the system, knowing that the equations admit existence and uniqueness of solutions is helpful to ensure the reliability of many numerical schemes; see Guermond et al. (2008) for examples of the potential pitfalls of simulating equations for which no existence and uniqueness results are available, and Bemfica et al. (2018) for a complementary discussion. While for most of standard physical theories existence and uniqueness of solutions had long been settled Choquet-Bruhat (2009), for theories of relativistic fluids with viscosity very few results are available Bemfica et al. (2018); Disconzi (2014, 2019), and none so far for the generalized IS theory in the nonlinear regime.

Eqs (2)–(5) can be written as a FOSH system with $\Phi=(\varepsilon,u^{i},n,\Pi)$, a suitable $\mathcal{B}$, $F=0$, and $\mathcal{A}^{0}$ and $\mathcal{A}^{k}$ given, respectively, by

[TABLE]

and

[TABLE]

where $A_{ij}=g_{ij}-g_{i0}\frac{u_{j}}{u_{0}}-\frac{u_{i}}{u_{0}}g_{0j}+g_{00}\frac{u_{i}u_{j}}{u_{0}^{2}}$.

Using the above formulation of the generalized IS equations as a FOSH, we establish existence and uniqueness of solutions to the initial-value formulation for the generalized EIS system under suitable assumptions. More precisely, given an equation of state $P=P(\varepsilon,n)$, a bulk viscosity $\zeta=\zeta(\varepsilon,n)$, and a relaxation time $\tau_{\Pi}=\tau_{\Pi}(\varepsilon,n)$, consider initial conditions $\mathring{\varepsilon}\equiv\varepsilon(0),\mathring{u}^{i}\equiv u^{i}(0),\mathring{\Pi}\equiv\Pi(0),\mathring{n}\equiv n(0),g_{ij}(0),\partial_{0}g_{ij}(0)$ satisfying (1) along an initial surface $\Sigma\equiv\{x^{0}=0\}$. Assume that $\mathring{\varepsilon}+P(\mathring{\varepsilon},\mathring{n})+\mathring{\Pi}$, $\tau_{\Pi}(\mathring{\varepsilon},\mathring{n})$, $\zeta(\mathring{\varepsilon},\mathring{n})>0$, $\frac{\partial P}{\partial\varepsilon}(\mathring{\varepsilon},\mathring{n})+\frac{\partial P}{\partial n}(\mathring{\varepsilon},\mathring{n})\mathring{n}/(\mathring{\varepsilon}+P(\mathring{\varepsilon},\mathring{n})+\mathring{\Pi})\geq 0$, and that $\mathring{n}$, $\frac{\partial P}{\partial\varepsilon}(\mathring{\varepsilon},\mathring{n})$, and $\frac{\partial P}{\partial n}(\mathring{\varepsilon},\mathring{n})$ are nonzero. Then there exists a spacetime $M$ containing $\Sigma$ such that a unique solution to the generalized EIS equations exists on $M$. A similar statement, i.e., existence and uniqueness of solutions for initial conditions as above, holds for the generalized IS equations in a fixed background.

We refer to Theorem 2 in the Appendix for a formal statement of existence and uniqueness of solutions, its proof, and for the generalized IS as a FOSH. Here, we make some relevant remarks.

Assumptions $\mathring{\varepsilon}+P(\mathring{\varepsilon},\mathring{n})+\mathring{\Pi}$, $\tau_{\Pi}(\mathring{\varepsilon},\mathring{n})$, $\zeta(\mathring{\varepsilon},\mathring{n})>0$, and $\frac{\partial P}{\partial\varepsilon}(\mathring{\varepsilon},\mathring{n})+\frac{\partial P}{\partial n}(\mathring{\varepsilon},\mathring{n})\mathring{n}/(\mathring{\varepsilon}+P(\mathring{\varepsilon},\mathring{n})+\mathring{\Pi})\geq 0$ are physically natural (our results are easily adapted to the case when $n$ is absent), while $\frac{\partial P}{\partial\varepsilon}(\mathring{\varepsilon},\mathring{n})$, $\frac{\partial P}{\partial n}(\mathring{\varepsilon},\mathring{n})\neq 0$ are very mild requirements. Note that these can be viewed as conditions on $\mathring{\varepsilon}$ and $\mathring{n}$, with the form of the equation of state remaining rather general. It is well-known that because of the covariance of Einstein’s equations only $g_{ij}$ and $\partial_{0}g_{ij}$, rather than $g_{\alpha\beta}$ and $\partial_{0}g_{\alpha\beta}$, are given as initial data Hawking and Ellis (1975); Wald (2010). Similarly, only a vector field on $\Sigma$, i.e., $\mathring{u}^{i}$ rather than $\mathring{u}^{\alpha}$, is given along $\Sigma$, with the full four-velocity on $\Sigma$ determined from $u^{\mu}u_{\mu}=-1$.

The proof of Theorem 2 follows from known arguments once the generalized IS equations are formulated as a FOSH system. The main difficulty to achieve the latter is that there is no method to determine whether a given set of equations can in principle be written as a FOSH system. Here, we relied on the following two ingredients to achieve this. First, computing the characteristics of the system (which is used for establishing the causality of the equations), we find them to be a combination of the flow lines of $u^{\alpha}$ and of null-cones with respect to an acoustical metric. The former is associated with transport equations (which are the prototypical example of a FOSH system); the latter is associated with wave equations (which can be rewritten as FOSH). This suggests that a FOSH formulation might be possible. Second, Eqs. (2), (3), and (5) resemble the equations of an ideal fluid, for which a FOSH formulation is known Anile (1990); Anile and Pennisi (1987). This suggests trying to adapt the procedure that works for an ideal fluid. In this regard, it is crucial that the viscous contributions due to bulk viscosity are given by a scalar: if vector (baryon diffusion) or tensor (shear viscosity) viscous contributions are included, the characteristics become very complicated and Eqs. (2), (3), and (5) no longer resemble those of an ideal fluid. It is not clear whether causality (in the nonlinear regime) or a formulation as a FOSH system can be obtained in these cases following the approach pursued here.

We assumed that all transport coefficients depend only on $\varepsilon$ and $n$ (however, see Sec. 2.3). Transport coefficients also generally depend on some characteristic microscopic variables Denicol et al. (2011). In many applications, however, these microscopic variables are either neglected, treated as parameters that can be independently estimated and thus considered given (such as the viscosities in Navier-Stokes theory), or directly eliminated. Otherwise, one would not be dealing with a purely hydrodynamic theory Romatschke and Romatschke (2017). Therefore, our assumptions fall well within the scope of applicability of hydrodynamic models. See also Sec. 2.3 for a more general setting where the dependence on the microscopic dynamics can be decoupled upon simple assumptions.

2.3. Other theories. Here we point out that our methods give, with no change, causality, existence, and uniqueness of solutions to Einstein’s equations coupled to other generalized IS theories (without shear viscosity or baryon diffusion). A quick inspection in our proofs provided in the Appendix reveals that they remain true, with no change, if $\tau_{\Pi}$ and $\zeta$ are allowed to depend on $\Pi$. In particular, consider Eq. (63) in Denicol et al. (2012) with $\delta_{\Pi\Pi}$ not necessarily zero. Defining $\bar{\zeta}=\zeta+\delta_{\Pi\Pi}\Pi$, we obtain Eq. (4) with $\zeta$ replaced by $\bar{\zeta}$. Hence, all our results go through with $\bar{\zeta}$ in place of $\zeta$. Our assumptions in this case, however, seem less natural since, for example, there is no reason to expect $\bar{\zeta}$ to be positive.

Finally, the methods we employ depend only on the principal part of the system222I.e., the subsystem comprised only of the terms involving derivatives in the equations.. Thus all our conclusions remain valid if Eq. (4) is replaced by

[TABLE]

where $f$ is an arbitrary function of $\varepsilon$, $n$, $\Pi$, $\tau_{\Pi}$, and $\zeta$ (but not depending on their derivatives) and $\tau_{\Pi}$ and $\zeta$ are allowed to depend on $\Pi$ (so that, in particular, $\zeta$ can be replaced by $\bar{\zeta}$ here as well).

Consider now the case when $\zeta$ and $\tau_{\Pi}$ depend on a microscopic variable $\ell$, and suppose that $\ell$ cannot be neglected or treated as a parameter but needs to be found by solving its own evolution equation which couples to the hydrodynamic variables. If one assumes, as it is often done Ryu et al. (2015), that the ratio $\zeta/\tau_{\Pi}$ depends only $\varepsilon$, $n$, and $\Pi$, then all our conclusions remain valid as long as $f$ is such that $\frac{1}{\tau_{\Pi}}f(\varepsilon,n,\Pi,\tau_{\Pi},\zeta)$ does not involve dynamical microscopic variables. Thus, our results are applicable to many cases of interest in heavy ion collisions (e.g., Denicol et al. (2018)) and neutron star mergers. In this regard, note that our condition (1) depends only on the ratio $\zeta/\tau_{\Pi}$.

3. Conclusions. In this paper we proved that the generalized EIS equations for a relativistic fluid with bulk viscosity and nonzero baryon charge (without shear viscosity or baryon diffusion) are causal in the full nonlinear regime. We also proved that these equations admit a FOSH formulation, existence, and uniqueness of solutions. Our results hold under hypothesis typically satisfied by viscous relativistic fluids. In particular, we do not require any symmetry or near-equilibrium assumption. The most immediate application of these results is in the study of neutron star mergers. In this regard, it is crucial that our results hold under general assumptions on the equation of state since the latter is not yet known from first principles.

Given that theories of IS type are the most widely used formulation of relativistic viscous fluid dynamics Jeon and Heinz (2015), it is crucial to put them on a firm theoretical (and mathematical) foundation. Combining our results with known properties of the generalized IS equations Israel and Stewart (1979); Pu et al. (2010); Romatschke (2010); Denicol et al. (2012), the case considered here seems to be the first example in the literature of a theory of relativistic viscous fluids such that (i) causality, local existence, and uniqueness of solutions hold in the nonlinear regime with or without coupling to Einstein’s equations; (ii) linear stability around equilibrium holds; (iii) the equations of motion are derivable from microscopic approaches (such as kinetic theory); and (iv) solutions are guaranteed to exist in function spaces (the Sobolev spaces, see the Appendix) well-suited for the implementation of numerical codes. (The theory introduced in Bemfica et al. (2018) also satisfies (i)-(iii), but it is applicable only to conformal fluids and its solutions exist on function spaces more restrictive than Sobolev spaces.) This is the first time that all such properties are shown to hold since Eckart’s seminal work in 1940 Eckart (1940).

As Theorems 1 and 2 settle the basic question of causality and existence of solutions to the EIS system, it is now possible to generalize to relativistic viscous fluids key results known to hold for ideal fluids, such as the formation of shocks Christodoulou (2007), global-in-time results Speck (2012), or the problem of accurately describing the fluid-vacuum interface on stars Rendall (1992).

These results have two further applications: the hydrodynamic evolution of the quark-gluon plasma formed in heavy ion collisions Heinz and Snellings (2013) and cosmology. Current state-of-the-art modeling of quark-gluon plasma dynamics involves solving IS-like equations taking into account shear and bulk viscosities and, more recently, baryon diffusion Denicol et al. (2018). Differently than the case of neutron star mergers Alford et al. (2018), in heavy ion collisions shear and baryon diffusion effects are not negligible though bulk viscosity also plays an important role Noronha-Hostler et al. (2013, 2014); Ryu et al. (2015). In this regard, we note Theorem 2 guarantees existence of solutions even in the presence of cavitation where $P+\Pi\sim 0$ Torrieri et al. (2008); Rajagopal and Tripuraneni (2010) and $\varepsilon>0$, a phenomenon whose numerical description can be quite challenging Denicol et al. (2015). Another application of the EIS system studied here is in cosmology, where the inclusion of bulk viscosity has been widely studied Maartens (1995); Li and Barrow (2009); Disconzi et al. (2015); Zimdahl et al. (2001); Piattella et al. (2011). Indeed, our results do not make any symmetry assumptions, thus allowing the study of cosmological models with viscosity outside the symmetry class of homogeneous and isotropic solutions (e.g. Montani and Venanzi (2017); East et al. (2016)).

The equations here investigated suffice for many important applications (including neutron star mergers). Recent developments in EFTs, however, indicate that a more complete description of fluid phenomena should include further terms in the equations. While such contributions are negligible in many cases of interest, they consist of a genuine prediction of EFT, and their omission can lead to inconsistencies in 2${}^{\text{nd}}$ order Kubo formulas Romatschke and Romatschke (2017); Philipsen and Schäfer (2014). Generalizing our results to theories that include these effects is, therefore, important for a more complete description of relativistic viscous fluids. While our techniques do not apply directly to such cases, we expect them to provide a starting point for studying these scenarios.

Our results open the door for new in-depth studies of fluid dynamics under extreme conditions such as in relativistic turbulent phenomena Eyink and Drivas (2018) and neutron star mergers. In particular, the latter is expected to take center stage in the coming years and our causality and existence results, which remain valid in the full nonlinear regime, provide the necessary cornerstone for quantitative studies of non-equilibrium viscous fluid phenomena in strong gravitational fields.

Acknowledgements. We would like to thank the anonymous referees for a careful reading of the manuscript and several comments that improved the paper. FSB is partially supported by a Discovery Grant administered by Vanderbilt University. MMD is partially supported by a Sloan Research Fellowship provided by the Alfred P. Sloan foundation, NSF grant DMS-1812826, and a Discovery Grant administered by Vanderbilt University. JN is partially supported by CNPq grant 306795/2017-5 and FAPESP grant 2017/05685-2.

Appendix A Appendix

Appendix B Causality

Here we give a precise formulation of the causality of the generalized EIS system and its proof.

Theorem 1**.**

Let $(\varepsilon,u^{\alpha},\Pi,n,g_{\alpha\beta})$ be a solution to the generalized EIS equations defined on a globally hyperbolic spacetime $M$, and let $\Sigma$ be a Cauchy surface. Suppose $\frac{\zeta}{\tau_{\Pi}(\varepsilon+P+\Pi)}+\alpha_{1}+\frac{\alpha_{2}\,n}{\varepsilon+P+\Pi}\geq 0$ and that (1) holds. Then, for any $p\in M$ in the future of $\Sigma$, $(\varepsilon(p),u^{\alpha}(p),\Pi(p),n(p),g_{\alpha\beta}(p))$ depends only on $\left.(\varepsilon,u^{\alpha},\Pi,n,g_{\alpha\beta},\kappa_{\alpha\beta})\right|_{\Sigma\cap J^{-}(p)}$, where $J^{-}(p)$ is the causal past of $p$ and $\kappa$ is the extrinsic curvature of $\Sigma$ in $M$.

Theorem 1 gives a precise mathematical formulation of the statement that in the generalized EIS system faster-than-light signals are absent and the future cannot influence the past. The concept of causality is properly formulated in globally hyperbolic spacetimes since even solutions to vacuum Einstein’s equations can otherwise fail to be causal Hawking and Ellis (1975); Wald (2010). A Cauchy surface plays the role of the set $\{x^{0}=\text{constant}\}$ when coordinates are introduced, with $x^{0}$ playing the role of time. However, Theorem 1 is formulated with no mention to coordinates because a precise statement of causality in general relativity has to be coordinate free. The causal past of $p$, $J^{-}(p)$, is the general relativistic analogue of the past light-cone at a point Hawking and Ellis (1975); Wald (2010). The conclusion of the theorem says that the values of the fields $\varepsilon,u_{\alpha},\Pi$, and $g_{\alpha\beta}$ at $p$ are determined only by the the dynamics of such fields in the (causal) past of $p$. The curvature $\kappa_{\alpha\beta}$ can be thought of as the normal derivative of $g_{\alpha\beta}$ along $\Sigma$. To understand the role of $\kappa_{\alpha\beta}$ note that, in the most common situation, $\Sigma$ is a surface where initial data for the generalized EIS system is provided. In this case the metric and its first derivatives have to be prescribed on $\Sigma$ since Einstein’s equations are of second order.

Proof of Theorem 1: To prove Theorem 1, we need to compute the characteristics of the generalized EIS system Courant and Hilbert (1991). First, we write equations (2), (3), (4), and (5) in matrix form as $\mathcal{M}^{\mu}\nabla_{\mu}\Psi+\mathcal{N}\Psi=0$, where $\Psi=(\varepsilon,u^{\nu},n,\Pi)$ and the matrices $\mathcal{M}^{\mu}$ are given by

[TABLE]

and $\mathcal{N}$ is given by

[TABLE]

Fix a point $p$ in spacetime. Choose coordinates $\{x^{\mu}\}$ near $p$. Expanding the covariant derivative, the system reads $\mathcal{M}^{\mu}\partial_{\mu}\Psi+\mathcal{M}_{\Gamma}\Psi+\mathcal{N}\Psi=0$, where $\mathcal{M}_{\Gamma}$ denotes the matrix with terms involving the Christoffel symbols coming from expanding $\nabla_{\mu}$.

Letting $\xi_{\mu}$ be an arbitrary covector in the cotangent space at $p$, the characteristic determinant of the system is $\mathcal{P}(\xi)=\det(\mathcal{M}^{\mu}\xi_{\mu})$. Computing, we find

[TABLE]

The system’s characteristics are determined by the $\xi_{\mu}$’s such that $P(\xi)=0$. Assuming condition (1), we have $\varepsilon+P+\Pi\neq 0$, $\tau_{\Pi}\neq 0$, and the characteristics are thus determined by $u^{\alpha}\xi_{\alpha}=0$ or by the solutions coming from setting the term in square brackets in (7) to zero. The characteristics corresponding to $u^{\alpha}\xi_{\alpha}=0$ are simply the flow lines of $u^{\alpha}$, which are causal as long as $u^{\alpha}$ remains timelike. The characteristics associated with the vanishing of the bracket in (7) are precisely the characteristics of an acoustical metric Christodoulou (2007); Disconzi and Speck (2018) with an effective sound speed squared given by $c_{\text{acoustical}}^{2}=\frac{\zeta}{\tau_{\Pi}(\varepsilon+P+\Pi)}+\alpha_{1}+\frac{\alpha_{2}\,n}{\varepsilon+P+\Pi}\geq 0$. Thus, the system will be causal as long as $u^{\alpha}$ remains timelike (and normalized) and $c_{\text{acoustical}}^{2}\leq 1$, which is precisely condition (1).

The discussion above has not taken into account Einstein’s equations (6). We consider now the full system (2), (3), (4), (5), and (6), and assume that our coordinates about $p$ are wave coordinates Choquet-Bruhat (2009), in which case Einstein’s equations can be written as

[TABLE]

where $\mathcal{Z}_{\alpha\beta}$ involves only $g_{\alpha\beta}$ and its first derivatives. Since $T_{\alpha\beta}$ involves no derivatives of $\Psi$ and equations (2), (3), (4), (5), involve at most one derivative of $g_{\alpha\beta}$ (via the Christoffel symbols that have been grouped in the matrix $M_{\Gamma}$), the characteristic determinant of the system is now $\mathcal{P}(\xi)\mathcal{Q}(\xi)$, where $\mathcal{P}(\xi)$ is as above and $\mathcal{Q}(\xi)=(g^{\mu\nu}\xi_{\mu}\xi_{\nu})^{10}$. $\mathcal{Q}(\xi)=0$ corresponds to the characteristics associated with the light-cones, hence causal. Standard arguments from general relativity show that the characteristics computed in this fashion are geometric or coordinate independent Wald (2010). ∎

Appendix C Proof of Theorem 2

Denote by $H_{\text{ul}}^{s}$ the uniformly local Sobolev space (functions whose distributional derivatives up to order $s$ exist and are square-integrable in a partition), and by $C^{s}$ the space of continuously diffentiable functions up to order $s$ Choquet-Bruhat (2009); Kato (1975). We have:

Theorem 2**.**

Let $\mathcal{I}=(\Sigma,\mathring{\varepsilon},\mathring{u}^{i},\mathring{\Pi},\mathring{n},\mathring{g}_{ij},\mathring{\kappa}_{ij})$ be an initial data set for the generalized EIS system, with an equation of state $P=P(\varepsilon,n)$, a bulk viscosity $\zeta=\zeta(\varepsilon,n)$, and a relaxation time $\tau_{\Pi}=\tau_{\Pi}(\varepsilon,n)$. Assume that $\mathring{\varepsilon}+P(\mathring{\varepsilon},\mathring{n})+\mathring{\Pi}$, $\tau_{\Pi}(\mathring{\varepsilon},\mathring{n})$, $\zeta(\mathring{\varepsilon},\mathring{n})>0$, $\frac{\partial P}{\partial\varepsilon}(\mathring{\varepsilon},\mathring{n})+\frac{\partial P}{\partial n}(\mathring{\varepsilon},\mathring{n})\mathring{n}/(\mathring{\varepsilon}+P(\mathring{\varepsilon},\mathring{n})+\mathring{\Pi})\geq 0$, and that $\mathring{n}$, $\frac{\partial P}{\partial\varepsilon}(\mathring{\varepsilon},\mathring{n})$, and $\frac{\partial P}{\partial n}(\mathring{\varepsilon},\mathring{n})$ are nonzero. Suppose that $\mathring{g}_{ij}\in H^{s+1}_{\text{ul}}(\Sigma)$, $\mathring{\varepsilon},\mathring{u}^{i},\mathring{\Pi},\mathring{n},\mathring{\kappa}_{ij}\in H^{s}_{\text{ul}}(\Sigma)$, and that $P,\zeta,\lambda,\tau_{\Pi}\in C^{s}(\mathbb{R}^{2})$, where $s\geq 3$. Suppose that (1) holds for $\mathcal{I}$. Then, there exists a globally hyperbolic development of $\mathcal{I}$.

An initial data set for the generalized EIS system involves a three-dimensional manifold $\Sigma$ that plays the analogue of a slice of “constant time” (this is just an analogy since constant time is a coordinate dependent notion). The initial data along $\Sigma$ cannot be arbitrary but need to satisfy the constraint equations Hawking and Ellis (1975); Wald (2010).

Since a solution to Einstein’s equations is not a set of fields on spacetime but spacetime itself, existence and uniqueness of solutions is properly stated in the context of general relativity as the existence of a globally hyperbolic development. Also, for Einstein’s equations, uniqueness is meant in a geometric sense Hawking and Ellis (1975); Wald (2010).

We begin showing that equations (2), (3), (4), and (5) can be written as a FOSH system. We use that $\nabla_{\alpha}u^{0}=-\frac{u_{i}}{u_{0}}\nabla_{\alpha}u^{i}$, multiply (3) by $C_{i}^{\alpha}=g_{i}^{\alpha}-\frac{u_{i}}{u_{0}}\delta^{\alpha}_{0}$, and use that $C_{j}^{\alpha}u_{\alpha}=0$. Assuming that $\zeta,\tau_{\Pi},(\varepsilon+P+\Pi)>0$ and $n,\left(\frac{\partial P}{\partial\varepsilon}\right)_{n},\left(\frac{\partial P}{\partial n}\right)_{\varepsilon}\neq 0$, we readily see that equations (2), (3), (4), and (5) assume the form $\mathcal{A}^{\mu}\nabla_{\mu}\Phi+\mathcal{C}\Phi=0$, where the matrices $\mathcal{A}^{\mu}$ and $\Phi$ are as in Section 2.2, and $\mathcal{C}$ is the matrix that has all entries equal to zero except for $\mathcal{C}_{6,6}=\frac{1+\lambda\Pi}{\zeta\left(\frac{\partial P}{\partial\varepsilon}\right)_{n}}$. Expanding the covariant derivative we have $\mathcal{A}^{\mu}\partial_{\mu}\Phi+\mathcal{A}_{\Gamma}\Phi+\mathcal{C}\Phi=0$, where $\mathcal{A}_{\Gamma}$ denotes the matrix with terms involving Christoffel symbols. Thus, we have $\mathcal{A}^{\mu}\partial_{\mu}\Phi+\mathcal{B}=0$, with $\mathcal{B}=\mathcal{A}_{\Gamma}\Phi+\mathcal{C}\Phi$. Note that $\mathcal{B}$ involves no derivatives of $\Phi$. Below, we show that $\mathcal{A}^{0}$ is positive definite.

Embed $\Sigma$ into $\mathbb{R}\times\Sigma$ and fix a point $p\in\Sigma$. We consider Einstein’s equations written in wave coordinates in a neighborhood of $p$, i.e., equation (8). Without loss of generality we can assume that the coordinates are such that $\{x^{i}\}$ are coordinates on $\Sigma$ and that $g_{\alpha\beta}$ is the Minkowski metric at $p$ Wald (2010). It has long been known that Einstein’s equations can be written as a FOSH system Fischer and Marsden (1972). Indeed, setting $\Theta=(g_{\alpha\beta},\partial_{i}g_{\alpha\beta},\partial_{0}g_{\alpha\beta})$, (8) reads $\widetilde{A}^{\mu}\partial_{\mu}\Theta+\widetilde{\mathcal{B}}=0$, where $\widetilde{\mathcal{A}}^{0}$ and $\widetilde{\mathcal{A}}^{k}$ are given by, respectively,

[TABLE]

and

[TABLE]

where $1_{10}$ and $0_{10}$ are the $10\times 10$ identity and zero matrices, respectively, and $\widetilde{\mathcal{B}}$ contains the terms on the right-hand side of (8). The full system thus reads $\mathfrak{A}^{\mu}\partial_{\mu}\Xi+\mathfrak{B}=0$, where $\Xi=(\Phi,\Theta)$, $\mathfrak{B}=(\mathcal{B},\widetilde{\mathcal{B}})$, and

[TABLE]

The matrices $\mathfrak{A}^{\mu}$ are symmetric and $\mathfrak{B}$ does not involve any derivative of $\Xi$. In $\mathfrak{B}$ we have already eliminated $u^{0}$ appearing on the right-hand side of (8) in terms of $u^{i}$, so that $\mathfrak{B}$ depends only on $\Xi$. To see that such elimination is indeed possible, note that at $p$ the normalization $u^{\mu}u_{\mu}=-1$ reads $(u^{0})^{2}=1+\delta_{ij}u^{i}u^{j}$, and we define $u^{0}$ as the positive solution (since $u^{\mu}$ is future directed). Below we argue that the fields vary continuously with the spacetime coordinates. Thus, by continuity, the equation $u^{\mu}u_{\mu}=-1$ has precisely two solutions $u^{0}$ in a neighborhood of $p$, one, and only one, of which determines a future pointing vector.

Since $\widetilde{\mathcal{A}}^{0}$ is known to be positive definite Fischer and Marsden (1972), to conclude that $\mathfrak{A}^{0}$ is positive definite it suffices to analyze $\mathcal{A}^{0}$. We do this by verifying that all the leading principal minors of $\mathcal{A}^{0}$ are positive. The $5\times 5$ upper left part of $\mathcal{A}^{0}$ has all minors positive, since it corresponds to the ideal fluid which is already known to form a symmetric hyperbolic system for $(\varepsilon,u^{i},n)$.

It remains to compute the determinant of $\mathcal{A}^{0}$. For this, we first compute $\det\mathcal{A}^{0}$ at $p$, in which case $u_{\mu}=1/\sqrt{1-v^{2}}(-1,\mathbf{v})$ (recall that $g_{\alpha\beta}$ is the Minkowski metric at $p$). Then

[TABLE]

where $v^{2}=\mathbf{v}^{2}$,with $v\in[0,1)$, and $K=\alpha_{1}+\frac{\alpha_{2}\,n}{\varepsilon+P+\Pi}$. Since

[TABLE]

if $K\geq 0$, we see that $\det\mathcal{A}^{0}>0$ if

[TABLE]

but this later condition is the same as (1). Because the fields are continuous functions of the spacetime coordinates, we conclude that $\det\mathcal{A}^{0}>0$ will also hold in a neighborhood of $p$.

In order to invoke theorems of hyperbolic differential equations to establish Theorem 2, we need to check that the system is FOSH when the fields take the initial data (the solution $\Xi$ has yet to be shown to exist; the only fields at our disposal at this point are coming from the initial data). But this is immediate since $\mathfrak{A}$ and $\mathfrak{B}$ do not depend on derivatives of $\Xi$ and the foregoing conditions hold for the initial data by hypothesis. Furthermore, the continuity required for the foregoing arguments is true because functions in $H^{s}$ are in particular continuous when $s\geq 3$ in view of the standard Sobolev embedding theorem Evans (2010). By standard existence and uniqueness theorems for FOSH systems Choquet-Bruhat (2009); Majda (1984) we obtain a solution in a neighborhood of $p$. It is known that a solution to Einstein’s equations in wave coordinates yields a solution to the full Einstein equations if and only if the constraints are satisfied, which is the case by hypothesis. Moreover, defining $u^{0}$ via the relation $u^{\mu}u_{\mu}=-1$ (which uniquely defines $u^{0}$ in a neighborhood of $p$, as above) and working backwards, we see that the original equations (2), (3), (4), (5), and (6) are satisfied in a neighborhood of $p$. Finally, a solution valid in a neighborhood of $\Sigma$ is obtained by a standard argument of constructing solutions as above in the neighborhood of different points $p$ and gluing them via diffeomorphisms Wald (2010). Standard properties of Einstein’s equations give the existence of the desired globally hyperbolic development Choquet-Bruhat (2009).