# First-Order Viscous Relativistic Hydrodynamics on the Two-Sphere

## TL;DR

This paper presents numerical solutions to the BDNK first-order viscous relativistic hydrodynamics equations on a spherical surface, demonstrating their behavior and potential singularity formation in finite time for conformal fluids.

## Contribution

It develops a finite difference numerical method on cubed-sphere coordinates for BDNK equations and explores their solutions, including singularity formation, on a spherical geometry.

## Key findings

- Solutions can develop finite-time singularities from smooth initial data.
- Numerical methods are extendable to full 3D simulations including radial variations.
- Flow diverges from equilibrium under certain viscous conditions.

## 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.