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.

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.

Tables

Table 1. Table 1: Simulation parameters for the four sets of initial data considered in Sec.  IV and App.  A . For each set of initial data, we list the total number of spatial grid points N tot N_{\mathrm{tot}} in the highest resolution run, the amplitude coefficient α KO \alpha_{\mathrm{KO}} of the added Kreiss-Oliger dissipation, and the CFL factor C ≡ Δ ​ t / h C\equiv\Delta{t}/h used for time stepping, where Δ ​ t \Delta{t} is the time step and h h is the grid spacing. References to the total number of grid points on the two-sphere corresponds to a total of 6 ​ N 2 6N^{2} points across all six patches.

Initial Data	$N_{\mathrm{tot}}$	$C$	${\alpha }_{\mathrm{KO}}$
$1$D Gaussian	$2^{14}$	$0.1$	$0.5$
$2$D Fluid Perturbations	$2^{12}$	$0.2$	$0.5$
$2$D Gaussian	$2^{14}$	$0.2$	$0.5$
$2$D Kelvin-Helmholtz Instability	$2^{16}$	$0.1$	$0.5$

Table 2. Table 2: Percentage residuals between the predicted and fitted oscillation frequencies and damping times of the even perturbations ( 27 , 31a ) to the Euler and BDNK equations. For the highest resolution fits, the standard error of the fitted Euler frequency is ∼ 10 − 4 % \sim 10^{-4}\% and, for the BDNK simulations I and II, are { 10 − 4 % , 10 − 3 % } \{10^{-4}\%,10^{-3}\%\} and { 10 − 2 % , 10 − 2 % } \{10^{-2}\%,10^{-2}\%\} for the oscillation frequency and damping time, respectively. From the linearized equations of motion, we expect an oscillation frequency ω osc = 2 / 3 \omega_{\rm{osc}}=\sqrt{2/3} in the Euler fluid, and, for BDNK simulations I and II, frequencies ω osc = { 0.82 , 0.73 } \omega_{\textrm{osc}}=\{0.82,0.73\} and damping times τ d = { 20.37 , 0.66 } \tau_{\mathrm{d}}=\{20.37,0.66\} , respectively, where we have quoted numerical values to two decimal places. The fitted parameters are obtained by fitting the time series of the energy density at the north pole to purely oscillatory or exponentially damped sinusoids. We find better agreement between the fitted and predicted parameters as the grid resolution is increased, although, as seen in the rightmost two columns, the overdamped oscillations in simulation II make the extraction of the oscillation frequency more difficult relative to the other two simulations.

	$\frac{\eta }{s}=\frac{0}{4\pi }$	$\frac{\eta }{s}=\frac{1}{4\pi }$, Sim I		$\frac{\eta }{s}=\frac{1}{4\pi }$, Sim II
$N$	$\Delta {\omega }_{\mathrm{osc}}$	$\Delta {\omega }_{\mathrm{osc}}$	$\Delta {\tau }_d$	$\Delta {\omega }_{\mathrm{osc}}$	$\Delta {\tau }_d$
$2^3+1$	$0.004\%$	$0.6\%$	$20\%$	$80\%$	$37\%$
$2^4+1$	$0.0005\%$	$0.03\%$	$0.9\%$	$1.1\%$	$1.9\%$
$2^5+1$	$0.0003\%$	$0.001\%$	$0.03\%$	$1.1\%$	$0.9\%$

Table 3. Table 3: Fitted values of the time t d t_{d} obtained by fitting the time series of max ⁡ { | ϵ ′ | } \max\{|\epsilon^{\prime}|\} for η / s = 20 / ( 4 ​ π ) \eta/s=20/(4\pi) (shown in Fig.  16 ) to the functional form ( 46 ) in order to estimate the time at which spatial gradients in the continuum solution become infinite.

$N$	Fitted $t_d$
$2^{12}+1$	$45.9$
$2^{13}+1$	$43.4$
$2^{14}+1$	$42.5$

Equations

$${\nabla }_{\mu }{T}^{\mu \nu }=0,{\nabla }_{\mu }{J}^{\mu }=0.$$

$$T^{\mu \nu }$$

$$J^{\mu }$$

$$\mathcal{E}\equiv u_{\mu }{u}_{\nu }{T}^{\mu \nu },\mathcal{P}\equiv \frac{1}{d}{\Delta }_{\mu \nu }{T}^{\mu \nu },\mathcal{N}\equiv -u_{\mu }{J}^{\mu }$$

$${\mathcal{Q}}^{\mu }\equiv -{\Delta }^{\mu }{}_{\alpha }{u}_{\beta }{T}^{\alpha \beta },{\mathcal{J}}^{\mu }\equiv {\Delta }^{\mu }{}_{\alpha }{J}^{\alpha },$$

$${\mathcal{T}}^{\mu \nu }\equiv \frac{1}{2}\left({\Delta }^{\mu }{}_{\alpha }{\Delta }^{\nu }{}_{\beta }+{\Delta }^{\nu }{}_{\alpha }{\Delta }^{\mu }{}_{\beta }-\frac{2}{d}{\Delta }^{\mu \nu }{\Delta }_{\alpha \beta }\right)T^{\alpha \beta },$$

$$T^{\mu \nu }=T_{(0)}^{\mu \nu }+T_{(1)}^{\mu \nu }+T_{(2)}^{\mu \nu }+\dots ,$$

$$J^{\mu }=J_{(0)}^{\mu }+J_{(1)}^{\mu }+J_{(2)}^{\mu }+\dots ,$$

$$T_{(0)}^{\mu \nu }$$

$$J_{(0)}^{\mu }$$

$$T^{\mu \nu }$$

$$+{\mathcal{Q}}^{\mu }{u}^{\nu }+u^{\mu }{\mathcal{Q}}^{\nu }-2\eta {\sigma }^{\mu \nu },$$

$$J^{\mu }$$

$$\mathcal{A}$$

$$\Pi$$

$${\mathcal{Q}}^{\mu }$$

$${\sigma }^{\mu \nu }$$

$${\beta }_{ϵ}$$

$${\beta }_n$$

$$P(ϵ,n)=\frac{ϵ}{3},{\tau }_P=\frac{{\tau }_{ϵ}}{3},\zeta =0,$$

$${\beta }_{ϵ}=\frac{{\tau }_Q}{3},{\beta }_n=0,$$

$${\tau }_{ϵ}=\frac{3\chi }{4ϵ},{\tau }_Q=\frac{3\lambda }{4ϵ},$$

$$\eta >0,\chi =a_1\eta ,\lambda \ge \frac{3\eta a_1}{a_1-1}{a}_1\ge 4.$$

$$\eta ={\eta }_0{ϵ}^{3/4},\chi ={\chi }_0{ϵ}^{3/4},\lambda ={\lambda }_0{ϵ}^{3/4},$$

$$\left({\lambda }_0,{\chi }_0\right)$$

$${\eta }_0=\frac{4{\epsilon }_0^{1/4}}{3}\cdot \frac{\eta }{s},$$

$${\left(\frac{\eta }{s}\right)}_{\min }=\frac{\hslash }{4\pi k_B},$$

$$\mathrm{d}{s}^2=-\mathrm{d}{t}^2+\mathrm{d}{r}^2+R^2\mathrm{d}{\theta }^2+R^2{\sin }^2\theta \mathrm{d}{\varphi }^2,$$

$$u^{\mu }={\left[u^t,0,u^{\theta },u^{\varphi }\right]}^T,$$

$$u^t=\sqrt{1+R^2(u^{\theta }{)}^2+R^2(u^{\varphi }{)}^2{\sin }^2\theta }.$$

$$f(\theta ,\varphi )=\sum _{l=0}^{\infty }\sum _{m=-l}^l{A}_{lm}{Y}_l^m(\theta ,\varphi ),$$

$${\int }_{S^2}{Y}_l^m{Y}_{l^{\prime }}^{m^{\prime }}\mathrm{d}\Omega =4\pi {\delta }_{l{l}^{\prime }}{\delta }_{m{m}^{\prime }}.$$

$$\frac{1}{4\pi }{\int }_{S{}^2}|f(\theta ,\varphi ){|}^2\mathrm{d}\Omega$$

$$C_l$$

$${\omega }_{\mu \nu }$$

$${\int }_{S{}^2}{W}^0\mathrm{d}\Omega$$

Full text

I Introduction

Relativistic hydrodynamics can be understood as a long-wavelength effective theory describing the evolution of conserved macroscopic densities of an underlying microscopic system [1, 2, 3]. It has been successfully used to model a wide range of high-energy systems spanning several disciplines in modern day physics. Both viscous and ideal hydrodynamics have been successfully used to describe the bulk dynamics of the quark-gluon plasma (QGP), an exotic state of matter produced in ultra-relativistic heavy-ion collisions [4, 5, 3, 6, 7]. Hydrodynamical theories are often used to model the matter content of neutron stars, with applications ranging from the stability of neutron stars to perturbations [8, 9, 10, 11, 12, 13, 14], to calculating the emission of gravitational waves from binary mergers involving neutron stars [15, 16, 17], with more recent studies beginning to investigate the potential importance of viscous effects on merger dynamics [18, 19, 20, 21, 22]. Further, relativistic hydrodynamics can also be used as a mirror to elucidate fundamental properties (e.g., turbulence) of purely gravitational systems through the fluid-gravity correspondence [23, 24, 25, 26, 27, 28, 29, 30], which establishes a duality between the Einstein equations in $(d+1)$-dimensional spacetimes with negative cosmological constant to the equations of relativistic hydrodynamics in $d$ dimensions. The ubiquity of relativistic hydrodynamics across such a wide variety of astrophysical and physical systems underscores the importance of its study as a stand-alone theory.

One approach to constructing theories of hydrodynamics is to assume the underlying flow is sufficiently close to local thermodynamic equilibrium such that it can be described using the typical thermodynamic variables (e.g., energy, pressure, temperature etc) plus gradients thereof, which are assumed to be small and are treated systematically order by order. Zeroth-order in this “gradient expansion” corresponds to local thermodynamic equilibrium in which the hydrodynamic variables are uniquely defined as the expectation value of underlying quantum operators. The zeroth-order hydrodynamical equations of motion—commonly referred to as the Euler equations—describe “ideal” or “perfect” fluid dynamics in which entropy is conserved (except in the presence of shocks). Since, for example, viscous dissipation and heat conduction arise due to velocity and temperature gradients, entropy-generating effects only enter the theory at first order, where one keeps terms linear in gradients of the hydrodynamic variables. A defining feature of first and higher-order theories is that the hydrodynamic variables are no longer uniquely defined—one has freedom in how to locally define the hydrodynamic variables without changing the physical content of the theory, so long as the same equilibrium state is approached in the limit of zero gradients [1, 2]. A given choice of definition of the hydrodynamic variables is typically referred to as a “choice of (hydrodynamic) frame.”

Constructing out-of-equilibrium theories is a highly non-trivial task because one’s choice of frame can profoundly impact the structure and mathematical properties of the resulting equations of motion. The early formulations of first-order dissipative theories resulting from the frame choices made by Eckart [31] and Landau and Lifschitz [32] were later shown to possess linearly unstable equilibrium states and acausal solutions [33, 34]. Second-order theories, such as that of Müller, Israel, and Stewart (MIS) [35, 36, 37, 38], became, and have largely since remained, popular alternatives to first-order approaches. MIS-type theories cure the pathologies of the first-order Eckart and Landau-Lifschitz theories by extending the state space to include additional, second-order degrees of freedom with their own evolution equations that are designed to guarantee non-negative entropy production [5]. Such second-order theories have been shown to possess linearly stable equilibrium states [39, 40]. At the nonlinear level, necessary conditions for causality with bulk and shear viscosities at vanishing chemical potential were derived in Ref. [41], nonlinear causality constraints with only energy or number diffusion were derived in Ref. [42], and, with only bulk viscosity present, causality and local well-posedness was established in Ref. [43]. These mathematical properties have motivated the use of MIS-type theories as the standard model of viscous relativistic fluids, having had particularly great success in accurately modeling various properties of the QGP in agreement with experimental data [3, 44].

