# An exact Riemann Solver for multidimensional special relativistic hydrodynamics

## 1 Introduction

The decay of a discontinuity separating two constant initial states (Riemann problem) has played a very important role in the development of numerical hydrodynamic codes in classical (Newtonian) hydrodynamics after the pioneering work of Godunov . Nowadays, most modern high-resolution shock-capturing methods  are based on the exact or approximate solution of Riemann problems between adjacent numerical cells and the development of efficient Riemann solvers has become a research field in numerical analysis in its own (see, e.g., ).

Riemann solvers began to be introduced in numerical relativistic hydrodynamics at the beginning of the nineties and, presently, the use of high-resolution shock-capturing methods based on Riemann solvers is considered as the best strategy to solve the equations of relativistic hydrodynamics in nuclear physics (heavy ion collisions) and astrophysics (stellar core collapse, supernova explosions, extragalactic jets, gamma-ray bursts). This fact has caused a rapid development of Riemann solvers for both special and general relativistic hydrodynamics .

The main idea behind the solution of a Riemann problem (defined by two constant initial states, $`L`$ and $`R`$, at left and right of their common contact surface) is that the self-similarity of the flow through rarefaction waves and the Rankine-Hugoniot relations across shocks allow one to connect the intermediate states $`I_{}`$ ($`I=L,R`$) with their corresponding initial states, $`I`$. The analytical solution of the Riemann problem in classical hydrodynamics  rests on the fact that the normal velocity in the intermediate states, $`v_I_{}^n`$, can be written as a function of the pressure $`p_I_{}`$ in that state (and the flow conditions in state $`I`$). Thus, once $`p_I_{}`$ is known, $`v_I_{}^n`$ and all other unknown state quantities of $`I_{}`$ can be calculated.

In the case of relativistic hydrodynamics the same procedure can be followed , the major difference with classical hydrodynamics stemming from the role of tangential velocities. While in the classical case the decay of the initial discontinuity does not depend on the tangential velocity (which is constant across shock waves and rarefactions), in relativistic calculations the components of the flow velocity are coupled through the presence of the Lorentz factor.

## 2 The equations of relativistic hydrodynamics

The equations of relativistic hydrodynamics admit a conservative formulation which has been exploited in the last decade to implement high-resolution shock-capturing methods. In Minkowski space time the equations in this formulation read

$$_t𝐔+_i𝐅^{(i)}=0$$
(1)   
where $`𝐔`$ and $`𝐅^{(i)}(𝐔)`$ ($`i=1,2,3`$) are, respectively, the vectors of conserved variables and fluxes

$$𝐔=(D,S^1,S^2,S^3,\tau )^T$$
(2)   
$$𝐅^{(i)}=(Dv^i,S^1v^i+p\delta ^{1i},S^2v^i+p\delta ^{2i},S^3v^i+p\delta ^{3i},S^iDv^i)^T.$$
(3)   
The conserved variables (the rest-mass density, $`D`$, the momentum density, $`S^i`$, and the energy density $`\tau `$) are defined in terms of the primitive variables,$`(\rho ,v^i,\epsilon )`$, according to

$`D=\rho W,S^i=\rho hW^2v^i,\tau =\rho hW^2pD`$ (4)   
where $`W=(1v^2)^{1/2}`$ is the Lorentz factor and $`h=1+\epsilon +p/\rho `$ the specific enthalpy.

In the following we shall restrict our discussion to an ideal gas equation of state with constant adiabatic exponent, $`\gamma `$, for which the specific internal energy is given by

$$\epsilon =\frac{p}{(\gamma 1)\rho }.$$
(5)   
## 3 Relation between the normal flow velocity and pressure behind relativistic rarefaction waves

Choosing the surface of discontinuity to be normal to the $`x`$-axis, rarefaction waves are self-similar solutions of the flow equations depending only on the combination $`\xi =x/t`$. Getting rid of all the terms with $`y`$ and $`z`$ derivatives in equations (1) and substituting the derivatives of $`x`$ and $`t`$ in terms of the derivatives of $`\xi `$, the system of equations can be reduced to just one ordinary differential equation (ODE) and two algebraic conditions

$`\rho hW^2(v^x\xi )dv^x+(1\xi v^x)dp=0`$ (6)   
$`hWv^y=\text{constant}`$ (7)   
$`hWv^z=\text{constant},`$ (8)   
with $`\xi `$ constrained by

$`\xi `$ $`=`$ $`{\displaystyle \frac{v^x(1c_s^2)\pm c_s\sqrt{(1v^2)[1v^2c_s^2(v^x)^2(1c_s^2)]}}{1v^2c_s^2}},`$ (9)   
because non-trivial similarity solutions exist only if the Wronskian of the original system vanish. We have denoted by $`c_s`$ the speed of sound, provided by the equation of state. The plus and minus sign correspond to rarefaction waves propagating to the right ($`_{}`$) and left ($`_{}`$), respectively. The two solutions for $`\xi `$ correspond to the maximum and minimum eigenvalues of the Jacobian matrix associated to $`𝐅^{(x)}(𝐔)`$ , generalizing the result found for vanishing tangential velocity .

