First-Order Viscous Relativistic Hydrodynamics on the Two-Sphere

Abstract

A few years ago, Bemfica, Disconzi, Noronha, and Kovtun (BDNK) formulated the first causal, stable, strongly hyperbolic, and locally well-posed theory of first-order viscous relativistic hydrodynamics. Since their inception, there have been several numerical and analytic studies of the BDNK equations which have revealed their promise in modeling relativistic flows when viscous, first-order corrections to ideal hydrodynamics are important. In this paper, we present numerical solutions to the BDNK equations for a $4$D conformal fluid in Minkowski spacetime constrained to the surface of a geometric sphere. We numerically solve the underlying equations of motion by use of finite difference methods applied in cubed-sphere coordinates -- a multi-block grid structure which regularly and continuously covers the surface of a sphere. We present three test cases of our code: linearized fluid perturbations of equilibrium states, a smooth, stationary initial Gaussian pulse of energy density, and Kelvin-Helmholtz-unstable initial data. In the Gaussian test case with sufficiently large entropy-normalized shear viscosity, the flow, though initialized in equilibrium, dynamically diverges away from equilibrium and the regime of validity of first-order hydrodynamics as very steep gradients form in the solution, causing convergence to be lost in the numerical simulation. This behavior persists at all grid resolutions we have considered, and also occurs at much higher resolutions in simulations of planar-symmetric ($1+1$)D conformal flows. These solutions provide numerical evidence that singularities in solutions to the BDNK equations can form in finite time from smooth initial data. The numerical methods we employ on the two-sphere can be readily extended to include variations in the radial direction, allowing for full ($3+1$)D simulations of the BDNK equations in astrophysical applications.

Figures (20)

• Figure 1: Illustration of the cubed-sphere grid on which we discretize and solve the Euler and BDNK equations. Left: a geometric illustration of covering the two-sphere with six regular, non-overlapping patches. Each patch is constructed by projecting onto the sphere the solid angle of each face of a cube located at the center of the sphere. Right: the resulting numerical grid structure on which the Euler and BDNK equations are discretized and solved. Points in each of the six regions on the sphere can be projected onto the corresponding face of the cube, allowing one to transform to local coordinates $X$ and $Y$ [defined in Eqs. (\ref{['eq:CS1']}-\ref{['eq:CSMetric']})] on the cubed-sphere grid and there discretize and solve the underlying equations of motion. We take the grid to be uniform and discretize the equations of motion in space using centered finite-difference stencils.
• Figure 2: Energy density at the north pole from numerical simulations of the even-perturbed initial data (\ref{['eq:Mode:Eps']}, \ref{['eq:Mode:EvenVel']}) to the Euler (blue) and BDNK equations with $\eta/s=1/(4\pi)$, where the red and green curves correspond to the under and overdamped solutions (\ref{['eq:BDNKEvenPertSol1']}, \ref{['eq:BDNKEvenPertSol2']}), respectively. As illustrated in Tab. \ref{['tab:ModeAnalysis']}, the oscillation and damping response in these numerical simulations agree well with the predictions from the linearized equations of motion.
• Figure 3: Energy density and polar velocity profiles of an inviscid [$\eta/s=0/(4\pi)$] flow with Gaussian initial data \ref{['eq:2DGaussian:InitialData']}. Upper left: energy density initial data. Upper right: snapshot of the energy density just before the peak has arrived at the south pole. Lower left (right): energy density (polar velocity) just before once again peaking at the north pole. As the energy density peak traverses from pole to pole, increasingly steep gradients form in the solution, culminating in the formation of a shock at $t\approx 8$ as the energy density is focused into an extremely narrow peak at the north pole (see left panel of Fig. \ref{['fig:2DGaussian:EulerShock']}).
• Figure 4: Energy density at $t=8$ evolved from the Gaussian initial data \ref{['eq:2DGaussian:InitialData']} with $\eta/s=0/(4\pi)$ (left) and $\eta/s=1/(4\pi)$ (right). At this time, convergence is lost in the numerical simulations of the inviscid flow due to the formation of a shock as the energy density is focused into a very narrow peak at the north pole. As shown in the right panel of this figure, viscosity smooths the shock that forms in the inviscid solution. Convergence is maintained throughout the $\eta/s=1/(4\pi)$ simulation, in which a steady state is reached at late times.
• Figure 5: Mode decomposition (up to $l=4$) of the energy density for flows with Gaussian initial data \ref{['eq:2DGaussian:InitialData']} and $\eta/s\in\{0,1,3,10\}$ (upper left, upper right, lower left, lower right). In the $\eta/s=1/(4\pi)$ solution, viscosity serves to damp the subdominant, higher-frequency modes present in the inviscid flow, preventing the formation of a shock and allowing a steady state to be reached at late times in the numerical simulation. As $\eta/s$ is further increased, the relative power of the higher-frequency modes grows as steeper gradients form in the solution, causing convergence to be lost in the $\eta/s\in\{10,20\}/(4\pi)$ solutions at all grid resolutions we have considered. We note that the $\eta/s\in\{0,10\}/(4\pi)$ simulations lose convergence at $t\approx8$ and $t\approx6$, respectively, so the corresponding $C_{l}$ coefficients should be treated with caution at later times. Further, we have not performed a convergence test on the $C_{l}$, so we restrict to only a qualitative discussion of the mode structure.