Despite the success of second-order theories, there is sound motivation for seeking causal, stable, and locally well-posed generalizations of the first-order theories of Eckart and Landau and Lifschitz. The complicated partial differential equation (PDE) structure of the equations of motion in second-order theories has restricted rigorous proofs of causality, stability, and local well-posedness either to the linearized theory or when including only a subset of the full range of first and second-order transport coefficients. Further, even those proofs which do exist provide state-dependent, dynamic conditions that one must check are satisfied at each step in a numerical simulation, as opposed to state-independent conditions which can be satisfied a priori through an appropriate choice of frame. Further, it has been shown that MIS-type theories do not generically admit arbitrarily strong viscous shock solutions [40, 45]111In analogy to similar issues that arise in BDNK theory with a poor choice of hydrodynamic frame [46], we surmise this particular problem with the class of MIS theories analyzed in [40, 45] is a frame-induced hyperbolicity issue, where the maximum characteristic speed of the underlying PDEs is below the upstream flow velocity., and that solutions to such theories can develop singularities from smooth initial data in finite time [47].

In recent years, a promising potential alternative to second-order theories which addresses some of these limitations has emerged in the form of a first-order theory formulated by Bemfica, Disconzi, Noronha, and Kovtun (BDNK) over a series of works [48, 49, 50, 51, 2]. Owing to BDNK theory’s simpler PDE structure in comparison to second-order theories, causality, strong hyperbolicity, and local well-posedness in the full nonlinear regime and the linear stability of equilibrium states was rigorously proved in Ref. [50] with all physical first-order dissipative degrees of freedom present and without simplifying symmetry assumptions. However, to guarantee that these properties hold a suitable choice of hydrodynamic frame must be chosen, which amounts to satisfying a set of nonlinear constraints on the transport coefficients [50].

Following its formulation, there have been several analytic and numerical studies of solutions to the BDNK equations which have demonstrated its ability to reproduce the physical behavior expected of a dissipative relativistic fluid. Solutions to the BDNK equations in $4$D Minkowski spacetime for ($1+1$)D and ($2+1$)D conformal flows were obtained numerically in Refs. [46, 52], which considered various types of initial data ranging from the description of viscous shock waves to the Kelvin-Helmholtz instability. Complementing the numerical study in Ref. [46], the existence of arbitrarily strong shock wave solutions in suitably-chosen frames was rigorously proved in Ref. [53], with a further study of viscous BDNK shock wave solutions presented in Ref. [54]. A numerical analysis of solutions to the BDNK equations with non-trivial ideal gas microphysics was carried out in Ref. [55], while the relationship between BDNK theory and holography has also been explored mathematically in Refs. [56, 57] and numerically in Ref. [58].

In this paper, we present the first numerical solutions to the BDNK equations for a ($2+1$)D flow in spherical polar coordinates. In particular, we consider a fluid with an underlying conformal symmetry in $4$D Minkowski spacetime restricted by an unspecified radial force to a spherical shell of thickness $\delta$ and radius $R$. We take the shell to be infinitely thin ($\delta/R\ll 1$) and treat the flow as effectively constrained to the surface of a geometric sphere of radius $R$.

The main challenge of solving PDEs on the two-sphere is that one cannot cover the manifold $S^{2}$ with a single coordinate chart. Though convenient, spherical polar coordinates $(\theta,\phi)$ are singular at the poles $\theta\in\{0,\pi\}$. There are several ways of overcoming this difficulty—see, for example, Refs. [59, 60, 61, 29] for previous studies of relativistic hydrodynamics in curvilinear coordinates. In this work, we avoid the irregularity of spherical coordinates at the poles by instead adopting cubed-sphere coordinates [62, 63, 29], which cover $S^{2}$ continuously and regularly with six singularity-free coordinate charts. We numerically solve the equations of motion on the cubed-sphere multi-block grid using a fourth-order-accurate method of lines with a uniform finite-difference spatial discretization.

We consider three test cases of our code. First, we consider odd and even-parity perturbations of equilibrium states of the Euler and BDNK equations. We derive analytic predictions for the oscillation frequency and, for the BDNK fluid, characteristic damping time of these perturbations from the linearized equations of motion. We demonstrate that the analytic predictions agree with the oscillation and damping response present in the numerical solutions. We also consider Kelvin-Helmholtz-unstable initial data and qualitatively demonstrate how viscosity gives rise to vorticity diffusion and shears the characteristic rolls which form in the solutions to the Euler equations.

The third test case we consider is the evolution of a smooth, stationary Gaussian pulse in the energy density which is initially in local thermodynamic equilibrium. For small values of the entropy-normalized shear viscosity [$\eta/s\sim 1/(4\pi)$], viscosity has the expected effect of damping high-frequency components in the solution, and preventing formation of a shock that would otherwise form (i.e. the solution to the Euler equations beginning with the same initial conditions). Surprisingly, for sufficiently large values of the $\eta/s$ [$\sim 10/(4\pi)$], viscosity stops “working” as expected, and very steep gradients form in the numerical solution, causing the flow to diverge away from local equilibrium and the regime of validity of first-order hydrodynamics before convergence is lost in the numerical simulations at all grid resolutions we have considered.

We observe similar divergent behavior in ($1+1$)D numerical simulations of a planar-symmetric conformal flow in Cartesian coordinates [with $\eta/s=20/(4\pi)$], as presented in App. A. In this case, the reduced dimensionality of the flow makes achieving much higher grid resolutions computationally feasible, allowing us to provide a more detailed analysis of how this behavior is persistent at very high resolutions and hence is most likely not a numerical artifact. Rather, the trend in the behavior of the convergence tests with increasing grid resolution suggests a discontinuity forms in finite time in the corresponding continuum solution, providing numerical evidence that singularities can develop in solutions to the BDNK equations from smooth initial data. Though achieving such high resolutions on the two-sphere is too computationally expensive, the similar qualitative behavior to the ($1+1$)D solutions suggests that discontinuities likely also form in the corresponding ($2+1$)D continuum solutions.

We emphasize that this work focuses on the mathematical well-posedness of the BDNK equations and their numerical integration on a spherical manifold. We refer the reader to, for example, Refs. [4, 6, 7, 64], Refs. [15, 16, 17, 65, 66] and Refs. [67, 68, 69, 70], respectively, for extensive discussions on the application of relativistic hydrodynamics to modeling physical systems such as the QGP, neutron stars and general-relativistic magnetohydrodynamical flows.

The remainder of this paper is organized as follows. We begin by reviewing first-order hydrodynamics and BDNK theory in Sec. II. We describe the tensor components, equations of motion, and our numerical methods on the two-sphere in Sec. III. Our results are presented in Sec. IV, with the study of linear fluid perturbations presented in Sec. IV.1, Gaussian initial data in Sec. IV.2, and Kelvin-Helmholtz-unstable initial data in Sec. IV.3. We conclude with a summary and discussion of our results in Sec. V. We relegate various technical details to the appendix. We present in App. A our ($1+1$)D numerical simulations. We relegate lengthy expressions for the coordinate BDNK tensor components and equations of motion to App. B. We consider in App. C general spherical harmonic perturbations of equilibrium states, deriving dispersion relations for the harmonic frequency for odd and even-parity perturbations and demonstrating linear mode stability for a class of hydrodynamic frames in the high-frequency ($l\to\infty$) limit. We relegate additional details about cubed-sphere coordinates to App. D and to App. E convergence tests for the numerical simulations presented in this work.

We work in natural units where $\hbar=c=k_{B}=1$ and report all dimensionful quantities in code units. We use the metric signature $(-,+,+,+)$.

II First-Order Relativistic Hydrodynamics

In this section, we first review the gradient expansion approach to dissipative relativistic hydrodynamics. We then provide a brief review of BDNK theory and discuss our choice of hydrodynamic frame and the regime-of-validity diagnostics we keep track of in our numerical simulations.

II.1 The Gradient Expansion

Hydrodynamics is a long-wavelength effective theory describing the evolution of conserved, macroscopic densities. These densities are ultimately determined by the underlying microscopic theory describing a given system, and are obtained by averaging out the microphysical degrees of freedom with length scales much shorter than that of interest in the effective theory. In relativistic hydrodynamics, these conserved densities are the stress-energy tensor $T^{\mu\nu}$ and the baryon current tensor $J^{\mu}$, which, given a suitable set of initial data, are evolved forward in time by the conservation equations

[TABLE]

The underlying microphysics informs the macroscopic content in the conserved densities through an additional relation between state variables required to close the system (1) and the dissipative degrees of freedom which arise outside of thermodynamic equilibrium.

To construct a theory of relativistic hydrodynamics, one typically begins by decomposing the conserved densities transverse and longitudinal with respect to a timelike four-velocity $u^{\mu}$  [31, 1, 2]:

[TABLE]

where, in $d+1$ dimensions,

[TABLE]

and $\Delta^{\alpha\beta}=g^{\alpha\beta}+u^{\alpha}u^{\beta}$ is a projector onto the space orthogonal to $u^{\alpha}$. The four-vectors $\mathcal{Q}^{\mu}$, $\mathcal{J}^{\mu}$ are transverse to $u^{\mu}$ and the rank-two tensor $\mathcal{T}^{\mu\nu}$ is transverse, symmetric, and traceless.

The decomposition (2) is general. The gradient expansion approach to hydrodynamics assumes that, for a flow sufficiently close to local thermodynamic equilibrium, the scalars and tensors (3) can be expressed in terms of a chosen set of hydrodynamic variables and their gradients, forming a set of so-called constitutive relations. This expansion is carried out order by order,

[TABLE]

where $T^{\mu\nu}_{(k)}$ and $J^{\mu}_{(k)}$ contain $\mathcal{O}(\partial^{k})$ terms, products of $\mathcal{O}(\partial^{k-1})$ and $\mathcal{O}(\partial)$ terms, and so on. Truncating the gradient expansion at some order $k$ and specifying a set of constitutive relations gives rise to $k$th-order hydrodynamics.

At zeroth order in the gradient expansion, one obtains the perfect-fluid densities,

[TABLE]

where $\epsilon$, $P$, and $n$ denote the energy density, pressure, and baryon number density, respectively, that an observer co-moving with the fluid would measure. The hydrodynamic equations of motion (1) in this case are commonly referred to as the Euler equations, which constitute a locally well-posed initial value problem [71, 5] and which exactly conserve entropy (in the absence of shocks). From Eqs. (1, 5), one can see that in four spacetime dimensions, there are only five evolution equations for the six variables $\{\epsilon,p,n,u^{\mu}\}$. To close the system, an additional relation $P=P(\epsilon,n)$ between the hydrodynamic variables is required, and is typically determined from the underlying theory which describes the microphysics of the system.

Entropy-generating effects, such as heat flow due to thermal gradients and the viscous response of fluid elements to shear and bulk stresses, enter the theory at first order in the gradient expansion. However, several complications of the theory also arise outside of local thermodynamic equilibrium. In equilibrium, the hydrodynamic variables are uniquely defined as expectation values of underlying microscopic operators; out-of-equilibrium, the hydrodynamic variables no longer have unique definitions which follow from microscopic theory—sufficiently close to equilibrium, they can be viewed merely as parametrizations of the densities $T^{\mu\nu}$ and $J^{\mu}$ which do retain their unique microscopic definitions [3, 1, 2].

