### 4.3 Relativistic beam scheme

Sanders and Prendergast [254] proposed an explicit scheme to solve the equilibrium limit of the
non-relativistic Boltzmann equation, i.e., the Euler equations of Newtonian fluid dynamics. In their
so-called beam scheme the Maxwellian velocity distribution function is approximated by several Dirac delta
functions or discrete beams of particles in each computational cell, which reproduce the appropriate
moments of the distribution function. The beams transport mass, momentum, and energy into adjacent
cells, and their motion is followed to first-order accuracy. The new (i.e., time advanced) macroscopic
moments of the distribution function are used to determine the new local non-relativistic Maxwell
distribution in each cell. The entire process is then repeated for the next time step. The CFL stability
condition requires that no beam of gas travels farther than one cell in one time step. This beam scheme,
although being a particle method derived from a microscopic kinetic description, has all the desirable
properties of modern characteristic-based wave propagating methods based on a macroscopic continuum
description.
The non-relativistic scheme of Sanders and Prendergast [254] has been extended to relativistic flows by
Yang et al. [303]. They replaced the Maxwellian distribution function by its relativistic analogue, i.e., by
the more complex Jüttner distribution function, which involves modified Bessel functions. For
three-dimensional flows the Jüttner distribution function is approximated by seven delta functions or
discrete beams of particles, which can viewed as dividing the particles in each cell into seven distinct
groups. In the local rest frame of the cell these seven groups represent particles at rest and particles moving
in , , and directions, respectively.

Yang et al. [303] show that the integration scheme for the beams can be cast into the form of an upwind
conservation scheme in terms of numerical fluxes. They further show that the beam scheme not only splits
the state vector but also the flux vectors, and has some entropy-satisfying mechanism embedded as
compared with an approximate relativistic Riemann solver [74, 257] based on Roe’s method [248]. The
simplest relativistic beam scheme is only first-order accurate in space, but can be extended to higher-order
accuracy in a straightforward manner. Yang et al. consider three high-order accurate variants (TVD2,
ENO2, ENO3) generalizing their approach developed in [304, 305] for Newtonian gas dynamics, which is
based on the essentially non-oscillatory (ENO) piecewise polynomial reconstruction scheme of Harten et
al. [121].

Yang et al. [303] present several numerical experiments including relativistic one-dimensional shock tube
flows and the simulation of relativistic two-dimensional Kelvin–Helmholtz instabilities. The shock tube
experiments consist of a mildly relativistic shock tube, relativistic shock heating of a cold flow, the
relativistic blast wave interaction of Woodward and Colella [300] (see Section 6.2.3), and the perturbed
relativistic shock tube flow of Shu and Osher [261].