From equations (7) and (8) it follows that $`v^y/v^z=`$ constant, i.e., the tangential velocity does not change direction across rarefaction waves. Notice that, in a kinematical sense, the Newtonian limit ($`v^i1`$) leads to $`W=1`$, but equations (7) and (8) do not reduce to the classical limit $`v^{y,z}=`$ constant, because the specific enthalpy still couples the tangential velocities. Thus, even for slow flows, the Riemann solution presented in this paper must be employed for thermodynamically relativistic situations ($`h1`$). The same result can be deduced from the Rankine-Hugoniot relations for shock waves (see next section).

Using (9) and the definition of the sound speed, the ODE (6) can be written as

$$\frac{dv^x}{dp}=\pm \frac{c_s}{W^2\gamma p}\frac{1}{\sqrt{1+g(\xi _\pm ,v^x,v^t)}}$$
(10)   
where $`v^t=\sqrt{(v^y)^2+(v^z)^2}`$ is the absolute value of the tangential velocity and

$$g(\xi _\pm ,v^x,v^t)=\frac{(v^t)^2(\xi _\pm ^21)}{(1\xi _\pm v^x)^2}.$$
(11)   
Considering that in a Riemann problem the state ahead of the rarefaction wave is known, equation (10) can be integrated with the constraint $`hWv^t=`$ constant, allowing to connect the states ahead ($`a`$) and behind ($`b`$) the rarefaction wave. The solution is only a function of $`p_b`$ and can be stated in compact form as

$$v_b^x=_{}^a(p_b).$$
(12)   
## 4 Relation between post-shock flow velocities and pressure for relativistic shock waves.

The Rankine-Hugoniot conditions relate the states on both sides of a shock and are based on the continuity of the mass flux and the energy-momentum flux across shocks. Their relativistic version was first obtained by Taub  (see also ).

Considering the surface of discontinuity as normal to the $`x`$-axis, the invariant mass flux across the shock can be written as

$$jW_sD_a(V_sv_a^x)=W_sD_b(V_sv_b^x).$$
(13)   
where $`V_s`$ is the coordinate velocity of the hyper-surface that defines the position of the shock wave and $`W_s`$ is the correspondent Lorentz factor, $`W_s=(1V_s^2)^{1/2}`$. According to our definition, $`j`$ is positive for shocks propagating to the right. In terms of the mass flux, $`j`$, the Rankine-Hugoniot conditions are

$$[v^x]=\frac{j}{W_s}\left[\frac{1}{D}\right],$$
(14)   
$$[p]=\frac{j}{W_s}\left[\frac{S^x}{D}\right],$$
(15)   
$$\left[hWv^y\right]=0,$$
(16)   
$$\left[hWv^z\right]=0,$$
(17)   
$$[v^xp]=\frac{j}{W_s}\left[\frac{\tau }{D}\right].$$
(18)   
Equations (16) and (17) imply that the quantity $`hWv^{y,z}`$ is constant across a shock wave and, hence, that the orientation of the tangential velocity does not change. The latter result also holds for rarefaction waves (see §3). Equations (14), (15) and (18) can be manipulated to obtain $`v_b^x`$ as a function of $`p_b`$, $`j`$ and $`V_s`$. Using the relation $`S^x=(\tau +p+D)v^x`$ and after some algebra, one finds

$$v_b^x=\left(h_aW_av_a^x+\frac{W_s(p_bp_a)}{j}\right)\left(h_aW_a+(p_bp_a)\left(\frac{W_sv_a^x}{j}+\frac{1}{\rho _aW_a}\right)\right)^1.$$
(19)   
The final step is to express $`j`$ and $`V_s`$ as a function of the post-shock pressure. First, from the definition of the mass flux we obtain

$$V_s^\pm =\frac{\rho _a^2W_a^2v_a^x\pm |j|\sqrt{j^2+\rho _a^2W_a^2(1v_{a}^{x}{}_{}{}^{2})}}{\rho _a^2W_a^2+j^2}$$
(20)   
where $`V_s^+`$ ($`V_s^{}`$) corresponds to shocks propagating to the right (left).

Second, from the Rankine-Hugoniot relations and the physical solution of $`h_b`$ obtained from the Taub adiabat  (the relativistic version of the Hugoniot adiabat), that relates only thermodynamic quantities on both sides of the shock, the square of the mass flux $`j^2`$ can be obtained as a function of $`p_b`$. Using the positive (negative) root of $`j^2`$ for shock waves propagating towards the right (left), the desired relation between the post-shock normal velocity $`v_b^x`$ and the post-shock pressure $`p_b`$ is obtained. In a compact way the relation reads