At first order, the most general definition of the out-of-equilibrium variables consistent with thermodynamics is parametrized by fourteen transport coefficients, only six independent combinations of which are invariant under field redefinitions that leave $T^{\mu\nu}$ and $J^{\mu}$ invariant (to first order) [2]. Of these six, three correspond to the “classical” transport coefficients of shear viscosity, bulk viscosity and charge conductivity ($\eta$, $\zeta$ and $\sigma$, respectively), while, in BDNK theory, the remaining three constitute relaxation times which act as causal regulators. Since these dissipative processes ultimately arise from microphysical interactions, one must again resort to the underlying microscopic theory to determine the values of the transport coefficients for a given system of interest.

One’s choice of definition of the out-of-equilibrium hydrodynamic variables (also referred to as a choice of “hydrodynamic frame”) plays a crucial role in the mathematical properties of the resulting equations of motion. The first formulations of first-order viscous theories were proposed by Eckart [31] and Landau and Lifschitz [32], who proposed frame choices in which a co-moving observer measures no out-of-equilibrium corrections to the energy density ($\mathcal{E}=\epsilon$) and, respectively, no out-of-equilibrium charge flow ($\mathcal{J}^{\mu}=0$) or no heat flux ($\mathcal{Q}^{\mu}=0$). Despite making these frame choices which have some physical consequence, these theories were later shown to violate causality and possess unstable equilibrium states [33, 34]. The first causal, stable and locally well-posed first-order theory of viscous hydrodynamics (without assumptions on certain transport coefficients vanishing or restrictively smooth initial data) was formulated by Bemfica, Disconzi, Noronha and Kovtun (BDNK) over a series of works [48, 49, 50, 51, 2]. In the next section, we present a review of BDNK theory and specialize to the case of a fluid with an underlying conformal symmetry.

II.2 BDNK Theory

The BDNK tensor for a non-conformal fluid with all six transport coefficients present was introduced in Ref. [50], in which causality, linear stability of equilibrium states, local well-posedness, and strong hyperbolicity when coupled to the Einstein equations were proven over a range of hydrodynamic frames. In BDNK theory, all out-of-equilibrium corrections to the baryon current density are set to zero (i.e., $\mathcal{N}=n$ and $\mathcal{J}^{\mu}=0$ in Eq. 2) and the hydrodynamic variables $\epsilon$ and $n$ are used to parametrize the conserved densities. In particular, the BDNK tensors are given by [50]

[TABLE]

where the out-of-equilibrium corrections to the energy density and pressure are given by

[TABLE]

and the heat-flow vector and shear tensor take the forms

[TABLE]

where the $\beta_{i}$ are given by

[TABLE]

In this work, we restrict to a fluid with an underlying conformal symmetry (which implies $T^{\mu}{}_{\mu}=0$) and zero chemical potential. These choices imply [52]

[TABLE]

so that only $\eta$, $\tau_{\epsilon}$ and $\tau_{Q}$ remain as non-trivial transport coefficients. Defining the related quantities $\chi$ and $\lambda$ by

[TABLE]

causality, stability and local well-posedness of conformal BDNK theory is guaranteed in frames which satisfy [48]

[TABLE]

Following Refs. [46, 52], we choose to utilize the dimensional relations for a conformal fluid $\epsilon\propto T^{4}$ and $\{\eta,\chi,\lambda\}\propto T^{3}$ to parametrize the transport coefficients in terms of the energy density:

[TABLE]

where $\lambda_{0}$, $\chi_{0}$, and $\eta_{0}$ are dimensionless constants. We work exclusively in the class of frames (frame B in Ref. [46])

[TABLE]

which satisfies (12) and has a maximum characteristic speed equal to the speed of light [46].

Throughout this work, we will specify the remaining degree of freedom in our choice of frame by either specifying $\eta_{0}$ directly, or implicitly through a choice of the entropy-normalized shear viscosity $\eta/s$, following the conventions of Ref. [46]. Using the relations $s=(\epsilon+P)/T$ and $\epsilon=\epsilon_{0}T^{4}$ for conformal fluids [46, 48], one obtains

[TABLE]

where, as in Ref. [46], we take $\varepsilon_{0}=10$, which is motivated by $\varepsilon_{0}=\epsilon/T^{4}\sim{10}$ for the quark-gluon plasma [72, 73].

Given Eqs. (14, 15), a choice of $\eta_{0}$ or $\eta/s$ controls how much dissipation there is in the fluid. In particular, taking $\eta_{0}\to 0$ recovers the perfect fluid densities and Euler equations. Throughout, we consider values of $\eta/s$ which are multiples of the Kovtun-Son-Starinets (KSS) bound: a minimum value of $\eta/s$ derived using string-theory methods for a wide class of strongly-interacting quantum field theories and conjectured to hold for all fluids [74],

[TABLE]

where we have restored factors of $\hbar$ and $k_{B}$. Experimental measurements of $\eta/s$ for several strongly-interacting systems have respected the KSS bound: for finite nuclear matter, measurements have yielded $\eta/s\sim(2.5-6.5)\cdot(\eta/s)_{\mathrm{min}}$ [75]; for the QGP, $\eta/s\sim 0.1$ (which is $\sim 25\%$ larger than the KSS bound) [76, 77]. For standard laboratory fluids, $\eta/s$ is typically much larger than the KSS bound, e.g., water at STP has $\eta/s\sim 380\cdot(\eta/s)_{\mathrm{min}}$ [74] (see, for example, Ref. [76] for additional examples). In this work, we consider values of $\eta/s$ ranging from $(1-20)\cdot(\eta/s)_{\mathrm{min}}$.

II.3 Regime of Validity Diagnostics

When truncating the gradient expansion (4) at some finite order, one assumes that the higher-order corrections are negligible. Such an assumption is not guaranteed to hold at all times during an evolution—one may have constructed far-from-equilibrium initial data or there could be dynamic evolution towards a regime where higher-order viscous corrections become important. In a gradient expansion approach, it is thus important to keep track of diagnostic quantities during an evolution that provide information about whether the flow is instantaneously within the regime of validity of the underlying viscous theory.

We keep track of two diagnostics in our BDNK simulations: the inner product $u_{\mu}u_{\nu}T^{\mu\nu}$ to monitor the violation of the weak-energy condition as suggested in Ref. [48], and the relative size of the first-order corrections to the BDNK tensor components via the pointwise ratio $|{T^{\mu\nu}_{(1)}}|/|T^{\mu\nu}_{(0)}|$ as done in Ref. [46]. Since $u_{\mu}u_{\nu}T^{\mu\nu}=\epsilon+\mathcal{A}$ for the BDNK tensor (6a), a severe violation of the weak-energy condition indicates the first-order dissipative correction $\mathcal{A}$ may be sufficiently large such that the flow has exited the regime of validity of the first-order theory. Values $|{T^{\mu\nu}_{(1)}}|/|T^{\mu\nu}_{(0)}|\gtrsim{1}$ indicate the first-order viscous terms dominate the ideal ones and are not small corrections as assumed when truncating the gradient expansion at first order, and that higher-order terms also likely constitute sizable corrections.

Even if such diagnostic quantities strongly suggest the flow has exited the regime of validity of the theory, one still has valid solutions to the given equations of motion. However, in order to accurately model the physical evolution of the system under consideration, the equations of motion must be replaced with those from an alternative out-of-equilibrium theory. If the gradient expansion still converges, then a higher-order theory may accurately model the flow, otherwise a non-perturbative approach would be necessary. In this work, we use the diagnostic quantities to help interpret the out-of-equilibrium behavior of the numerical evolution of Gaussian initial data discussed in Sec. IV.2 (and in App. A).

III Equations of Motion and Numerical Methods on the Two-Sphere

In this section, we describe the hydrodynamical equations of motion, numerical methods, and analysis tools we use on the two-sphere. In Sec. III.1, we describe the background geometry of the flow, equations of motion and choices we make regarding initial data. We describe our numerical methods in Sec. III.2 and define some useful analysis tools in Sec. III.3 that we make use of in later parts of this work.

III.1 Equations of Motion

We consider a conformal fluid in $4$D Minkowski spacetime constrained to a geometric sphere of radius $r=R$. In spherical coordinates, the line element, $\mathop{}\!\mathrm{d}s^{2}=\eta_{\mu\nu}\mathop{}\!\mathrm{d}x^{\mu}\mathop{}\!\mathrm{d}x^{\nu}$, is given by

[TABLE]

and the four-velocity of the flow has components

[TABLE]

where $u^{t}$, $u^{\theta}$ and $u^{\phi}$ are functions of $(t,\theta,\phi)$. We take as independent components $u^{\theta}$ and $u^{\phi}$, with $u^{t}$ then determined by the normalization $u_{\mu}u^{\mu}=-1$,

[TABLE]

The perfect fluid (5) and BDNK (6) densities can be expressed in terms of the fields $\{\epsilon,n,u^{\theta},u^{\phi}\}$, which, upon substitution into the equations of motion (1), give a coupled set of PDEs that we evolve forward in time from a given set of initial data. We relegate to App. B.2 the (lengthy) expressions for the tensor components and equations of motion for the Euler and BDNK systems.

For a conformal fluid, $n$ does not feature in the BDNK stress-energy tensor (since $\beta_{n}=0$) and the first-order PDE governing the evolution of $n$ consequently decouples from the rest of the system. We are thus free to consider flows both with or without a baryon number density $n$. For the most part we ignore $n$, though in some cases it is useful to evolve $n$ as a tracer of the underlying fluid dynamics.

The remaining fields $\{\epsilon,u^{\theta},u^{\phi}\}$ are evolved forward in time from their specified initial data using either the Euler equations or the BDNK equations. Evolutions using the BDNK equations, which are second order PDEs, additionally require initial data for the time derivatives of the fields. We construct this initial data implicitly by setting the out-of-equilibrium terms $T_{(1)}^{tt}=T_{(1)}^{t\theta}=T_{(1)}^{t\phi}=0$ at $t=0$, which then determines the initial time derivatives of the hydrodynamic fields. Since the out-of-equilibrium stress-energy tensor is uniquely defined, unlike the hydrodynamic variables, this is the natural choice which ensures the flow is initialized in equilibrium. When considering time-dependent fluid perturbations in Sec. IV.1, we will instead choose to directly specify the initial time derivative of the hydrodynamic fields instead of setting the out-of-equilibrium $T^{t\mu}_{(1)}$ components to zero.

III.2 Numerical Methods

One of the main challenges of solving PDEs on the two-sphere is that it cannot be covered by a single coordinate chart that is regular everywhere. We circumvent this complication by solving the equations of motion on a multi-block grid in cubed-sphere coordinates [62, 63, 29]. The essence of the cubed-sphere approach is to cover the surface of the sphere with six non-overlapping coordinate patches, each of which is the projection of a face of a concentrically-placed cube onto the sphere. The cube can then be “unfolded”, giving rise to a regular, singularity-free grid structure on which the underlying equations can be discretized and solved. We show in Fig. 1 a geometric illustration of cubed-sphere coordinates and relegate details about coordinate transformations and the metric on the cubed-sphere grid to App. D.

