Two numerical routes for simulating relativistic viscous fluids

Viscosity complicates the description of fluids moving at relativistic speeds. Some older relativistic viscous theories allow disturbances to propagate faster than light or develop unphysical instabilities. BDNK hydrodynamics, named for Bemfica, Disconzi, Noronha and Kovtun, instead permits causal propagation and linear stability when its parameters satisfy appropriate conditions.

In a preprint submitted to arXiv on November 6, 2026, Teerthal Patel and colleagues develop two numerical formulations of this theory and establish mathematical results supporting their use. This account is based on the supplied complete extracted text of arXiv version 2611.09582v1 and its original metadata; it reports the manuscript’s proofs and simulations, not an independent verification.

Two ways to evolve the same fluid

BDNK is called a first-order hydrodynamic theory because its constitutive relations include first derivatives of fluid quantities. Nevertheless, conservation of energy and momentum produces second-order differential equations. Turning those equations into a practical simulation requires care: reformulating them must preserve their physical and mathematical properties.

The researchers take two routes. Their fully first-order formulation promotes spatial derivatives to independently evolved variables. These variables must remain consistent with the underlying fluid fields. Additional “cleaning” fields help control numerical violations of related constraints. The general system, including the particle-current sector and cleaning fields, evolves 44 variables.

The mixed-order formulation reduces the time derivatives to first order while retaining second-order spatial derivatives. It evolves nine variables and needs no independently evolved spatial-gradient fields or associated cleaning fields. Spatial derivatives enter its fluxes and sources, requiring different numerical treatment.

The authors prove that both formulations recover the original BDNK equations under the relevant conditions. For the first-order system, this requires compatible initial gradient data and zero initial cleaning fields. Cleaning does not, by itself, guarantee that every possible constraint error decays.

Mathematical guarantees have boundaries

The paper establishes local existence, uniqueness and continuous dependence on initial data in specified mathematical settings. “Local” means that the results guarantee evolution for some time interval, not indefinitely. Continuous dependence is established in Sobolev topologies below the highest controlled regularity.

These results require characteristic-speed conditions, sufficient initial-data regularity and analytic dependence of the relevant coefficients. The principal fluid theorems use prescribed analytic spacetime backgrounds and periodic spatial domains. They do not establish the corresponding result for a fluid evolving together with Einstein’s gravitational equations.

A broader mathematical contribution gives sufficient structural conditions for solutions of certain mixed-order systems to satisfy strictly hyperbolic equations. This removes repeated characteristic factors while retaining distinct characteristic roots, providing a route to the existence and uniqueness results.

The authors also simplify the selection of admissible hydrodynamic frames: prescriptions defining temperature and velocity away from equilibrium. Their practical parameterization covers a broad class of equations of state, but assumes vanishing thermal conductivity. These frames are definitions of fluid variables, not simply different observer viewpoints.

Tests from expanding plasmas to black-hole flows

Both implementations use finite-volume methods, Kurganov–Tadmor numerical fluxes and second-order time stepping. Tests reproduce viscous broadening of shock profiles and suppression of Kelvin–Helmholtz vortices as viscosity increases. Bjorken and Gubser flows provide comparisons with symmetry-reduced BDNK solutions; refinement studies show approximately second-order convergence in the reported benchmarks.

Fewer evolved variables did not automatically mean better numerical performance. In the Gubser tests, the mixed-order implementation required smoothing near the outer boundary and showed larger errors at comparable resolutions. Its stability also depended strongly on the reconstruction scheme. The first-order implementation generally performed better in those tests.

On fixed Schwarzschild and Kerr backgrounds, the team initialized radial accretion and a rotating torus using ideal-fluid configurations, then examined departures under finite viscosity. Boundary effects limited the accretion evolution. The torus required artificial interface smoothing, and its comparatively large changes may partly reflect limitations of the numerical setup. Neither test includes evolving spacetime or fluid self-gravity.

Why two viscous theories diverge

For heavy-ion applications, the researchers used an equation of state connecting lattice-QCD results to a hadron resonance gas, with Gaussian or fluctuating TRENTo initial profiles. They compared BDNK with a minimal Israel–Stewart model, initializing both with the same full energy-momentum tensor.

The descriptions agree closely at small shear viscosity. At larger viscosity, differences grow, especially with early initialization and rapid longitudinal expansion. The authors connect this behavior to Israel–Stewart’s independently evolving dissipative stresses, which respond over relaxation times, unlike BDNK’s derivative-defined stresses.

This comparison does not establish which theory better describes experiments. It concerns a particular Israel–Stewart variant and selected initial conditions. General initialization from an arbitrary target energy-momentum tensor remains unresolved, and the work does not compare final observables with collision data. Its contribution is a mathematically supported computational framework, with demonstrated benchmarks and explicitly identified limits.