$$v_b^x=𝒮_{}^a(p_b).$$
(21)   
We refer to the interested reader to references  for further details.

## 5 The solution of the Riemann problem with arbitrary tangential velocities.

The time evolution of a Riemann problem can be represented as:

$$IL𝒲_{}L_{}𝒞R_{}𝒲_{}R$$
(22)   
where $`𝒲`$ denotes a simple wave (shock or rarefaction), moving towards the initial left ($``$) or right ($``$) states. Between them, two new states appear, namely $`L_{}`$ and $`R_{}`$, separated from each other through the third wave $`𝒞`$, which is a contact discontinuity. Across the contact discontinuity pressure and normal velocity are constant, while the density and the tangential velocity exhibits a jump.

As in the Newtonian case, the compressive character of shock waves (density and pressure rise across the shock) allows us to discriminate between shocks ($`𝒮`$) and rarefaction waves ($``$):

$$𝒲_{()}=\{\begin{array}{ccc}\hfill _{()}& ,& p_bp_a\hfill \\ \hfill 𝒮_{()}& ,& p_b>p_a\hfill \end{array}$$
(23)   
where $`p`$ is the pressure and subscripts $`a`$ and $`b`$ denote quantities ahead and behind the wave. For the Riemann problem $`aL(R)`$ and $`bL_{}(R_{})`$ for $`𝒲_{}`$ and $`𝒲_{}`$, respectively.

The solution of the Riemann problem consists in finding the positions of the waves separating the four states and the intermediate states, $`L_{}`$ and $`R_{}`$. The functions $`𝒲_{}`$ and $`𝒲_{}`$ allow one to determine the functions $`v_R^x(p)`$ and $`v_L^x(p)`$, respectively. The pressure $`p_{}`$ and the flow velocity $`v_{}^x`$ in the intermediate states are then given by the condition

$$v_R^x(p_{})=v_L^x(p_{})=v_{}^x.$$
(24)   
which is an implicit algebraic equation in $`p_{}`$ and can be solved by means of an iterative method. When $`p_{}`$ and $`v_{}^x`$ have been obtained, the equation of state gives the specific internal energy and the remaining state variables of the intermediate state $`I_{}`$ can be calculated using the relations between $`I_{}`$ and the respective initial state $`I`$ given through the corresponding wave.

Notice that the solution of the Riemann problem depends on the modulus of $`v^t`$, but not on the direction of the tangential velocity. Figure 1 shows the solution of a particular Riemann problem  for different values of the tangential velocity $`v^y=0,0.5,0.9,0.99`$. The crossing point of any two lines in the upper panel gives the pressure and the normal velocity in the intermediate states. Whereas the pressure in the intermediate state can take any value between $`p_L`$ and $`p_R`$, the normal flow velocity can be arbitrarily close to zero in the case of an extremely relativistic tangential flow.

To study the influence of tangential velocities on the solution a Riemann problem, we have calculated the solution of a standard test involving the propagation of a relativistic blast wave produced by a large jump in the initial pressure distribution for different combinations of tangential velocities . Although the structure of the solution remains unchanged for different tangential velocities, the values in the constant state may change by a large amount.

## 6 Discussion and conclusions.

We have obtained the exact solution of the Riemann problem in special relativistic hydrodynamics with arbitrary tangential velocities. Unlike in Newtonian hydrodynamics, tangential velocities are coupled with the rest of variables through the Lorentz factor, present in all terms in all equations. It strongly affects the solution, especially for ultra-relativistic tangential flows. In addition, the specific enthalpy also acts as a coupling factor and modifies the solution for the tangential velocities in thermodynamically relativistic situations (energy density and pressure comparable to or larger than the proper rest-mass density), rendering the classical solution incorrect in slow flows with very large internal energies.

Our solution has interesting practical applications. First, it can be used to test the different approximate relativistic Riemann solvers and the multi-dimensional hydrodynamic codes based on directional splitting. Second, it can be used to construct multi-dimensional relativistic Godunov-type hydro codes. As an example, we have simulated a relativistic tube , in a $`100\times 100`$ Cartesian grid, where the initial discontinuity was located in a main diagonal. The initial states were $`\rho _L=10`$, $`\rho _R=1`$, $`p_L=13.3`$, $`p_R=0.66E3`$, $`v_L=0`$, $`v_R=0`$, and the adiabatic index is $`\gamma =5/3`$. Spatial order of accuracy was set to second order by means of a monotonic piecewise linear reconstruction procedure and second order in time is obtained by using a Runge-Kutta method for time advancing. The exact solution of the Riemann problem is used at every interface to calculate the numerical fluxes. The results are shown in Figure 2, and are comparable to those obtained with other HRSC methods. Profiles of all variables are stable and discontinuities are well resolved without excessive smearing. An efficient implementation of this exact Riemann solver in the context of multidimensional relativistic PPM  is in progress and will be reported elsewhere.