There are two main challenges in solving hyperbolic PDEs on a multi-block grid structure. First, one must allow propagation of information across the distinct coordinate patches. Second, if one employs a uniform grid structure, the points along the boundary of each sub-grid are shared with the neighboring sub-grids, so one must also ensure the solution is not multi-valued at these points. As described in Ref. [63], one approach to solving these problems is to first decouple the six sub-grids by employing off-centered finite difference stencils near the boundaries, but then reintroduce the coupling through the use of so-called penalty terms. The latter are designed to drive the solution at the shared boundary points to a common value, and allow information to propagate between sub-grids. Here, we employ a different approach in which we use centered finite difference stencils at every point in the computational domain. When solving the discretized equations at or near the boundary of a given sub-grid such that a finite difference stencil extends off the grid, we map the off-grid point(s) in the stencil back onto the sphere and back again onto the cubed-sphere grid such that it falls inside the computational domain of a neighboring sub-grid. We then evaluate a cubic-spline interpolant of the solution in this adjacent sub-grid to set the corresponding variable values in the finite difference stencil, performing the necessary Jacobian transformations for the four-velocity components.

With the aforementioned uniform grid structure and treatment of boundaries, we numerically solve the Euler and BDNK equations using a fourth-order-accurate method of lines, augmented with Kreiss-Oliger dissipation [78]. That is, we first discretize the equations in space using fourth-order-accurate finite difference stencils, then perform a time integration using a fourth-order-accurate Runge-Kutta scheme; Kreiss-Oliger dissipation is added as a filtering operation between time steps.

III.3 Analysis Tools

III.3.1 Spherical Harmonic Expansion of Functions on the Two-sphere

The scalar spherical harmonics form an orthonormal basis for scalar functions on the surface of a sphere. Any square-integrable scalar function $f(\theta,\phi)$ on the sphere can be expanded in the form

[TABLE]

where we consider real-valued functions $f:S^{2}\to\mathbb{R}$ and real spherical harmonics $Y_{l}^{m}(\theta,\phi)$ with normalization

[TABLE]

The total power of $f$ on the unit sphere is related to its spectral coefficients $A_{lm}$ by

[TABLE]

where

[TABLE]

The coefficients $C_{l}$ contain information about the frequency-domain structure of the function $f$. We compute $C_{l}$ for $l\in\{0,\ldots,12\}$ from the spectral coefficients $A_{lm}$, which we evaluate by making use of the orthogonality relation (21) applied to (20). The resulting integral for the $A_{lm}$ is numerically evaluated in spherical coordinates by, for each function call of $f$, transforming the coordinates $(\theta,\phi)$ to cubed-sphere coordinates and evaluating a cubic-spline interpolant of the scalar quantity $f$ from our numerical solution on the cubed-sphere grid. In all figures, we plot the $C_{l}$ normalized by the total power of $f$ as given by similarly numerically evaluating the left-hand side of (22).

III.3.2 Vorticity

An important quantity in studying flows is the vorticity, given by the curl of the flow velocity for non-relativistic fluids. For a relativistic fluid, the vorticity two-form $\omega_{\mu\nu}$ is given by the exterior derivative of the flow velocity,

[TABLE]

For a ($2+1$)D inviscid flow (where in the remainder of this subsection we use Greek indices to range over the time ([math]) and two space $(\theta,\phi)$ coordinates), a conserved quantity associated with $\omega_{\mu\nu}$ is given by the timelike component of its Hodge dual, $W^{\alpha}\equiv\epsilon^{\alpha\mu\nu}\omega_{\mu\nu}$, where

[TABLE]

is conserved [29]. Here, $\epsilon^{\alpha\mu\nu}$ denotes the Levi-Civita tensor (the Levi-Civita symbol normalized by the square root of the determinant of the metric), which, in combination with (24), gives the following expression for the conserved quantity $W^{0}$:

[TABLE]

In subsequent sections, we will refer to $W^{0}$ as the vorticity density.

IV Results

In this section, we present numerical solutions to the Euler and BDNK equations on the two-sphere for three types of initial data. In Sec. IV.1, we compare the inviscid and viscous numerical evolution of perturbed equilibrium states to the analytic oscillation and damping response predicted by the linearized equations of motion. We present in Sec. IV.2 the evolution of smooth, stationary Gaussian initial data and in Sec. IV.3 the evolution of a Kelvin-Helmholtz-unstable equatorial jet of over-density.

In Tab. 1 we list the maximum number of grid points, the Courant–Friedrichs–Lewy (CFL) factor used for time stepping and the Kreiss-Oliger dissipation coefficient used in our numerical simulations. Convergence tests for all the simulations presented in this work can be found in App. E.

IV.1 Linear Normal Mode Analysis

As our first test case, we consider the oscillation and, for the BDNK fluid, damping response to perturbations of uniform equilibrium states of the form $\epsilon=\mathrm{constant}$ and $u^{\theta}=u^{\phi}=0$. In Sec. IV.1.1, we explicitly write down the form of the fluid perturbations we consider using scalar and vector spherical harmonics. In Sec. IV.1.2, we extract from the linearized Euler and BDNK equations the oscillation frequency and (BDNK) damping timescale of the perturbations, which we compare to numerical solutions of the full nonlinear equations in Sec. IV.1.3.

In this section, we specialize to $(l=1,m=0)$ spherical harmonic perturbations. We present in App. C a more comprehensive analysis for general $(l,m)$ and demonstrate linear mode stability in an explicit class of hydrodynamics frames for odd and, in the high-frequency $l\to\infty$ limit, even-parity perturbations.

IV.1.1 Perturbed Fluid States

We perturb the energy density with a scalar spherical harmonic,

[TABLE]

where $\epsilon_{0}$ is the constant background energy density and $\delta_{\epsilon}\ll\epsilon_{0}$. We perturb the four-velocity components $u^{\theta}$ and $u^{\phi}$ using vector spherical harmonics. We adopt the conventions of Ref. [79], in which two types of vector spherical harmonics with different parities are defined:

[TABLE]

where $i,j\in\{2,3\}$, $\epsilon_{2}{}^{2}=\epsilon_{3}{}^{3}=0$, $\epsilon_{2}{}^{3}=-\csc{\theta}$ and $\epsilon_{3}{}^{2}=\sin{\theta}$. Under inversion about the origin, $(\theta,\phi)\to(\pi-\theta,\phi+\pi)$, the scalar spherical harmonics obey $Y_{l}{}^{m}\to(-1)^{l}Y_{l}{}^{m}$, i.e., they have parity $(-1)^{l}$. The vector spherical harmonics $V_{l}{}^{m}$ and $W_{l}{}^{m}$ have parity $(-1)^{l}$ and $(-1)^{l+1}$, respectively [79]. We refer to $V_{l}{}^{m}$ and $W_{l}{}^{m}$, and their respective vector perturbations, as having even and odd parity, respectively.

We consider $l=1$, $m=0$ fluid perturbations, for which the corresponding scalar spherical harmonic is given by

[TABLE]

and the non-trivial contravariant components of the vector spherical harmonics are

[TABLE]

The even and odd-perturbed four-velocity components thus assume the form

[TABLE]

where $\delta_{\theta,\phi}\ll 1$ and, in the even case, normalization of the four-velocity gives

[TABLE]

and similarly for the odd perturbations.

IV.1.2 Linearized PDE Analysis

Throughout this section, we state the $(l=1,m=0)$ instances of the general $(l,m)$ equations of motion derived in App. C.

The non-trivial linearized equations of motion governing odd perturbations of the Euler fluid reduce to

[TABLE]

where we have taken $\epsilon_{0}>0$. Equation (33) is solved by either $\delta_{\phi}\neq 0$ with $\omega=0$ (i.e., rigid rotation), or the trivial solution $\delta_{\phi}=0$. Linearized odd perturbations of the BDNK fluid obey

[TABLE]

Equation (34) shares the same trivial solutions as the odd-perturbed Euler equations (33), in addition to the solution

[TABLE]

Since $\rm{Re}(\omega)=0$ and $\rm{Im}(\omega)<0$ in (35), this solution corresponds to pure damping of the odd perturbation (31b) with a characteristic time scale $\tau_{\mathrm{d}}=3\lambda_{0}\epsilon_{0}^{-1/4}/4$. However, since the solution (35) has $\delta_{\epsilon}=0$ and $W^{\phi}$ is constant on the two-sphere, the underlying perturbation introduces only rigid rotation (i.e., it contains no gradients) and we therefore do not expect any physical damping to ensue. This expectation is reconciled with the damping of $u^{\phi}$ because the solution (35) corresponds to a frame mode: the non-uniquely-defined out-of-equilibrium component $u^{\phi}$ is damped, but this is purely a result of the choice of frame and is not a physical effect. This is evident from the stress-energy tensor components of the perturbed fluid state corresponding to the solution (35):

[TABLE]

Absent from (36) (and, in particular, the component $T^{t\phi}=0$) are any damping terms, meaning that there is no physical damping of the energy density or momentum.

The linearized equations of motion governing the even perturbations of the Euler fluid are

[TABLE]

which have solutions

[TABLE]

Linearized even perturbations of the BDNK fluid obey

[TABLE]

Fixing $\eta/s=(4\pi)^{-1}$, $R=1$, and $\epsilon_{0}=1$, Eq. (39) has four solutions (reported to two decimal places)

[TABLE]

where $k\in\{0,1\}$.

The six solutions in Eqs. (38, 40) correspond to only three distinct sets of initial data, which are given by the real parts of Eqs. (27, 31a) and their time derivatives. In the next section, we compare the oscillations and damping in the numerical evolution of this initial data to the analytic predictions from the linear analysis of even perturbations to the Euler and BDNK equations.

IV.1.3 Numerical Simulations

In addition to the choices $\eta/s=(4\pi)^{-1}$, $R=1$, and $\epsilon_{0}=1$, we also take $\delta_{\epsilon}/\epsilon_{0}=0.01$. We refer to the numerical evolution of the initial data corresponding to the BDNK solutions in Eqs. (40a, 40b) as simulations I and II, respectively.

We estimate oscillation frequencies and damping times in the numerical solutions from the time series of the energy density at the north pole, which we plot in Fig. 2 for all three simulations. As expected in the inviscid case, the flow is purely oscillatory with the peak in the energy density migrating from pole to pole. As also manifest in Fig. 2, BDNK simulations I and II undergo under and overdamped oscillations, respectively. We fit the Euler time series to a sinusoidal functional form and the BDNK time series to an exponentially damped sinusoid to estimate the oscillation frequencies and damping times, finding good agreement with the corresponding predicted values from the linearized equations of motion. This is illustrated in Tab. 2, wherein we list the percentage (absolute) residuals between the predicted and fitted oscillation frequencies and damping times as a function of grid resolution.

IV.2 Gaussian Initial Data

In this section, we study the evolution of smooth, stationary Gaussian initial data for increasing values of $\eta/s$. In particular, we consider the following initial data:

[TABLE]

where $\epsilon_{0}=0.1$, $A=0.4$, $w=45^{\circ}$. As discussed in Sec. III.1, the initial data {$\dot{\epsilon},\dot{u}^{\theta},\dot{u}^{\phi}$} for the BDNK equations is set such that the out-of-equilibrium components $T^{tt}_{(1)}$, $T^{t\theta}_{(1)}$, and $T^{t\phi}_{(1)}$ are initially zero, both for this set of initial data and that considered in the following section.

The initial data (41), shown in the upper left panel of Fig. 3, corresponds to an initially-stationary energy density profile with a peak centered on the north pole. The subsequent dynamics of the inviscid [$\eta/s=0/(4\pi)$] flow consists of the energy density spreading out into an axisymmetric azimuthal band-like structure as it migrates towards the south pole in the polar direction, where it once again concentrates into a narrow peak before returning to the north pole. We plot in the upper right of Fig. 3 a snapshot of the energy density just before it reaches the south pole for the first time.

Each pole-to-pole traversal is accompanied by the formation of increasingly steep gradients and a narrower, more ring-like distribution of the energy density. This is illustrated in the snapshots of the solution shown in the lower row of Fig. 3 just before the energy density peaks for the second time at the north pole. Thereafter, the energy density is focused into an extremely narrow peak, at which time a shock forms and convergence is lost in our numerical simulations. We plot in the left panel of Fig. 4 the energy density at the time when convergence is lost, wherein one can see that, relative to the initial configuration (upper left panel of Fig. 3), the energy density has been effectively focused into a point.

As expected, viscosity serves to smooth the shock that forms in the inviscid solution. We plot in the right panel of Fig. 4 the energy density with $\eta/s=1/(4\pi)$ at the same time convergence is lost for $\eta/s=0/(4\pi)$. There, we see that the peak in the energy density is much broader in the viscous compared to the inviscid solution, with the former’s amplitude being $\sim 0.5$ times that at time $t=0$, compared to $\sim 1.5$ for the latter (see the colorbars in Figs. 3-4). Thereafter, the $\eta/s=1/(4\pi)$ solution is further damped until a steady state is reached, with convergence being maintained throughout the simulation.

The behavior of the BDNK fluid differs from the pure damping effect observed in the $\eta/s=1/(4\pi)$ solution for larger values of $\eta/s$, for which increasingly steep gradients form in the solution. We illustrate this effect qualitatively in Fig. 5, wherein we show a mode decomposition of the energy density for $\eta/s\in\{0,1,3,10\}/(4\pi)$. Relative to the inviscid flow, the viscosity in the $\eta/s=1/(4\pi)$ case significantly damps higher frequency modes in the solution, thereby preventing the formation of a shock. As $\eta/s$ is further increased, the relative magnitude of the higher-frequency components also increases. Though the $\eta/s\in\{2,3\}/(4\pi)$ solutions remain qualitatively similar to the $\eta/s=1/(4\pi)$ case, we require higher grid resolutions to maintain fourth-order convergence in the numerical simulations.

The behavior of the $\eta/s\in\{10,20\}/(4\pi)$ solutions, however, markedly differ from that of the less viscous flows. Sufficiently steep gradients form in these more viscous flows that numerical convergence is lost at all grid resolutions we consider, which is to be expected if a discontinuity forms in the underlying continuum solution, since our numerical method assumes the solution is smooth. The onset of these steep gradients coincides with the diagnostics we monitor indicating the flow has dynamically diverged away from the regime of validity of first-order hydrodynamics.

Unlike in the $\eta/s\in\{0,1\}/(4\pi)$ solutions, a secondary pulse forms in the larger $\eta/s$ solutions which approaches the south pole as the main peak in the energy density returns to the north pole. Though the trailing part of the secondary peak remains spread out, the front of the secondary pulse develops increasingly sharp features and becomes larger in amplitude (relative to the primary pulse) as $\eta/s$ is increased. In the vicinity of the secondary pulse’s traveling front, the gradients in the solution grow steeper with $\eta/s$, causing the viscous contributions in the stress-energy tensor to increasingly dominate the inviscid contributions, resulting in the the regime-of-validity diagnostics assuming more extreme values.

We plot in Fig. 6 the $\eta/s=20/(4\pi)$ profile and regime-of-validity diagnostics when convergence is lost in the numerical simulation. At this time, the fronts of the secondary pulse have developed a sharp, ring-like structure, in the vicinity of which the regime-of-validity diagnostics take on extreme values. As shown in the lower row of Fig. 6, $|T^{tt}_{(1)}|/|T^{tt}_{(0)}|\sim 20$ and $u_{\mu}u_{\nu}T^{\mu\nu}\sim-3$ at the location of these fronts when convergence is lost, strongly indicating that the regime of validity of BDNK theory has been exited. The behavior of the regime-of-validity diagnostics for different values of $\eta/s$ is illustrated in Fig. 7, wherein we show the maximum and minimum values of $|T^{tt}_{(1)}|/|T^{tt}_{(0)}|$ and $u_{\mu}u_{\nu}T^{\mu\nu}$, respectively, over the spatial domain as a function of time. For $\eta/s\in\{10,20\}/(4\pi)$, the diagnostics rapidly take on increasingly extreme values as convergence is lost at $t\approx 6$.

The viscosity in the numerical simulations presented in this section thus appears to operate within two regimes. As one “turns on” the viscosity, the high-frequency components in the inviscid solution are damped, preventing the formation of a shock. As $\eta/s$ is increased to sufficiently large values, the flow rapidly enters an out-of-equilibrium state, and the formation of steep gradients prevents a hydrodynamic-attractor-like return to local equilibrium at late times. Instead, at all grid resolutions we have considered (up to a total of $\sim 10^{14}$ grid points—see Tab. 1), the formation of steep gradients is accompanied by a loss of convergence, and the regime-of-validity diagnostics assuming extreme values as the flow diverges away from equilibrium.

We observe similar qualitative behavior in numerical simulations of planar-symmetric ($1+1$)D conformal flows in Minkowski spacetime with smooth, stationary Gaussian initial data, as presented in App. A. For $\eta/s=1/(4\pi)$, the viscosity smooths a shock that forms in the inviscid flow and we converge to a steady state at late times. We converge to a late-time steady state for simulations with up to $\eta/s=10/(4\pi)$. For $\eta/s=20/(4\pi)$ we observe a similar steepening of gradients until convergence is lost and the monitored diagnostics assume extreme values which strongly indicate the flow exits the regime of validity of first-order hydrodynamics. The loss of numerical convergence and the blow up of the diagnostics are related insofar as the diagnostics indicate the regime of validity of BDNK theory is exited before the loss of convergence as very steep gradients form, and since these steep gradients appear to evolve towards a discontinuity, the loss of convergence is inevitable since our numerical method is not designed to capture shocks. However, all the evidence from our numerical simulations is that there is no causal relationship between the loss of numerical convergence and the blow up of the diagnostics. Since the ($1+1$)D simulations are less computationally expensive than the ($2+1$)D ones, we are able to achieve much higher grid resolutions, allowing us to present a more comprehensive analysis which suggests that the discontinuities are not numerical artifacts that can be resolved with sufficiently high resolution, nor are they a result of numerical instabilities. Rather, with increasing grid resolution, we maintain convergence incrementally closer to some finite time at which discontinuities appear to form and numerical convergence is lost, indicating a break down in smoothness in the continuum solution that our numerical method is, by design, unable to capture, thereby resulting in the loss of numerical convergence.

Taken together, the simulations presented in this section and in App. A provide numerical evidence that discontinuities can develop in finite time in continuum solutions to the BDNK equations with smooth initial data.

IV.3 Kelvin-Helmholtz Instability

Equilibrium configurations consisting of a continuous inviscid fluid with a velocity shear, or two inviscid fluids with a shear in the flow velocity parallel to their separating interface are unstable to perturbations. This class of instabilities is generally referred to as the Kelvin-Helmholtz instability. In non-relativistic flows, gravity, surface tension effects and viscosity can act as regulators of the instability and prevent it from growing without bound. In this section, we consider the effect of viscosity on the evolution of Kelvin-Helmholtz-unstable initial data by qualitatively comparing its evolution by the Euler and BDNK equations.

A comparison between the evolution of Kelvin-Helmholtz-unstable initial data by the Euler and BDNK equations for a $4$D conformal fluid with variations in two dimensions in Cartesian coordinates and Minkowski spacetime was carried out in Ref. [52], explicitly demonstrating the shearing of the characteristic Kelvin-Helmholtz rolls by the viscosity in the BDNK equations. Here, we demonstrate the same shearing effect in addition to the diffusion of vorticity in the viscous evolution of Kelvin-Helmholtz-unstable initial data on the two-sphere.

We consider the following initial data adapted from Refs. [52, 80]:

[TABLE]

where $c_{s}=1/\sqrt{3}$, $a=0.08$, $\sigma=0.2$, $\theta_{1}=70^{\circ}$, $\theta_{2}=110^{\circ}$, and we consider the two values $A\in\{0,0.01\}$.

The initial data (42) corresponds to a jet of number density $n=2$ and angular velocity $c_{s}/4$ centered on the equator with two interfaces which border regions with number density $n=1$ and angular velocity $-c_{s}/4$—the parameter $a$ in (42) sets the sharpness of the transition in $n$ and $u^{\phi}$ across the interfaces. The $A=0.01$ case introduces a small, sinusoidally varying component of the polar velocity along each of the interfaces which seeds the instability by both speeding up the dynamics and causing the jet to deform as the Kelvin-Helmholtz rolls begin to form. The values of $\eta/s\sim(4\pi)^{-1}$ we have considered thus far are sufficiently large to completely prevent the onset of the instability. Since we must consider sufficiently small values of the viscous parameters for the instability to be only partially regulated, we find it more convenient to directly specify $\eta_{0}=5\times 10^{-4}$.

We plot in Fig. 8 the number density initial data and a snapshot of the inviscid solution with $A=0$ (i.e., zero initial polar velocity) at a later time when the Kelvin-Helmholtz rolls have formed. At this time, the characteristic rolls form in an ace-of-spades-like pattern in both the number and vorticity densities, while troughs in the energy density form at the locations of the eight rolls. As the number of winds grows in each roll, the jet thins and pinches off in between each of the four ace-of-spades shapes. The corresponding evolution with $A=0$ by the BDNK equations almost completely shears the rolls—sheared rolls are just about visible in the right panel of Fig. 9, wherein we plot the number density in region I of the cubed-sphere grid at $t=50$.

For a fixed value of $\eta_{0}$, the viscosity regulates the instability to a lesser degree when the initial polar velocity is nonzero. We show in Fig. 10 snapshots of the inviscid and viscous number and vorticity density solutions with $A=0.01$. The inclusion of a non-zero initial polar velocity causes the jet to deform in shape and allows the rolls in the viscous solution to reach a further stage of their growth than in the $A=0$ case, with the shearing of the inviscid rolls apparent in the upper row. The lower row of Fig. 10 also demonstrates how viscosity smears out both the rolls in the vorticity density and the thin lines which connect them in the inviscid solution (also note the extrema of the vorticity density colorbar for $\eta_{0}=5\times 10^{-4}$ is roughly half that for $\eta_{0}=0$). That is, we see that viscosity has the expected effect of giving rise to the phenomenon of vorticity diffusion.

We plot in Fig. 11 the integral of the vorticity density over $S^{2}$ (for the $A=0.01$ simulations), which is a conserved quantity under evolution by the Euler equations. There, we see that vorticity is conserved at early times in our inviscid numerical simulation at different resolutions until $t\approx 30$, by which time the turbulent cascade to increasingly short length scales is at an advanced stage (see Fig. 10) and the convergence factor in the numerical simulations is rapidly falling (see Fig. 20). In contrast, the integrated vorticity density is more stable throughout the duration of the BDNK evolution, in which viscosity regulates the turbulent cascade present in the inviscid flow.

V Discussion

We have presented the first numerical solutions to the BDNK equations for a ($2+1$)D flow on a spherical manifold, specializing to the case of a $4$D conformal fluid constrained to the surface of a geometric sphere in Minkowski spacetime. Our numerical scheme consists of solving the equations of motion discretized in cubed-sphere coordinates using a fourth-order-accurate method of lines with a centered finite-difference spatial discretization and an explicit Runge-Kutta time integration. We applied this numerical implementation to three distinct sets of initial data designed to test the correctness of our code and to study out-of-equilibrium dynamics governed by the BDNK equations.

We compared the oscillation frequencies and damping timescales present in our numerical simulations of perturbed equilibrium states to analytic predictions obtained from the linearized equations of motion. We found good agreement between the results of the linear analysis and the numerical simulations. In another test case, we compared the inviscid and viscous evolution of Kelvin-Helmholtz-unstable initial data, qualitatively demonstrating the viscous shearing of the characteristic rolls that form in the inviscid solution and the viscous diffusion of vorticity.

We also studied the ensuing dynamics from smooth, stationary Gaussian initial data initialized in local thermodynamic equilibrium. Viscosity in the BDNK solution with $\eta/s=1/(4\pi)$ damps the high-frequency components present in the inviscid evolution by the Euler equations, smoothing the shock that forms in the inviscid flow and allowing the BDNK numerical simulation to converge to a late-time equilibrium state. For larger values of $\eta/s$, the relative power of the high-frequency components in the solution grows as increasingly steep gradients form in the solution. In the $\eta/s\in\{10,20\}/(4\pi)$ BDNK solutions, sufficiently steep gradients form such that convergence is lost in finite time at all grid resolutions we have considered. As the time at which convergence is lost is approached in the numerical simulations, the BDNK regime-of-validity diagnostics rapidly assume extreme values, indicating divergence away from local equilibrium and the regime of validity of first-order hydrodynamics. In particular, when convergence is lost in the $\eta/s=20/(4\pi)$ solution, the first-order terms in the $T^{tt}$ component dominate the perfect-fluid terms with $|T^{tt}_{(1)}|\sim 20|T^{tt}_{(0)}|$, and the weak-energy condition is severely violated with $u_{\mu}u_{\nu}T^{\mu\nu}\sim-3$ at the location of steepest gradients in the solution.

We observe a similar onset of steep gradients and divergence away from local equilibrium in numerical solutions to the BDNK equations for a planar-symmetric $4$D conformal fluid in Minkowski spacetime, as presented in App. A. These effects are also present at all resolutions considered in the ($1+1)$D simulations, for which we are able to achieve much higher grid resolutions due to the reduced dimensionality of the flow, providing strong evidence that discontinuities form in the corresponding continuum solutions. This appears to also be the case for the ($2+1$)D flows, though verifying this by achieving similarly high grid resolutions would be too computationally expensive. These solutions constitute numerical evidence that solutions to the BDNK equations can develop singularities in finite time from smooth initial data.

The development of increasingly steep gradients with viscosity is surprising behavior to observe in a dissipative theory. In principle, the fact that singularities can develop from smooth initial data in a viscous relativistic theory of hydrodynamics is not new—this has been established for MIS-type theories with only bulk viscosity present [47]. However, the out-of-equilibrium behavior of the BDNK and MIS equations are very different, so this behavior being possible in the latter does not imply it is possible in the former. Moreover, one intuitively expects increasing the shear viscosity to make the diffusion of momentum and thus the smoothing of gradients in the fluid more efficient. Though we observe this to indeed be the case for small viscosities [$\eta/s\sim 1/(4\pi)$], it is rather puzzling that large viscosities [$\eta/s\in\{10,20\}/(4\pi)$] are associated with the onset of increasingly steep gradients in the solution. Though it is not obvious what causes this behavior at large viscosities, it is likely an inherently non-linear effect—we explicitly demonstrate linear mode stability of equilibrium states in App. C, verifying that no instabilities arise at the linear level (as one expects from the sufficient frame conditions that should guarantee this in BDNK theory). The presence of this behavior in both the ($2+1$)D flows on the two-sphere and ($1+1$)D planar-symmetric flows suggest this behavior is not unique to the geometry of the two-sphere, which, in the case of the inviscid flow, has the effect of focusing the energy density at the poles until a shock forms in the solution.

If singularities do indeed develop in the underlying continuum BDNK solutions at large values of $\eta/s$, this implies one may need to place additional restrictions on transport coefficients, beyond those for well-posedness, to be able to use BDNK to model physical systems. As discussed at the end of Sec. II.2, values of $\eta/s$ from measurements of several strongly-interacting systems (such as finite nuclear matter and the QGP) nearly saturate the KSS bound $\eta/s\sim 1/(4\pi)$, and values of $\eta/s\sim\{10,20\}/(4\pi)$, where we start observing problems with the BDNK solutions, may be unphysically high222Note that we are considering relativistic scenarios with very high initial energy density/temperature, such that fluid velocities tens of percent the speed of light develop within the timescale of our simulations. Therefore that, as mentioned in the introduction, water at standard temperature and pressure (STP) has a much higher viscosity-to-entropy ratio relative to the KSS bound than scenarios studied here does not imply BDNK cannot be used to model water at STP. And, of course, water is not well approximated by a conformal fluid.. On the other hand, hydrodynamic frame independence of the stress tensor and baryon current only hold to linear order for perturbations about ideal fluid solutions, and the singularities we observe, which occur in regimes with large non-ideal contributions to the stress tensor, may again be due to a choice of a “bad” hydrodynamic frame. Thus, an important potential avenue for future work would be to investigate whether or not the observed large $\eta/s$ behavior is in any way generic for different hydrodynamical frames, as well as other classes of initial data. If singularities can form in BDNK hydrodynamics, it may also be interesting to investigate if they could be made sense of as weak-form solutions, in analogy with singularities forming in the Euler equations.

There are several additional avenues for future work. Generalizing our code to fluids without conformal symmetry would be important for considering astrophysically-motivated equations of state and studying the effects of the full range of first-order transport coefficients, such as bulk viscosity and heat conductivity. Coupling BDNK to curved spacetime in spherical symmetry would be a first step towards using BDNK theory to, for example, study viscous effects in gravitational collapse. A similar implementation of the BDNK equations as in this work, but for an intrinsically $3$D conformal fluid on the two-sphere, could be used to study how viscosity affects the late-time behavior of the inverse turbulent cascade observed in numerical solutions to the Euler equations studied in Ref. [29]. Lastly, a tensor product of cubed-spheres in the radial direction could extend our numerical scheme to full ($3+1$)D simulations of the BDNK equations, which would be applicable to a wide variety of astrophysical systems.

Acknowledgements.

We thank Marcelo Disconzi for useful discussions regarding this work. The simulations presented in this article were performed on computational resources managed and supported by Princeton Research Computing, a consortium of groups that includes the Princeton Institute for Computational Science and Engineering (PICSciE), and the Office of Information Technology’s High Performance Computing Center and Visualization Laboratory at Princeton University. FP acknowledges support from the NSF through the grants PHY-220728 and PHY-2512075.

Appendix A Numerical Simulations in Planar Symmetry

In this appendix, we present numerical solutions to the BDNK equations for a $4$D conformal fluid in Minkowski spacetime with variations in only a single spatial direction. The BDNK equations were numerically solved in this setting in Ref. [46], wherein the out-of-equilibrium dynamics of smooth, stationary Gaussian initial data was studied alongside viscous shock wave solutions. In that work, the numerical simulations with Gaussian initial data demonstrated the behavior of a universal hydrodynamic attractor; for increasing values of $\eta/s$, the solutions to the BDNK and a truncated version of MIS-type equations disagreed increasingly at early times as each underwent different transient dynamics before eventually settling down to similar late-time equilibrium states.

Here, we build on this study of out-of-equilibrium dynamics as governed by the BDNK equations with the same smooth, stationary Gaussian initial data. Our numerical solutions converge to a late-time equilibrium state for $\eta/s\in\{1,3,10\}/(4\pi)$, but, for $\eta/s=20/(4\pi)$, the energy density and velocity profiles form step-like discontinuities at all grid resolutions we have considered. As the resolution is increased, we converge to a finite value of time at which convergence is lost, suggesting the underlying continuum solution develops a discontinuity at this time, as opposed numerical instabilities. These simulations, in combination with the similar qualitative behavior in the ($2+1$)D flows presented in Sec. IV.2, provide evidence that solutions to BDNK equations can develop discontinuities from smooth initial data. Further, if continuum discontinuities do indeed develop, these simulations also suggest that for sufficiently large values of $\eta/s$, the hydrodynamic-attractor-like behavior observed in Ref. [46] can be violated, with the solution instead diverging away from local equilibrium before convergence is lost due to the formation of these continuum discontinuities, preventing an evolution of the flow to a late-time equilibrium state.

We have compared solutions obtained using our finite difference numerical scheme to those from a finite volume code333This code is publicly available on Github. written by the authors of Ref. [52]. For $\eta/s=20/(4\pi)$ we find good agreement between the codes up until the time at which discontinuities appear to form in the solution, after which neither code crashes but the solution profiles differ slightly near the propagating step-like discontinuities (with oscillations breaking out in $\dot{\epsilon}$ in the finite volume code). If a continuum discontinuity does indeed form, then one should not expect the two codes to agree at later times. Our finite difference code assumes a smooth continuum solution and is thus unable to resolve any continuum discontinuities. The finite volume code is based on high-resolution shock-capturing (HRSC) methods developed for handling shocks which can form in the Euler equations, in which case one is able to obtain physical weak-form solutions by requiring consistency with the conservation of stress-energy and the second law of thermodynamics. However, it is not clear whether the analogous Riemann problem for BDNK can be made sense of, in part because the underlying PDEs are second order. While it appears sensible to apply Euler-derived HRSC methods to the BDNK equations in the case of unresolved smooth shocks (i.e., when the length scale on which a shock is smoothed by viscosity is much smaller than the grid resolution—see Refs. [46, 52] for a further discussion), it is unclear whether applying these same methods to a continuum shock makes any physical sense, and thus whether the subsequent evolution is physical or not.

In the remainder of this section, we describe in further detail our numerical simulations summarized above. We first briefly describe the equations of motion, our initial data and numerical scheme in App. A.1. We then present our results in App. A.2. To illustrate the observed discontinuous behavior is not an unresolved numerical artifact, we also include in App. A.2 a discussion about convergence tests and estimates of the time when the continuum solution appears to develop a discontinuity.

A.1 Equations, Initial Data, and Numerical Scheme

The equations of motion governing a $4$D planar-symmetric conformal flow in Minkowski spacetime were derived in Ref. [46]. We thus only provide a brief review and refer the reader to Ref. [46] for further details.

We describe the flow using Cartesian coordinates, with respect to which the line element $\mathop{}\!\mathrm{d}{s}^{2}=\eta_{\mu\nu}\mathop{}\!\mathrm{d}{x}^{\mu}\mathop{}\!\mathrm{d}{x}^{\nu}$ reads

[TABLE]

We parametrize the four-velocity of the flow as

[TABLE]

where the Lorentz factor $W=\left(1-v^{2}\right)^{-1/2}$. Since we consider a conformal fluid with zero baryon number density, the flow is completely characterized by the energy density $\epsilon=\epsilon(t,x)$ and velocity $v=v(t,x)$. The coordinate components of the BDNK tensor (6a) are functions of only $\epsilon$, $v$ and their spatial and temporal derivatives, so that the equations of motion reduce to a coupled system of two second-order nonlinear PDEs in $\epsilon$ and $v$. We relegate the tensor component expressions and PDEs to App. B.1, which can also be found in Ref. [46].

We consider an initially stationary Gaussian profile in the energy density:

[TABLE]

with $\epsilon_{0}=0.1$, $A=0.4$, and $w=5$ for values of $\eta/s\in\{0,1,3,10,20\}/(4\pi)$. Since the BDNK system consists of two PDEs which are second order in time (and space), we must also prescribe initial data for the time derivatives of $\epsilon$ and $v$. Instead of explicitly specifying $\dot{\epsilon}$ and $\dot{v}$ at time $t=0$, we instead choose to initialize $T_{(1)}^{tt}=0=T_{(1)}^{tx}$ which determines $\dot{\epsilon}$ and $\dot{v}$ (see Eqs. A1-A2 in Ref. [46]).

We consider a spatial domain $x\in\left[-L,L\right]$ with $L=100$ and periodic boundary conditions. That is, if we label the grid points as $i=1,2,\ldots,N$, the value of the discretized hydrodynamic variables at grid location $i=N+1$ is given by that at $i=1$. We numerically solve the equations of motion via the method of lines with fourth-order-accurate finite-difference stencils in space, an explicit fourth-order-accurate Runge-Kutta time integration, and Kreiss-Oliger dissipation to damp high-frequency components in the solution which otherwise cause the numerical simulation to crash at early times.

Since our numerical scheme is fourth-order accurate in time and space, we expect a convergence factor (97) $Q_{N}(t)=2^{4}=16$. For both the solutions in planar symmetry and on the two-sphere, we find that the convergence factor is sometimes twice the value expected from the order of the numerical scheme, likely a result of symmetries in the initial data or equations of motion. This behavior is observed in the convergence factor for all the BDNK solutions presented in this section, shown in Figs. (13, 17). We show the $\eta/s\in\{0,20\}/(4\pi)$ convergence factors in the following section to illustrate how the behavior of $Q_{N}(t)$ in the solution to the Euler equations, which itself develops a shock (as also observed in Ref. [46]), is similar to that for the BDNK solution $\eta/s=20/(4\pi)$. More technical details about the computation of $Q_{N}(t)$ and convergence tests can be found in App. E.

Lastly, we note that we observe similar qualitative behavior in numerical solutions obtained using a second-order-accurate implicit Crank-Nicholson discretization/evolution scheme.

A.2 Results

In Fig. 12, we plot the energy density for the $\eta/s\in\{0,1,3\}/(4\pi)$ solutions at time $t=38$. The Gaussian profile, which is initially peaked at $x=0$, separates into two clumps which travel in opposite directions towards each boundary of the spatial domain. In the inviscid [$\eta/s=0/(4\pi)$] solution, there is a sharp steepening of the profile at locations $x\approx\pm 32$. At these locations, the upstream velocity at each front is larger than the downstream velocity, causing the energy and velocity profiles to steepen and form a shock in the absence of dissipation. Dissipation smooths these gradients in the $\eta/s\in\{1,3\}/(4\pi)$ solutions. In both of these viscous cases, the clumps are further damped as they propagate to either end of the spatial boundary, whereupon they re-enter on the opposite side (due to the periodic boundary conditions) and collide with the other clump, after which this process restarts and continues until a steady-state equilibrium solution is reached at late times with $v=0$ and $\epsilon=\mathrm{constant}$.

Since the numerical schemes employed throughout this work are based on finite differences, which assume a smooth underlying continuum solution, we are unable to resolve discontinuities such as those which form in the inviscid solution. If a discontinuity forms in the continuum solution at time $t=t_{d}$, we expect convergence to be lost at time $t_{c}<t_{d}$, and, with increasing resolution, we expect $t_{c}$ to approach $t_{d}$ from below, but never quite reach $t_{d}$ since the discontinuity is unresolvable with our chosen numerical scheme. We observe this behavior in convergence tests of the solution to the Euler equations. We plot in the left panel of Fig. 13 the convergence factor $Q_{N}(t)$ of the inviscid solution, which, with increasing resolution, sharply drops to unity (indicating a loss of convergence) closer and closer to some finite time $t\sim 40$.

The numerical BDNK solutions with $\eta/s\in\{1,3,10\}/(4\pi)$ exhibit qualitatively similar behavior: two clumps form which are progressively damped until an equilibrium state is reached at late times. The dynamics of the $\eta/s=20/(4\pi)$ solution markedly differ from the less viscous flows. As illustrated in Fig. 14, four clumps form in the energy density profile instead of two, and step-like discontinuities form in the energy density (and velocity) profile, at the location of which $\dot{\epsilon}$ blows up. As these steep gradients form in the solution, the monitored diagnostics suggest the flow exits the regime of validity of first-order hydrodynamics. We plot in Fig. 15 the maximum and minimum value of the diagnostic quantities $|T_{(1)}^{tt}/T_{(0)}^{tt}|$ and $u_{\mu}u_{\nu}T^{\mu\nu}$, respectively, over the spatial domain as a function of time until convergence is lost. Therein, one can see the weak-energy condition is never violated and $|T^{tt}_{(1)}/T^{tt}_{(0)}|$ never exceeds unity for $\eta/s=10/(4\pi)$, while, for $\eta/s=20/(4\pi)$, the diagnostics rapidly assume extreme values just before convergence is lost in the simulation, at which time $|T^{tt}_{(1)}|\sim 10|T^{tt}_{(0)}|$ and $u_{\mu}u_{\nu}T^{\mu\nu}\sim-1$.

The convergence factor $Q_{N}(t)$ for solutions with $\eta/s=20/(4\pi)$ at four different grid resolutions is shown in the right panel of Fig. 13. The behavior of the convergence factor is qualitatively similar to that of the inviscid solution: with increasing resolution, convergence in the $\eta/s=20/(4\pi)$ solution is maintained incrementally closer to some finite time $t_{d}$, before which $Q_{N}(t)$ rapidly drops to unity, indicating a loss of convergence. In our highest resolution run, which consists of $N=2^{14}+1$ points in the spatial domain, the independent residual is on the order of machine precision at early times (which is why the red curve in Fig. 13 is much lower than the nominal convergence factor at early times). That the behavior of the convergence factor near the time convergence is lost is persistent at such high grid resolutions suggests the development of these discontinuities is not an unresolved numerical artifact.

If we are indeed attempting to resolve a discontinuity which forms at some time $t_{d}$ in the underlying continuum solution, then, at this time, we expect the spatial derivatives to be infinite at the locations of the step-like transitions in the continuum solution. We now assume this behavior in our numerical solutions to obtain an estimate for $t_{d}$ and qualitatively show that, with increasing resolution, we converge to some finite $t_{d}$, providing further evidence that such a time $t_{d}$ exists.

As a measure of the assumed derivative blow up, we plot in (the background of) Fig. 16 the maximum absolute value of $\epsilon^{\prime}=\partial_{x}\epsilon$ over the spatial domain as a function of time. Since we are interested in the growth of the derivatives at times close to the onset of steep gradients while convergence is maintained, we mask the time series of the three highest resolutions to times $t\in[t_{1},t_{2}]$, where $t_{1}$ is the time when convergence is lost (which we define as when $Q_{N}(t)<2$) at the lowest resolution ($N=2^{11}+1$) and $t_{2}$ is the time at which the higher resolutions lose convergence respectively. These masked time series correspond to the curves in the foreground of Fig. 16. We fit each of these curves to the functional form

[TABLE]

to estimate $t_{d}$.

The behavior of the times $t_{c}$ at which convergence is lost and the fitted values of $t_{d}$, which are listed in Tab. 3, exhibit a sandwich-like behavior with increasing grid resolution. At the resolution is increased, $t_{c}$ is incremented by an increasingly small amount, while the fitted value of $t_{d}$ decreases by an increasingly small amount such that the two values appear to converge closer and closer to some middling value, presumably the true value of $t_{d}$.

Appendix B Tensor Coordinate Components and Equations of Motion

In this appendix we present the coordinate components of the BDNK tensor (6a) and the resulting equations of motion for a $4$D conformal fluid in Minkowski spacetime for (i) a planar-symmetric flow described in Cartesian coordinates (App. B.1) and (ii) a flow constrained to the surface of a geometric sphere of radius $R$ (App. B.2). The equations in planar symmetry were originally derived in Ref. [46].

B.1 Planar Symmetry

We consider a conformal fluid with four-velocity

[TABLE]

where the Lorentz factor $W=(1-v^{2})^{-1/2}$. The BDNK tensor for such a fluid has components

[TABLE]

where each component consists of a zeroth and first-order contribution, e.g., $T^{tt}=T_{(0)}^{tt}+T_{(1)}^{tt}$. The two non-trivial equations of motion are

[TABLE]

The zeroth-order contributions to the stress-energy tensor components are given by

[TABLE]

The first-order terms $T_{(1)}^{\mu\nu}$ are constructed from the out-of-equilibrium quantities $\mathcal{A}$, $\mathcal{Q}^{\mu}$ and $\mathcal{T}^{\mu\nu}$, which take on the form

[TABLE]

Orthogonality to $u^{\mu}$ implies $\mathcal{Q}^{t}=v\mathcal{Q}^{x}$ and planar symmetry gives $\mathcal{Q}^{y}=0=\mathcal{Q}^{z}$. Planar symmetry also requires the transverse-traceless tensor $\mathcal{T}^{\mu\nu}$ to assume the form

[TABLE]

with $\mathcal{T}^{yy}=-\partial_{\mu}u^{\mu}=\mathcal{T}^{zz}$. Orthogonality to $u^{\mu}$ requires $\mathcal{T}^{tt}=v\mathcal{T}^{tx}$, $\mathcal{T}^{tx}=v\mathcal{T}^{xx}$.

All the zeroth and first-order terms in the stress-energy tensor components are determined by (50, 51). Substitution of the tensor components into (49) then gives the equations of motion explicitly in terms of the two variables $\epsilon$ and $v$.

B.2 Two-Sphere

The four-velocity of a fluid constrained to a geometric sphere of radius $R$ is given by

[TABLE]

where each component is a function of $t$, $\theta$ and $\phi$. The BDNK tensor for a fluid with four-velocity (53) in the background of the metric (17) assumes the form

[TABLE]

where we again separate each component into terms zeroth and first order in gradients. With this notation, the three non-trivial equations of motion are

[TABLE]

The zeroth-order tensor components are given by

[TABLE]

The out-of-equilibrium correction to the energy density takes on the form

[TABLE]

The components of $\mathcal{Q}^{\mu}$ are given by

[TABLE]

The symmetric, transverse-traceless tensor $\mathcal{T}^{\mu\nu}$ has five independent components:

[TABLE]

The traceless condition $\mathcal{T}^{\mu}{}_{\mu}=0$ implies

[TABLE]

while the transverse condition $u_{\mu}\mathcal{T}^{\mu\nu}=0$ gives

[TABLE]

which determines another four components. The six remaining components are given by $\mathcal{T}^{\mu\nu}=\mathcal{T}^{\nu\mu}$.

The out-of-equilibrium tensor components $T^{\mu\nu}_{(1)}$ are determined by $\mathcal{A}$, $\mathcal{Q}^{\mu}$ and $\mathcal{T}^{\mu\nu}$ as given by (57-61) while the zeroth-order components are given by (56). Substitution of these components into (55) gives the equations of motion explicitly in terms of the three variables $\epsilon$, $u^{\theta}$, and $u^{\phi}$.

Appendix C Linear Normal Mode Analysis

In this appendix, we consider fixed-background, linear spherical harmonic perturbations of stationary, uniform equilibrium states of the Euler and BDNK equations. We explicitly demonstrate linear stability of this class of BDNK equilibrium states for hydrodynamic frames (13) with

[TABLE]

where $\eta_{0}>0$ and $a>b>0$. We perturb the four-velocity of the fluid using even and odd-parity vector spherical harmonics. We demonstrate mode stability in the high-frequency $l\to\infty$ limit for even perturbations and for all $l$ for odd perturbations.

In Sec. C.1, we describe the form of the perturbations we consider and derive the covariant BDNK stress-energy tensor to first order in the perturbations. We then consider separately the cases of even and odd perturbations in Secs. C.2-C.3, respectively.

C.1 Perturbed Stress-Energy Tensor

C.1.1 Hydrodynamic Perturbations

We consider perturbations of the form

[TABLE]

where $\epsilon_{0}=\rm{constant}$, $u^{\mu}_{0}=[1,0,0,0]^{\mathrm{T}}$, and $\delta{\epsilon}$ and $\delta u^{\mu}$ are expressed, respectively, in terms of scalar and vector spherical harmonics (see Eq. 28 for a definition of the latter):

[TABLE]

where $l\geq 0$, $|m|\leq l$, $j\in\{2,3\}$, $\delta_{\epsilon}/\epsilon_{0}\ll 1$ and $\delta_{u}\ll 1$. The even and odd-parity four-velocity perturbations are given explicitly by

[TABLE]

From (65), it follows that $\delta{u}=0$ for both parities if $l=m=0$.

C.1.2 Stress-Energy Components

The conformal BDNK tensor is given by [48]

[TABLE]

where

[TABLE]

and the transverse traceless symmetric tensor $\sigma^{\mu\nu}$ is defined in (8b). To derive the perturbed $T^{\mu\nu}$ components, we first consider the scalar and tensor quantities (67) perturbed to linear order.

Since we consider a background four-velocity $u_{0}^{\mu}=\delta^{\mu}{}_{t}$, we have

[TABLE]

where $\Gamma^{\nu}_{\mu\alpha}$ denote the Christoffel symbols associated with the Minkowski metric, which have $\Gamma^{\nu}_{\mu t}=0$ since the metric is stationary and has constant $g_{t\mu}$ components. Under the perturbations (63, 64), we thus have $\nabla_{\mu}u^{\nu}\to\nabla_{\mu}(\delta{u})^{\mu}$, while $\epsilon_{0}=\mathrm{constant}$ gives $\nabla_{u}\epsilon\to\nabla_{\mu}(\delta{\epsilon})$. The quantities $\mathcal{E}$, $\mathcal{Q}^{\mu}$ and $\pi^{\mu\nu}$ transform as $\mathcal{E}\to\epsilon_{0}+\delta\mathcal{A}$, $\mathcal{Q}^{\mu}\to\delta\mathcal{Q}^{\mu}$ and $\pi^{\mu\nu}\to\delta\pi^{\mu\nu}$, where, to first order,

[TABLE]

and we have defined the zeroth-order projection tensors

[TABLE]

The conformal BDNK tensor (67) then transforms as

[TABLE]

Parametrizing the transport coefficients $\left\{\eta,\lambda,\chi\right\}$ as in (13) means that, for example,

[TABLE]

Since the class of equilibrium states we consider are constant over the sphere, the transport coefficients in (69) always multiply $\mathcal{O}(\delta)$ terms, so they need only be kept to $\mathcal{O}(1)$. That is, we take $\eta\to\eta_{0}\epsilon_{0}^{3/4}$, $\lambda\to\lambda_{0}\epsilon_{0}^{3/4}$ and $\chi\to\chi_{0}\epsilon_{0}^{3/4}$ in (69).

Explicit coordinate components of the perturbed stress-energy tensor (71) can be obtained by substituting the perturbations (64, 65) into (69). The linearized equations of motion, $\nabla_{\mu}T^{\mu\nu}=0$, then provide constraints (dispersion relations) on the permissible values of the perturbation parameters $\{\delta_{\epsilon},\delta_{u},\omega\}$. In the next two sections, we analyze these constraints for even and odd perturbations of the Euler and BDNK equations for general values of $l$ and $m$.

C.2 Even Perturbations

C.2.1 Euler Fluid

The linearized even-perturbed Euler equations provide the following two constraints

[TABLE]

which have non-trivial solutions

[TABLE]

In the $l\to\infty$ limit, $\omega\sim l/(R\sqrt{3})$.

C.2.2 BDNK Fluid

The linearized equations of motion governing even perturbations of a conformal BDNK fluid reduce to the two constraints

[TABLE]

Setting $\delta_{\epsilon}=0$ trivially requires $\delta_{u}=0$. Assuming $\delta_{\epsilon}\neq 0$, if $l=m=0$, we have $\omega=0$ or $\omega=-4i\epsilon_{0}^{1/4}/(3\chi_{0})$ (the latter being a frame mode), otherwise if $l\geq l$ we obtain from (75a)

[TABLE]

The right-hand side of (76) is well-defined, for if the denominator is equal to zero, we can substitute the resulting solution for $\omega$ back into (75a) to obtain

[TABLE]

which cannot be satisfied since each term on the right-hand side is strictly positive. The substitution of (76) into (75b) yields the following equation for $\omega$

[TABLE]

where

[TABLE]

In the $l\to\infty$ limit, (78) becomes

[TABLE]

where

[TABLE]

We construct a solution to (80) by expanding $\omega$ in powers of $l$:

[TABLE]

We substitute the expansion (82) into (80) and solve the resulting equation order by order in powers of $l$ to determine the coefficients $\omega_{n}$. For $n>1$, the $\omega^{4}$ term in (80) strictly dominates the remaining terms as $l\to\infty$, so that one trivially recovers $\omega_{n}=0$. For $n=1$, the terms which scale as $l^{4}$ yield the asymptotic relation

[TABLE]

which has four real solutions

[TABLE]

The next-to-leading-order ($n=0$) behavior of $\omega(l)$ is determined by the equation

[TABLE]

which gives

[TABLE]

For a choice of frame $(\lambda_{0},\chi_{0})=(a\eta_{0},b\eta_{0})$ with $a>b>0$, it follows from (86) that $\mathrm{Im}(\omega_{0})<0$, thus proving linear mode stability in the $l\to\infty$ limit.

C.3 Odd Perturbations

C.3.1 Euler Fluid

The linearized equations of motion governing odd perturbations of the Euler equations reduce to the constraints

[TABLE]

which, for a non-trivial time dependence $\omega\neq 0$, only have a trivial solution $\delta_{\epsilon}=\delta_{u}=0$.

C.3.2 BDNK Fluid

The linearized BDNK equations provide the following constraints on the odd perturbations

[TABLE]

If $l=m=0$, we have $\omega=0$ or $\omega=-4i\epsilon_{0}^{1/4}/(3\chi_{0})$ (the latter being a frame mode), which is the same solution as in the even-parity $l=m=0$ case since the four-velocity perturbations vanish, so these perturbations are equivalent. For $l\geq 1$, the constraints (88) are satisfied for arbitrary $\delta_{u}$, $\delta_{\epsilon}=0$ and

[TABLE]

which has $\mathrm{Im}(\omega)<0$ provided $\eta_{0},\lambda_{0}\geq 0$ (since $l\geq 1$).

Appendix D Cubed-Sphere Coordinates

In this appendix, we describe in further detail the multi-block cubed-sphere grid we employ in our numerical scheme on the two-sphere.

We provide a geometric illustration of the cubed-sphere grid in Fig. 1. In the left panel therein, we depict the solid angle projection of one face of a cube onto a sphere (with the cube placed concentrically with respect to the sphere). This defines some region on the sphere, and the collection of the six projected regions corresponding to each face of the cube covers the surface of the sphere with no gaps or overlap between the regions. The inverse transformation allows one to project each of these six regions on the sphere onto its own face of the cube, giving rise to the grid structure shown in the right panel of Fig. 1. In this grid, region I is identified with $x=\mathrm{constant}>0$ face of the cube, region II with the $y>0$ face, region III with the $x<0$ face, region IV with the $y<0$ face, and regions V and VI with the $z>0$ and $z<0$ faces, respectively (so that the north and south poles of the sphere lie in regions V and VI, respectively). Each of these grids is regular, free from coordinate singularities and can be endowed with a simple uniform structure allowing for the use of simple finite-difference methods. In particular, we choose to make each patch on the cubed-sphere a uniform $N\times N$ grid, so that all the regions are characterized by the same discretization scale $h$. References to the number $N$ of grid points on the two-sphere in App. E corresponds to a total of $6N^{2}$ points across all six patches.

We denote the local coordinates on each grid by $X$ and $Y$ (see the right panel of Fig. 1), where the coordinate transformations for each region are given by

• •

Region I ($x=\mathrm{constant}>0$)

[TABLE]

• •

Region II ($y=\mathrm{constant}>0$)

[TABLE]

• •

Region III ($x=\mathrm{constant}<0$)

[TABLE]

• •

Region IV ($y=\mathrm{constant}<0$)

[TABLE]

• •

Region V ($z=\mathrm{constant}>0$)

[TABLE]

• •

Region VI ($z=\mathrm{constant}<0$)

[TABLE]

where $r=\sqrt{x^{2}+y^{2}+z^{2}}$, $-1\leq X\leq 1$, and $-1\leq Y\leq 1$ for all regions. By computing the Jacobian of the transformations (90-95) and performing a change of coordinates from either Cartesian or spherical coordinates to cubed-sphere coordinates, one can show that the Minkowski metric takes on the same form for all regions [63]:

[TABLE]

With the cubed-sphere metric (96) and coordinates (90-95), the Euler and BDNK equations can be written in terms of cubed-sphere coordinates. We discretize these equations in space (i.e., in the local coordinates $X$ and $Y$) using centered finite difference stencils (see the discussion in Sec. III.2 regarding our treatments of boundary points for a given patch) and integrate the discretized equations in time using a fourth-order Runge-Kutta scheme. We take the patches to be uniform and thus characterized by a single discretization scale $h$.

Appendix E Convergence Tests

In this appendix, we present convergence tests for all the ($1+1$)D and ($2+1$)D simulations presented in this work except the $\eta/s\in\{0,20\}/(4\pi)$ ($1+1$)D simulations, which are presented in Fig. 13.

In all our simulations, we employ a fourth-order-accurate numerical scheme which is characterized by a single discretization scale $h$ (or, equivalently, number of grid points $N$). The discretized solution obtained from such a scheme only satisfies the continuum equations of motion to within an $\mathcal{O}(h^{4})$ truncation error, which vanishes in the continuum limit $h\to 0$. If the discretized solution obtained from our numerical scheme approaches a solution of the continuum PDEs in the limit $h\to 0$, then the discretized solution should also give an $\mathcal{O}(h^{4})$ truncation error when plugged into a discretization of the equations of motion performed using a distinct numerical scheme—we denote this “independent residual” by $R_{N}$.

We assess the validity of our discretized solutions by monitoring the convergence factor

[TABLE]

where $||\cdot||_{1}$ denotes the vector one-norm and $R_{N}$ is computed at each point in the spatial domain by plugging the discretized solution obtained from the method of lines into the equations of motion discretized using a leapfrog scheme. For a fourth-order-accurate scheme and a sufficiently smooth continuum solution, we expect $Q_{N}(t)\sim 16$. As discussed in App. A, we sometimes observe a convergence factor greater than (and typically twice) the value expected value for a fourth-order scheme, likely due to underlying symmetries present in the initial data or equations of motion which cause the lowest-order error functions in the Richardson expansion of the discretized solution to vanish.

We plot in Fig. 17 the convergence factor for the $\eta/s\in\{1,3,10\}/(4\pi)$ ($1+1$)D simulations. The remaining figures show the convergence factor for the ($2+1$)D simulations: Fig. 18 corresponds to the numerical simulations of perturbed fluid states discussed in Sec. IV.1.3, Fig. 19 corresponds to the Gaussian initial data simulations presented in Sec. IV.2, and Fig. 20 corresponds to the evolution of Kelvin-Helmholtz-unstable initial data discussed in Sec. IV.3.