PAPER • OPEN ACCESS

# Non-linear Eulerian hydrodynamics of dark energy: Riemann problem and finite volume schemes

To cite this article: Linda Blot et al JCAP05(2023)001

View the [article online](#) for updates and enhancements.

You may also like

- Soliton gas in integrable dispersive hydrodynamics
  Gennady A EI

- Effects of Tangential Velocity in the Reactive Relativistic Riemann Problem
  A. Harpole and I. Hawke

- On viscous limit solutions of the Riemann problem for the equations of isentropic gas dynamics in Eulerian coordinates
  B P Andreianov
# Non-linear Eulerian hydrodynamics of dark energy: Riemann problem and finite volume schemes

Linda Blot,ª Pier Stefano Corasanitiᵇ and Fabian Schmidtª

E-mail: lblot@mpa-garching.mpg.de, pier-stefano.corasaniti@obspm.fr, fabians@mpa-garching.mpg.de

Received October 11, 2022

Revised March 28, 2023

Accepted April 13, 2023

Published May 2, 2023

**Abstract.** Upcoming large-scale structure surveys can shed new light on the properties of dark energy. In particular, if dark energy is a dynamical component, it must have spatial perturbations. Their behaviour is regulated by the speed of sound parameter, which is currently unconstrained. In this work, we present the numerical methods that will allow to perform cosmological simulations of inhomogeneous dark energy scenarios where the speed of sound is small and non-vanishing. We treat the dark energy component as an effective fluid and build upon established numerical methods for hydrodynamics to construct a numerical solution of the effective continuity and Euler equations. In particular, we develop conservative finite volume schemes that rely on the solution of the Riemann problem, which we provide here in both exact and approximate forms for the case of a dark energy fluid.

**Keywords:** cosmological simulations, dark energy theory

**ArXiv ePrint:** 2210.04800

<figure><img src="image_1.png" /></figure>

<footnote>ªMax Planck Institute for Astrophysics,
Karl-Schwarzschild str. 1, 85748 Garching, Germany</footnote>

<footnote>ᵇLaboratoire Univers et Théories (LUTh), UMR 8102 CNRS,
Observatoire de Paris, Université Paris Diderot,
5 Place Jules Janssen, 92190 Meudon, France</footnote>
# Contents

1 Introduction 1

2 Euler equations for dark energy fluids 3

3 The Riemann problem 5

3.1 Structure of the solution in 1D 6

3.2 Structure of the solution in 3D 9

4 Riemann solvers 11

4.1 Exact Riemann solver 11

4.2 Approximate Riemann solvers 14

4.2.1 Primitive variables Riemann solver 14

4.2.2 Harten, Lax and van Leer Riemann solver 16

5 Finite volume conservative methods 19

5.1 MUSCL-Hancock method 20

5.2 Piecewise Linear Method (PLM) 22

5.3 Piecewise Parabolic Method (PPM) 23

5.4 Advection test for higher-order schemes 25

6 Conclusions 26

A Elementary wave solutions in 1D 28

A.1 Rarefaction waves 28

A.1.1 Left rarefaction 28

A.1.2 Right rarefaction 29

A.2 Shock waves 29

A.2.1 Left shock 29

A.2.2 Right shock 30

B Elementary wave solutions in 3D 30

B.1 Rarefaction waves 30

B.2 Shock waves 31

# 1 Introduction

We have compelling evidence that the universe is currently undergoing an accelerated phase of expansion. In the standard cosmological model, this requires the existence of an exotic component, dubbed dark energy, that is characterized by negative pressure and dominates the present-day cosmic energy budget. A simple cosmological constant term (Λ) in Einstein's equation of General Relativity can play the role of dark energy since it behaves as a negative pressure component. Together with the Cold Dark Matter (CDM) hypothesis, the ΛCDM scenario currently reproduces all cosmological observations on the scales thus far available [1-4]. Nevertheless, the success of this model remains purely empirical as the physical origin of the cosmic dark components remains unknown.
Observations indicate that the energy density associated with $\Lambda$ is $\rho_{\Lambda} \approx 10^{-47}$ GeV$^4$. It is the smallness of this value that has so far posed puzzling problems to any attempt of finding a satisfying theoretical explanation of $\Lambda$ (for a review see [5]). On the observational side, the $\Lambda$CDM model is showing its limitations in the form of tensions between cosmological parameter values recovered from different kinds of cosmological probes [6, 7]. For these reasons, alternative hypotheses to the cosmological constant have been advanced. Beyond the cosmological constant, dark energy may result either from a modification of gravity on cosmic scales [8, 9] or be the effect of a new degree of freedom not accounted for by the Standard Model of particles [10].

Given the lack of new fundamental principles which may guide the dark energy theoretical model building, there is hope nevertheless that the next generation of cosmological observations may reveal hints of the new physics underlying the dark energy phenomenon. Ongoing and upcoming large-scale structure surveys have the potential to constrain dark energy properties using established probes such as galaxy clustering, cosmic shear, and galaxy cluster number counts as well as newly proposed probes such as voids, counts in cells, and others. The usual approach is to model dark energy as a perfect fluid characterized by homogeneous pressure and density that are related through an equation of state parameter, $w$, the cosmological constant case corresponding to $w = -1$. However, if dark energy is dynamical, i.e. its energy density is not constant in time, general covariance imposes that it has to vary in space as well. We then have an additional macroscopic quantity that specifies the properties of the dark energy fluid, namely the speed of sound of dark energy perturbations, $c_s$. This suggests that the detection of the clustering of dark energy through observations of the cosmic structures may provide smoking gun evidence of the dynamical nature of the dark energy phenomenon.

Bounds on the dark energy equation of state and speed of sound can be translated into constraints on the microscopic properties of dark energy since these quantities can be related to terms in a Lagrangian effective field theory of dark energy (see e.g. [11, 12]). For instance, quintessence models in which dark energy results from the dynamics of a minimally coupled scalar field are characterized by $c_s$ equal to the speed of light $c$ [13], while k-essence models in which dark energy is due to a scalar field with a non-canonical kinetic term have $c_s^2 \approx 0$ [14]. It is important to notice that a perfect fluid approach is not limited to the phenomenological study of dark energy. For instance models of unified dark matter and dark energy [15] as well as alternative Dark Matter scenarios (see e.g. [16]) are described as dark fluids characterized by an equation of state and a speed of sound. For this reason throughout the text we may refer indistinctly to Dark Energy or dark fluids.

The clustering properties of Dark Energy can be tested through observations of the cosmic distribution of matter at different scales and times. On large scales, linear perturbation theory can be used to predict the effect of dark energy inhomogeneities on the evolution of matter density perturbations [17–19]. However, since dark energy is a late-time phenomenon, the impact of dark energy perturbations on those scales is small and well within cosmic variance uncertainties, so that current data remains uninformative [19, 20].

This may not be the case at small scales and late times where the collapse process is highly non-linear. For quintessence-like fluids, dark energy perturbations on sub-horizon scales have been damped by free-streaming, so that only the background expansion is affected. Thus, in these models the non-linear clustering of matter can be investigated using standard N-body simulations [21–23]. Instead, in models with $c_s < c$ perturbations may cluster at small scales so that it is necessary to include them in the simulation to capture the fully non-linear process of structure collapse in the presence of dark energy inhomogeneities. Previous studies looked
at the fully clustered case $c_s^2 = 0$ in the context of the spherical collapse model [24–26] or higher-order perturbation theory [27–29]. Recently the first cosmological simulation code that includes the effects of clustering dark energy has been developed, k-evolution [30, 31]. This is a relativistic N-body code that treats the dark energy component as a scalar field. While it has the advantage of including all the relevant relativistic effects, it cannot deal with the large dynamic range of typical high-resolution cosmological simulations since it employs a fixed grid that cannot be refined. Moreover, ref.s [32, 33] reported instabilities in the scalar field solver in the context of cosmological simulations, providing strong motivation for an independent check of whether these instabilities are a physical property of clustering dark energy.

The work presented here is complementary to k-evolution in that it uses the Newtonian approximation but treats dark energy as a fluid and is meant to be implemented in fully non-linear hydrodynamics codes. In fact, we will show that the non-linear evolution of dark energy fluctuations is described by hyperbolic Euler equations which can be cast in the form of modified equations of hydrodynamics. Numerical schemes that are best suited to solve this type of equations are based on Finite Volume (FV) methods. Many of these schemes are already implemented in cosmological simulation codes such as ENZO [34], FLASH [35] or RAMSES [36] to follow the cosmological collapse of baryonic gas. These codes use high-order variants of the Godunov scheme originally introduced in [37] and have been shown to provide accurate solutions to a variety of problems in gas dynamics. The accuracy of these schemes relies on their ability to capture discontinuities in the fluid flow and correctly predict the velocity of propagating waves. The starting point of these methods is the solution to the Riemann problem. Exact and approximate solutions have been derived for several systems including real gases [38], inviscid flows of perfect gases [39], gases with generic equations of state [40] and compressible liquids [41].

Here, we present a detailed study of the Riemann problem for cosmic dark fluid as a first step toward a cosmological study of the non-linear structure formation in inhomogeneous dark energy models. We derive an exact solution to the Riemann problem and construct several approximated solvers that in combination with Godunov-type schemes can efficiently solve the non-linear fluid equations for dark energy.

The paper is organized as follows: in section 2 we present the Euler equations for a clustering dark energy fluid and cast them in the form that allows for an easier numerical implementation. We then present the solution of the Riemann problem for this fluid in section 3 and develop numerical Riemann solvers in section 4. Finally, in section 5 we present working examples of finite volume conservative schemes that can be used to solve the Euler equations for clustering dark energy and test them using standard hydrodynamical test cases. We draw our conclusions in section 6.

## 2 Euler equations for dark energy fluids

The equations describing the evolution of a dark energy (or any non-relativistic) fluid under the influence of Newtonian gravity in a Friedmann-Lemaitre-Robertson-Walker (FLRW) background have been derived in [27], using conformal-Newtonian gauge. In Cartesian comoving coordinates these read as:

$$ \frac{\partial \rho}{\partial \tau} + 3\mathcal{H} \left(\rho + \frac{p}{c^2}\right) + \vec{\nabla} \cdot \left[ \left(\rho + \frac{p}{c^2}\right) \vec{v} \right] = 0 \qquad (2.1) $$

$$ \frac{\partial \vec{v}}{\partial \tau} + \mathcal{H} \vec{v} + (\vec{v} \cdot \vec{\nabla}) \vec{v} = -\frac{1}{\rho + \frac{p}{c^2}} \left( \vec{\nabla} p + \vec{v} \frac{\partial}{\partial \tau} \frac{p}{c^2} \right) - \vec{\nabla} \Phi, \qquad (2.2) $$
where $\rho$ and $p$ are the dark energy energy density and pressure respectively, $\vec{v}$ is the dark energy peculiar velocity with respect to the Hubble flow, $\tau$ is the conformal time defined by $d\tau = dt/a(t)$, with $a$ the scale factor, $H = d\ln a/d\tau$ is the Hubble rate and $\Phi$ is the gravitational potential. The derivation of these equations assumes non-relativistic peculiar velocities ($v \ll c$), and scales much smaller than the horizon, which removes a number of relativistic contributions.

Since we aim to use conservative numerical methods we want to cast the fluid equations in the form that is closest to the usual conservative form for baryonic fluids. To achieve this we can start by re-writing eqs. (2.1) and (2.2) as:

$$ \frac{\partial}{\partial \tau} \left( \rho + \frac{p}{c^2} \right) + 3\mathcal{H} \left( \rho + \frac{p}{c^2} \right) + \vec{\nabla} \cdot \left[ \left( \rho + \frac{p}{c^2} \right) \vec{v} \right] = \frac{\partial p}{\partial \tau} \qquad (2.3) $$

$$ \frac{\partial}{\partial \tau} \left[ \left( \rho + \frac{p}{c^2} \right) \vec{v} \right] + 4\mathcal{H} \left( \rho + \frac{p}{c^2} \right) \vec{v} + \vec{\nabla} \cdot \left[ \left( \rho + \frac{p}{c^2} \right) \vec{v} \otimes \vec{v} \right] = - \left( \rho + \frac{p}{c^2} \right) \vec{\nabla} \Phi - \vec{\nabla} p. \qquad (2.4) $$

Given the non-barotropic effective description of the dark energy fluid, it is necessary to split the evolution of the background density $\bar{\rho}$ and pressure $\bar{p}$ from that of the fluctuations $\delta\rho$ and $\delta p$ such that

$$ \rho(\vec{x}, \tau) \equiv \bar{\rho}(\tau) + \delta\rho(\vec{x}, \tau) \qquad (2.5) $$

$$ p(\vec{x}, \tau) \equiv \bar{p}(\tau) + \delta p(\vec{x}, \tau), \qquad (2.6) $$

where $\bar{p}(\tau) = w(\tau)\bar{\rho}(\tau)c^2$ and the evolution of the background density is governed by

$$ \frac{\partial \bar{\rho}}{\partial \tau} = -3\mathcal{H}(1+w)\bar{\rho}. \qquad (2.7) $$

We then describe the dark energy pressure perturbations by relating them to dark energy density perturbations in the dark energy rest frame, denoted by a subscript rf:

$$ \delta p_{\text{rf}}(\vec{x}, \tau) = c_s^2(\tau) \delta \rho_{\text{rf}}(\vec{x}, \tau), \qquad (2.8) $$

where $c_s(\tau)$ is the effective sound speed, which in general can be time-dependent.

We can now introduce the new variable:

$$ \Pi \equiv 1 + w + \left( 1 + \frac{c_s^2}{c^2} \right) \delta, \qquad (2.9) $$

where $\delta = \delta\rho/\bar{\rho}$. We will restrict here to the case $\Pi > 0$ because the solution of the Riemann problem becomes ill-defined for $\Pi \le 0$ (see e.g. eq. (A.3) and (A.4)), but there is also a physical reason for restricting to this case. In fact, one can also write $\Pi$ as $\Pi = (\rho_{\text{rf}} + p_{\text{rf}}/c^2)/\bar{\rho}_{\text{rf}}$ and imposing that $(\rho + p/c^2)_{\text{rf}} > 0$ is equivalent to enforcing the null-energy condition. Notice that this translates to a floor in the physical dark energy underdensity:

$$ \delta > - \frac{1+w}{1+\frac{c_s^2}{c^2}}, \qquad (2.10) $$

meaning that the total dark energy density can never reach 0 if $w < 0$. Eq. (2.1) and eq. (2.2) can finally be written in the form:

$$ \frac{\partial \Pi}{\partial \tau} + \left( 1 + \frac{c_s^2}{c^2} \right) \vec{\nabla} \cdot \left( \Pi \vec{v} \right) = 3\mathcal{H} \left( w - \frac{c_s^2}{c^2} \right) \left( \Pi - 1 - w \right) \qquad (2.11) $$

$$ \frac{\partial(\Pi\vec{v})}{\partial\tau} + \vec{\nabla} \cdot (\Pi\vec{v} \otimes \vec{v}) + \frac{c_s^2}{\left(1 + \frac{c_s^2}{c^2}\right)} \vec{\nabla}\Pi = \mathcal{H}(3w - 1)\Pi\vec{v} - \Pi\vec{\nabla}\Phi, \qquad (2.12) $$
which we call quasi-conservative in analogy with the conservative form of the Euler equations for baryonic fluids. In fact, in the case of a vanishing equation of state $w$ and sound speed $c_s$, eqs. (2.11) and (2.12) reduce to the standard conservation laws of mass and momentum. Here, we have intentionally written the friction terms due to the cosmic expansion on the right-hand side to distinguish the advection part of the Euler equations from non-advectioning time-dependent source terms. We have also assumed that $w$ and $c_s$ are constant in time. Relaxing this assumption will introduce additional source terms. However, even if $c_s$ varies in time, this does not alter the form of the advection part of the Euler equations. Hence, the derivations presented hereafter remain valid also in this case.

This system of equations gives a description of the gravitational dynamics of dark energy perturbations in the presence of other matter components provided a closure equation for the gravitational potential $\Phi$. This is given by the Poisson equation [see e.g. 27]:

$$ \nabla^2 \Phi = 4\pi G a^2 \left( \rho_{\text{tot}} + 3 \frac{p_{\text{tot}}}{c^2} \right) \qquad (2.13) $$

where $\rho_{\text{tot}}$ and $p_{\text{tot}}$ are the total energy density and pressure of all gravitating components. Similar to what is done in the case of cosmological hydrodynamical simulations of dark matter and baryons, the idea of simulating the clustering of dark matter and dark energy is to use the N-body method to solve the dynamics of dark matter particles and the hydrodynamical description given by eq. (2.11) and (2.12) to solve that of dark energy fluctuations, where the same total potential, sourced by matter and dark energy fluctuations, is used for both components. These equations can be implemented for example in cosmological N-body/hydro codes with an Adaptive Mesh Refinement (AMR) of the simulation volume in which the energy density and pressure of dark energy can be added to that of the matter density field on the same AMR grid before numerically solving the Poisson equation.

### 3 The Riemann problem

The advection parts of eqs. (2.11) and (2.12) form a hyperbolic system of partial differential equations that can be solved numerically using finite volume methods. Here, we will focus on the Godunov method originally introduced in [37]. At the heart of the method, as we will see in detail in section 5, is the solution of the Riemann Problem (RP).

Let us briefly introduce here what is the Riemann Problem and how the solution is found. We can start by writing eqs. (2.11) and (2.12) in a compact state-vector form in Cartesian coordinates:

$$ \mathbf{U}_{\tau} + \mathbf{F}(\mathbf{U})_{x} + \mathbf{G}(\mathbf{U})_{y} + \mathbf{H}(\mathbf{U})_{z} = \mathbf{S}(\mathbf{U}), \qquad (3.1) $$

where the indices denote partial derivatives, $\mathbf{U}$ is a state vector of unknowns, $\mathbf{F}(\mathbf{U})$, $\mathbf{G}(\mathbf{U})$ and $\mathbf{H}(\mathbf{U})$ are vectors of fluxes and $\mathbf{S}(\mathbf{U})$ a vector of sources. Let us focus on the advection part of the equations, which is simply given by imposing $\mathbf{S}(\mathbf{U}) = 0$. Physically, this corresponds to the case of a non-expanding background without gravity, which is what we consider in the following. Since in most applications this system is solved along each Cartesian direction separately, let us simply focus on the $x$ direction.

The Riemann problem is an initial value problem (IVP) of the advection equation

$$ \mathbf{U}_{\tau} + \mathbf{F}(\mathbf{U})_{x} = 0, \qquad (3.2) $$
with initial conditions

$$ \mathbf{U}(x, 0) = \begin{cases} \mathbf{U}_L & \text{if } x < 0 \\ \mathbf{U}_R & \text{if } x > 0 \end{cases} \qquad (3.3) $$

In general, the solution of the RP for a $m \times m$ non-linear hyperbolic system consists of $m+1$ constant states separated by $m$ waves in the plane $x-\tau$. These can be either shocks, rarefaction or contact waves associated with the so called characteristic fields of the advection equation eq. (3.2) [see e.g. 42, 43, for a standard textbook presentation]. The characteristics are curves in the plane $x-\tau$ along which eq. (3.2) reduces to an ordinary differential equation. In order to determine these curves let us rewrite eq. (3.2) as

$$ \mathbf{U}_{\tau} + \mathbf{A}(\mathbf{U})\mathbf{U}_{x} = 0, \qquad (3.4) $$

where $\mathbf{A}(\mathbf{U}) = \partial\mathbf{A}(\mathbf{U})/\partial\mathbf{U}$ is the Jacobian matrix of the flux. Finding the eigenvalues $\lambda_i$ and associated right eigenvectors $\mathbf{R}_i$ of the Jacobian we can verify the nature of the waves by computing:

$$ \vec{\nabla}_{\mathbf{U}} \lambda_i \cdot \mathbf{R}_i, \qquad (3.5) $$

where $\vec{\nabla}_{\mathbf{U}}$ indicates the derivative with respect to the state variables $\mathbf{U}$. If $\vec{\nabla}_{\mathbf{U}}\lambda_i \cdot \mathbf{R}_i \neq 0$ the characteristics are said to be genuinely non-linear and the associated waves are either shocks or rarefaction waves, otherwise the wave is a contact discontinuity.

When $\lambda_i$ is associated to a rarefaction wave the states to the left and right of the wave $\mathbf{U}_L$ and $\mathbf{U}_R$ are connected through a smooth transition. The characteristics at the two sides of the wave diverge so that a new region emerges: the rarefaction fan. This region is enclosed between the head and the tail of the fan, whose speeds are given by $s_H = \lambda_i(\mathbf{U}_R)$ and $s_T = \lambda_i(\mathbf{U}_L)$ respectively. This is depicted schematically in figure 1(a).

The solution across rarefaction waves can be found using the generalized Riemann Invariants (see [44] for a detailed discussion). These are ordinary differential equations that relate the changes in the state variables to the respective components of the $i$-th eigenvector:

$$ \frac{du_1}{r_1^i} = \frac{du_2}{r_2^i} = \dots = \frac{du_m}{r_m^i}, \qquad (3.6) $$

where $u_i$ are the elements of the state vector $\mathbf{U}$.

When the wave is a shock the characteristics at the two sides of the wave converge and the left and right states $\mathbf{U}_L$ and $\mathbf{U}_R$ are connected through a single jump-discontinuity that moves with speed $s_s$, as shown in figure 1(b). In this case we can make use of the Rankine-Hugoniot conditions, that relate the state variables on the two sides of the discontinuity:

$$ \Delta\mathbf{F} = s\Delta\mathbf{U}. \qquad (3.7) $$

In case the wave is a contact discontinuity the states to the left and right of the wave are connected through a single jump-discontinuity, as shown in figure 1(c). The characteristics on the two sides run parallel to the wave, which moves with speed $s_c = \lambda_i(\mathbf{U}_L) = \lambda_i(\mathbf{U}_R)$. We can use either the Riemann Invariants or the Rankine-Hugoniot conditions to find the solution across contact waves.

## 3.1 Structure of the solution in 1D

Let us now focus on the specific case of the Dark Energy fluid and find the analytical solution of the Riemann problem. In section 4.1 we will provide the full solution, however it is useful
<figure><img src="image_2.png" /><figcaption>(a) Rarefaction wave: the left and right states are separated by the rarefaction fan, whose head and tail (thick lines) move with speed $s_H$ and $s_T$ respectively.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>(b) Shock wave: the left and right states are separated by a jump-discontinuity (thick line) that moves with speed $s_s$.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>(c) Contact wave: the left and right states are separated by a jump-discontinuity (thick dashed line) that moves with speed $s_c$.</figcaption></figure>

**Figure 1.** Structure of the solution for the three types of waves. The directions of the characteristics on the two sides of the waves are indicated by the arrows.

to first derive a set of equations that relate state variables across different wave patterns. Hereafter, we will extensively follow Toro's textbook [43] to which we refer interested readers for a detailed presentation of numerical methods in fluid dynamics.

The vectors in eq. (3.1) are given by:

$$ \mathbf{U} = \begin{bmatrix} u_1 \\ u_2 \end{bmatrix}, \quad \mathbf{F} = \begin{bmatrix} \left(1 + \frac{c_s^2}{c^2}\right) u_2 \\ \frac{u_2^2}{u_1} + \frac{c_s^2}{1 + \frac{c_s^2}{c^2}} u_1 \end{bmatrix}, \quad \mathbf{S} = \begin{bmatrix} 3\mathcal{H}\left(w - \frac{c_s^2}{c^2}\right) (u_1 - 1 - w) \\ \mathcal{H}(3w - 1)u_2 - u_1 \frac{\partial \Phi}{\partial x} \end{bmatrix}, \quad (3.8) $$

where $u_1 = \Pi$ and $u_2 = \Pi v$ are the conservative state variables. Since this is a 2 × 2 system of equations the solution consists of 3 states separated by 2 waves. The Jacobian matrix reads:

$$ \mathbf{A}(\mathbf{U}) = \frac{\partial \mathbf{F}}{\partial \mathbf{U}} = \begin{bmatrix} 0 & 1 + \frac{c_s^2}{c^2} \\ \frac{c_s^2}{1 + \frac{c_s^2}{c^2}} - v^2 & 2v \end{bmatrix}, \quad (3.9) $$

and its eigenvalues are given by

$$ \lambda_{\pm} = v \pm c_s \sqrt{1 - \frac{v^2}{c^2}}. \quad (3.10) $$

These are real and distinct if and only if $c_s > 0$. In such a case we have $\lambda_-(v) < \lambda_+(v)$ and the system is strictly hyperbolic. Moreover, since we are considering Dark Energy models for which peculiar velocities are non-relativistic ($v \ll c$) the wave structure of the system approximately reproduces that of an isothermal fluid with constant sound speed, for which $\lambda_{\pm} = v \pm c_s$. It is also worth noticing that the condition of hyperbolicity requires a limiting speed for the dark fluid velocity $v \le c$, which is direct consequence of including terms $p/c^2$ in eq. (2.1) and (2.2). The reader may wonder whether in the presence of gravitational interactions the speed of dark energy fluctuations may numerically become ultra-relativistic, thus leading to loss of hyperbolicity of the system. We believe this not to be the case, since these dark
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Structure of the solution of the Riemann problem for a dark energy fluid in 1D: two waves $\lambda_{\pm}$ (thick lines) generate from the discontinuity at $x_0$ and a new state $\mathbf{U}_*$ emerges in between them.</figcaption></figure>

fluids are cold, with negligibly small velocities, and cosmological simulations generally have gravitational potentials $\Phi \ll 1$.

The right eigenvectors associated with $\lambda_{\pm}$ are given by

$$ \mathbf{R}_{+} = \begin{bmatrix} 1 + \frac{c_s^2}{c^2} \\ v + c_s \sqrt{1 - \frac{v^2}{c^2}} \end{bmatrix}, \quad \mathbf{R}_{-} = \begin{bmatrix} 1 + \frac{c_s^2}{c^2} \\ v - c_s \sqrt{1 - \frac{v^2}{c^2}} \end{bmatrix}, \qquad (3.11) $$

while the left eigenvectors read as

$$ \mathbf{L}_{+} = \begin{bmatrix} \frac{1}{2\left(1 + \frac{c_s^2}{c^2}\right)} & \left(1 - \frac{v}{c_s\sqrt{1 - \frac{v^2}{c^2}}}\right), & \frac{1}{2c_s\sqrt{1 - \frac{v^2}{c^2}}} \\ \frac{1}{2\left(1 - \frac{c_s^2}{c^2}\right)} & \left(1 + \frac{v}{c_s\sqrt{1 - \frac{v^2}{c^2}}}\right), & -\frac{1}{2c_s\sqrt{1 - \frac{v^2}{c^2}}} \end{bmatrix}, \qquad (3.12) $$

and satisfy the relation $\mathbf{L}_{\pm} \cdot \mathbf{R}_{\pm} = \delta_{\pm\pm}$, where $\delta_{\pm\pm}$ is the Kronecker symbol.

Using eqs. (3.10) and (3.11) we obtain

$$ \vec{\nabla}_{\mathbf{U}} \lambda_{\pm} \cdot \mathbf{R}_{\pm} = \pm \frac{c_s}{\Pi \sqrt{1 - \frac{v^2}{c^2}}} \left[ 1 + \frac{v^2}{c^2} \left( 1 + \frac{c_s^2}{c^2} \right) \right], \qquad (3.13) $$

hence the characteristic fields are genuinely non-linear if and only if $c_s > 0$. In this case, the waves associated with $\lambda_{\pm}$ are either rarefaction or shock waves.

In figure 2 we sketch the structure of the solution to the Riemann problem. Given the data to the left $\mathbf{U}_L$ and the right $\mathbf{U}_R$ of the initial discontinuity at $x = x_0$, we want to determine the state $\mathbf{U}^*$ in the region enclosed between the left and right waves and derive the solution of the Riemann problem $\mathbf{U}(x, \tau)$. Once the nature of the left and right waves is known we can determine the value of $\mathbf{U}^*$ as a function of the left or right state using the Generalized Riemann Invariants or the Rankine-Hugoniot conditions.

Notice that eq. (3.2) can be rewritten also in terms of the primitive state variables $\mathbf{W} = (\Pi, v)$:

$$ \mathbf{W}_{\tau} + \tilde{\mathbf{A}}(\mathbf{W})\mathbf{W}_{x} = 0, \qquad (3.14) $$
where

$$ \tilde{\mathbf{A}}(\mathbf{W}) = \begin{bmatrix} (1 + \frac{c_s^2}{c^2}) v & (1 + \frac{c_s^2}{c^2}) \Pi \\ \frac{c_s^2}{\Pi} \left( \frac{1}{1 + \frac{c_s^2}{c^2}} - \frac{v^2}{c^2} \right) & (1 - \frac{c_s^2}{c^2}) v \end{bmatrix}. \quad (3.15) $$

Since conservative and primitive variables are related by a linear transformation $d\mathbf{U} = \Lambda d\mathbf{W}$ with

$$ \Lambda = \begin{bmatrix} 1 & 0 \\ v & \Pi \end{bmatrix}, \quad \Lambda^{-1} = \begin{bmatrix} 1 & 0 \\ -\frac{v}{\Pi} & \frac{1}{\Pi} \end{bmatrix} \quad (3.16) $$

it is easy to show that $\mathbf{A}(\mathbf{U}) = \Lambda \tilde{\mathbf{A}}(\mathbf{W}) \Lambda^{-1}$. This implies that the Jacobian matrix, whether written in terms of primitive variables or conservative ones, has identical eigenvalues, while the corresponding eigenvectors are related by a linear transformation. Therefore, the Riemann problem admits the same wave structure in the two formulations.

In this case the right eigenvectors are given by

$$ \tilde{\mathbf{R}}_+ = \begin{bmatrix} (1 + \frac{c_s^2}{c^2}) \Pi \\ -\frac{c_s^2}{c^2} v + c_s \sqrt{1 - \frac{v^2}{c^2}} \end{bmatrix}, \quad \tilde{\mathbf{R}}_- = \begin{bmatrix} (1 + \frac{c_s^2}{c^2}) \Pi \\ -\frac{c_s^2}{c^2} v - c_s \sqrt{1 - \frac{v^2}{c^2}} \end{bmatrix}, \quad (3.17) $$

while the left eigenvectors read as

$$ \begin{aligned} \tilde{\mathbf{L}}_+ &= \left[ \frac{1}{2 \left(1 + \frac{c_s^2}{c^2}\right) \Pi} \left( 1 + \frac{c_s^2 v}{c^2 \sqrt{1 - \frac{v^2}{c^2}}} \right), \frac{1}{2 c_s \sqrt{1 - \frac{v^2}{c^2}}} \right], \\ \tilde{\mathbf{L}}_- &= \left[ \frac{1}{2 \left(1 + \frac{c_s^2}{c^2}\right) \Pi} \left( 1 - \frac{c_s^2 v}{c^2 \sqrt{1 - \frac{v^2}{c^2}}} \right), \frac{-1}{2 c_s \sqrt{1 - \frac{v^2}{c^2}}} \right] \end{aligned} \quad (3.18) $$

In appendix A we derive the relations between the state variables across different wave patterns.

## 3.2 Structure of the solution in 3D

In the previous section we discussed the Riemann problem in 1D. As we will show here this provides the basis to solve the Riemann problem in 3D. In fact, let us consider eqs. (2.11) and (2.12) without source terms in Cartesian coordinates. In terms of the conservative variable these can be written in a state-vector form:

$$ \mathbf{U} + \mathbf{A}(\mathbf{U})\mathbf{U}_x + \mathbf{B}(\mathbf{U})\mathbf{U}_y + \mathbf{C}(\mathbf{U})\mathbf{U}_z = 0, \quad (3.19) $$
with $\mathbf{U} = (\Pi, \Pi v_x, \Pi v_y, \Pi v_z)$ where $v_x, v_y$ and $v_z$ are the components of the peculiar velocity vector along $x, y$ and $z$ respectively, and

$$ 
\begin{aligned}
\mathbf{A}(\mathbf{U}) &= \begin{bmatrix}
0 & 1 + \frac{c_s^2}{c^2} & 0 & 0 \\
\frac{c_s^2}{1+\frac{c_s^2}{c^2}} - v_x^2 & 2v_x & 0 & 0 \\
-v_x v_y & v_y & v_x & 0 \\
-v_x v_z & v_z & 0 & v_x
\end{bmatrix}, \quad
\mathbf{B}(\mathbf{U}) &= \begin{bmatrix}
0 & 0 & 1 + \frac{c_s^2}{c^2} & 0 \\
-v_x v_y & v_y & v_x & 0 \\
\frac{c_s^2}{1+\frac{c_s^2}{c^2}} - v_y^2 & 2v_y & 0 & 0 \\
-v_y v_z & v_z & 0 & v_y
\end{bmatrix}, \\
\mathbf{C}(\mathbf{U}) &= \begin{bmatrix}
0 & 0 & 0 & 1 + \frac{c_s^2}{c^2} \\
-v_x v_z & v_z & 0 & v_x \\
-v_y v_z & 0 & v_z & v_y \\
\frac{c_s^2}{1+\frac{c_s^2}{c^2}} - v_z^2 & 0 & 0 & 2v_z
\end{bmatrix}.
\end{aligned}
$$

The eigenvalues of these matrices are:

$$ 
\begin{aligned}
\lambda_1 &= v - c_s \sqrt{1 - \frac{v^2}{c^2}} \\
\lambda_{2,3} &= v \\
\lambda_4 &= v + c_s \sqrt{1 - \frac{v^2}{c^2}}
\end{aligned}
$$

with $v = v_x, v_y$ and $v_z$ for $\mathbf{A}, \mathbf{B}$ and $\mathbf{C}$ respectively. Hence, the 3D case can be seen as the composition of three advections with the same structure of the Riemann problem.

Let us consider the advection equations along the x-direction. The eigenvalues $\lambda_{1,4}$ are real and distinct and associated to the genuinely non-linear characteristic fields that correspond to shock and rarefaction waves, while $\lambda_{2,3}$ are linearly degenerate and describe a contact discontinuity. Moreover, the eigenvalues $\lambda_{1,4}$ are identical to the 1D case, so that the solution for $\Pi$ and $v_x$ will be the same as the one derived in appendix A. The right eigenvectors read as

$$ 
\begin{aligned}
\mathbf{R}_1 &= \begin{bmatrix}
1 + \frac{c_s^2}{c^2} \\
v_x - c_s \sqrt{1 - \frac{v_x^2}{c^2}} \\
\left(1 + \frac{v_x c_s}{c^2 \sqrt{1 - \frac{v_x^2}{c^2}}}\right) v_y \\
\left(1 + \frac{v_x c_s}{c^2 \sqrt{1 - \frac{v_x^2}{c^2}}}\right) v_z
\end{bmatrix}, \quad
\mathbf{R}_4 &= \begin{bmatrix}
1 + \frac{c_s^2}{c^2} \\
v_x + c_s \sqrt{1 - \frac{v_x^2}{c^2}} \\
\left(1 - \frac{v_x c_s}{c^2 \sqrt{1 - \frac{v_x^2}{c^2}}}\right) v_y \\
\left(1 - \frac{v_x c_s}{c^2 \sqrt{1 - \frac{v_x^2}{c^2}}}\right) v_z
\end{bmatrix}, \\
\mathbf{R}_2 &= \begin{bmatrix}
0 \\
0 \\
1 \\
0
\end{bmatrix}, \quad
\mathbf{R}_3 &= \begin{bmatrix}
0 \\
0 \\
0 \\
1
\end{bmatrix},
\end{aligned}
\quad (3.20) $$

As can be seen from the eigenvectors $\mathbf{R}_2$ and $\mathbf{R}_3$ only the transverse velocity components change across the contact wave. The wave structure is depicted in figure 3. In appendix B we derive the expressions for the transverse velocities across different wave patterns.
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Structure of the solution of the 3D Riemann problem along the $x$ dimension. $\lambda_{1,4}$ are genuinely non-linear fields and can be either shock or rarefaction waves (thick lines). All state variables change across these waves. $\lambda_{2,3}$ are linearly degenerate fields and are associated to a contact discontinuity (thick dashed line). Only the transverse velocities $v_y$ and $v_z$ change across this wave.</figcaption></figure>

# 4 Riemann solvers

## 4.1 Exact Riemann solver

Having derived relations between primitive state variables $\mathbf{W}_L$, $\mathbf{W}_R$ and $\mathbf{W}^*$ across the different wave patterns, we can derive an exact solution of the RP for dark fluids as given by eq. (3.2). Since the solution for the transverse velocities is trivial we will focus here on the 1D case.

The first step consists in determining the value of $\Pi^*$ in the star region of figure 2. We have expressed the velocity $v^*$ in terms of $\Pi^*$ and the known states $\mathbf{W}_L$, $\mathbf{W}_R$ for all possible wave patterns. Thus, equating these relations for the left and right waves, $\Pi^*$ is the root of the algebraic equation

$$f_L(\Pi^*, \mathbf{W}_L) - f_R(\Pi^*, \mathbf{W}_R) = 0, \quad (4.1)$$

where

$$f_L(\Pi, \mathbf{W}_L) = \begin{cases} f_{Ls}(\Pi, \mathbf{W}_L) & \text{if } \Pi > \Pi_L \text{ (shock)} \\ f_{Lr}(\Pi, \mathbf{W}_L) & \text{if } \Pi \le \Pi_L \text{ (rarefaction)} \end{cases} \quad (4.2)$$

with $f_{Ls}(\Pi, \mathbf{W}_L)$ and $f_{Lr}(\Pi, \mathbf{W}_L)$ given by eq. (A.18) and eq. (A.7) respectively, while

$$f_R(\Pi, \mathbf{W}_R) = \begin{cases} f_{Rs}(\Pi, \mathbf{W}_R) & \text{if } \Pi > \Pi_R \text{ (shock)} \\ f_{Rr}(\Pi, \mathbf{W}_R) & \text{if } \Pi \le \Pi_R \text{ (rarefaction)} \end{cases} \quad (4.3)$$

with $f_{Rs}(\Pi, \mathbf{W}_R)$ and $f_{Rr}(\Pi, \mathbf{W}_R)$ given by eq. (A.22) and eq. (A.12) respectively. Eq. (4.1) can be solved for $\Pi^*$ to the desired level of accuracy using standard numerical root-finder schemes. Solving for $\Pi^*$ completely determines the nature of the waves. The value of $v^*$ can then be computed from one of the functions $f_L(\Pi^*, \mathbf{W}_L)$ or $f_R(\Pi^*, \mathbf{W}_R)$. Once the state $\mathbf{W}^*$ is known, the solution of the Riemann problem across the entire spatial interval can be sampled at any given time as following:
• in the region to the left of the discontinuity, corresponding to $x - x_0 < v^*\tau$:

– if $\Pi_L > \Pi^*$ the left wave is a rarefaction, then

$$ \mathbf{W}(x, \tau) = \begin{cases} (\Pi_L, v_L) & \text{if } x - x_0 < s_{\text{HLr}}\tau \\ (\Pi_{\text{Lfan}}, v_{\text{Lfan}}) & \text{if } s_{\text{HLr}}\tau < x - x_0 < s_{\text{TLr}}\tau \\ (\Pi^*, v^*) & \text{if } x - x_0 > s_{\text{TLr}}\tau \end{cases} \quad (4.4) $$

where $v_{\text{Lfan}}$ is given by solving eq. (A.8) and $\Pi_{\text{Lfan}}$ is given by eq. (A.9).

– if $\Pi_L < \Pi^*$ the left wave is a shock, then

$$ \mathbf{W}(x, \tau) = \begin{cases} (\Pi_L, v_L) & \text{if } x - x_0 < s_{\text{Ls}}\tau \\ (\Pi^*, v^*) & \text{if } x - x_0 > s_{\text{Ls}}\tau \end{cases} \quad (4.5) $$

where $s_{\text{Ls}}$ is given by eq. (A.17).

• in the region to the right of the discontinuity corresponding to $x - x_0 > v^*\tau$:

– if $\Pi_R > \Pi^*$ the right wave is a rarefaction, then

$$ \mathbf{W}(x, \tau) = \begin{cases} (\Pi^*, v^*) & \text{if } x - x_0 < s_{\text{TRr}}\tau \\ (\Pi_{\text{Rfan}}, v_{\text{Rfan}}) & \text{if } s_{\text{TRr}}\tau < x - x_0 < s_{\text{HRr}}\tau \\ (\Pi_R, v_R) & \text{if } x - x_0 > s_{\text{HRr}}\tau \end{cases} \quad (4.6) $$

where $v_{\text{Rfan}}$ is given by solving eq. (A.13) and $\Pi_{\text{Rfan}}$ is given by eq. (A.14).

– if $\Pi_R < \Pi^*$ the right wave is a shock, then

$$ \mathbf{W}(x, \tau) = \begin{cases} (\Pi^*, v^*) & \text{if } x - x_0 < s_{\text{Rs}}\tau \\ (\Pi_R, v_R) & \text{if } x - x_0 > s_{\text{Rs}}\tau \end{cases} \quad (4.7) $$

where $s_{\text{Rs}}$ is given by eq. (A.21).

We construct an exact Riemann solver as in [43]. We first solve eq. (4.1) using the Newton-Raphson method. The value of $\Pi^*$ is obtained to a given level of accuracy $\epsilon$ through an iteration procedure

$$ \Pi_{(i)}^* = \Pi_{(i-1)}^* - \frac{f[\Pi_{(i-1)}^*]}{f'[\Pi_{(i-1)}^*]}, \quad (4.8) $$

where $f(\Pi^*) = f_L(\Pi^*, \mathbf{W}_L) - f_R(\Pi^*, \mathbf{W}_R)$ and $f'(\Pi^*) = df/d\Pi|_*$ (the derivatives can be computed analytically). The iteration continues until the desired accuracy is reached

$$ \frac{|\Pi_{(i)}^* - \Pi_{(i-1)}^*|}{|\Pi_{(i)}^* + \Pi_{(i-1)}^*|} < \frac{\epsilon}{2}. \quad (4.9) $$

An initial guess value is necessary to start the iteration. Thus, the computation of the solution requires several iterations if the initial guess is too far off the solution. To address this point we use an adaptive scheme to optimize the initial guess using approximate Riemann solvers described in section 4.2.
<table><caption>**Table 1.** Initial data for three test problems with exact solution.</caption><thead><tr><th>Test</th><th>Π<sub>L</sub></th><th>v<sub>L</sub></th><th>Π<sub>R</sub></th><th>v<sub>R</sub></th></tr></thead><tbody><tr><td>1</td><td>1.0</td><td>0.0</td><td>0.125</td><td>0.0</td></tr><tr><td>2</td><td>1.0</td><td>-0.0003</td><td>1.0</td><td>0.0003</td></tr><tr><td>3</td><td>0.3</td><td>0.0001</td><td>0.3</td><td>-0.0001</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Exact solution of Test 1 for Π (left panel) and v (right panel) at time τ = 600 for $c_s = 10^{-3}$ (blue solid line) and $10^{-4}$ (red dashed line) respectively.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 5.** As in figure 4 for Test 2.</figcaption></figure>

In table 1 we quote the initial data for three standard test cases of the RP in arbitrary units. Notice that the equation of state of the fluid only enters via the source terms in the equations for Π and u (equation (2.11)), and thus does not appear in the RP test cases. Test 1 is the standard Sod test [45], whose solution consists of a left rarefaction and a right shock. Test 2 has a solution consisting of two rarefaction waves, while Test 3 consists of two shocks. We consider two different values of the speed of sound: $c_s = 10^{-3}$ and $10^{-4}$. We set the initial discontinuity at $x_0 = 1$ and sample the solution in the interval $0 < x < 2$ at time $τ = 600$. We set the accuracy of the Newton-Raphson method to $\epsilon = 10^{-6}$.

The exact solutions to Test 1, 2 and 3 are shown in figures 4, 5 and 6 respectively. In the case of Test 1 we can clearly see the opening of the rarefaction fan moving from right to
<figure><img src="image_9.png" /><figcaption>Figure 6. As in figure 4 for Test 3.</figcaption></figure>

left and the shock wave propagating to the right. For Test 2 we notice the opening of the two rarefaction fans propagating to the left and to the right respectively, while in the case of Test 3 we can see the density and velocity profile of the colliding shock waves. It is worth noting that these solutions only depend on the value of the sound speed, and are independent of the equation of state.

## 4.2 Approximate Riemann solvers

While the Riemann problem can be solved using an exact solver as the one described in section 4.1, for any practical application such as cosmological simulations this is computationally too expensive since it requires solving iteratively algebraic equations at each cell of the spatial domain of integration. For this reason, in this section, we present some approximate solvers. Let us notice here that in cosmological simulations one can combine different Riemann solvers in an adaptive fashion, where the more accurate and more expensive solvers are only used in case of large gradients, while the cheaper ones are used for smooth flow configurations, which represent the vast majority of the field.

### 4.2.1 Primitive variables Riemann solver

Let us consider the characteristic equations of the advection equation eq. (3.14) in primitive variables:

$$ \tilde{\mathbf{L}}_{\pm} \cdot d\mathbf{W} = 0, \qquad (4.10) $$

where $\tilde{\mathbf{L}}_{\pm}$ are the left eigenvectors eq. (3.18) and $d\mathbf{W} = (d\Pi, dv)$. These read:

$$ \frac{d\Pi}{\Pi} + \frac{1 + \frac{c_s^2}{c^2}}{c_s} \frac{dv}{\sqrt{1 - \frac{v^2}{c^2} + \frac{vc_s}{c^2}}} = 0 \quad \text{along } dx/d\tau = \lambda_+ \qquad (4.11) $$

$$ \frac{d\Pi}{\Pi} - \frac{1 + \frac{c_s^2}{c^2}}{c_s} \frac{dv}{\sqrt{1 - \frac{v^2}{c^2} - \frac{vc_s}{c^2}}} = 0 \quad \text{along } dx/d\tau = \lambda_- \qquad (4.12) $$

Following the derivation presented in [43], we can connect the star state to the left and right states by integrating these equations along the direction of the characteristics. We can
simplify the solution by defining the quantities:

$$C^+ = \frac{1 + \frac{c_s^2}{c^2}}{c_s} \frac{\Pi}{\sqrt{1 - \frac{v^2}{c^2} + \frac{vc_s}{c^2}}}, \quad (4.13)$$

$$C^- = \frac{1 + \frac{c_s^2}{c^2}}{c_s} \frac{\Pi}{\sqrt{1 - \frac{v^2}{c^2} - \frac{vc_s}{c^2}}}, \quad (4.14)$$

and approximating them as constant in time. This means that we can compute their value at the foot of the characteristic by only using the known left and right state variables. We can then connect the star state to the left (right) state by integrating eq. (4.11) (eq. (4.12)) along the characteristic of speed $\lambda_+$ ($\lambda_-$). This gives:

$$\Pi_L + C_{+,L}v_L = \Pi_* + C_{+,L}v_*, \quad (4.15)$$

$$\Pi_* - C_{-,R}v_* = \Pi_R - C_{-,R}v_R, \quad (4.16)$$

Solving the linear system of eqs. (4.15) and (4.16) we finally obtain the approximate state variables in the star region:

$$\Pi_* = \frac{C_{+,L}\Pi_R + C_{-,R}\Pi_L + C_{-,R}C_{+,L}(v_L - v_R)}{C_{-,R} + C_{+,L}} \quad (4.17)$$

$$v_* = \frac{\Pi_L - \Pi_R + C_{+,L}v_L + C_{-,R}v_R}{C_{-,R} + C_{+,L}} \quad (4.18)$$

We refer to these solutions as the Linear Primitive Variable Riemann Solver (LPVRS).

An alternative approximation can be found by defining:

$$C^+ = \frac{1 + \frac{c_s^2}{c^2}}{c_s} \frac{1}{\sqrt{1 - \frac{v^2}{c^2} + \frac{vc_s}{c^2}}}, \quad (4.19)$$

$$C^- = \frac{1 + \frac{c_s^2}{c^2}}{c_s} \frac{1}{\sqrt{1 - \frac{v^2}{c^2} - \frac{vc_s}{c^2}}}. \quad (4.20)$$

Integrating eq. (4.11) (eq. (4.12)) along the characteristic of speed $\lambda_+$ ($\lambda_-$) in this case gives:

$$\ln \Pi_L - \ln \Pi^* = -C_L^+ v_L + C_L^+ v^*, \quad (4.21)$$

$$\ln \Pi^* - \ln \Pi_R = C_R^- v^* - C_R^- v_R, \quad (4.22)$$

and the approximate state variables in the star region are:

$$v^* = \frac{C_L^+ v_L + C_R^- v_R + \ln \Pi_L - \ln \Pi_R}{C_L^+ + C_R^+}, \quad (4.23)$$

$$\Pi^* = \exp \left[ \frac{C_L^+ \ln \Pi_R + C_R^- \ln \Pi_L + C_L^+ C_R^- (v_L - v_R)}{C_L^+ + C_R^+} \right]. \quad (4.24)$$

We refer to this second approximation as the Primitive Variable Riemann Solver (PVRS). Notice that, unlike the LPVRS case, we no longer approximate $\Pi$ as being constant along the characteristic.
To find a simple solution for the state variables inside the rarefaction fan we make the approximation that $\sqrt{1 - \frac{v^2}{c^2}} \approx 1$ and expand the functions $g_{Lr}$ and $g_{Rr}$ given by equations (A.5) and (A.10) to first order. This gives:

$$v_{L}^{\text{fan}} = \frac{x - x_{0}}{\tau} + c_{s} \qquad (4.25)$$

$$\Pi_{L}^{\text{fan}} = \exp \left( \log \Pi_{L} + \frac{v_{L}}{c_{s}} + \frac{v_{L}c_{s}}{c^{2}} - \frac{v_{L}^{\text{fan}}}{c_{s}} - \frac{v_{L}^{\text{fan}}c_{s}}{c^{2}} \right) \qquad (4.26)$$

$$v_{R}^{\text{fan}} = \frac{x - x_{0}}{\tau} - c_{s} \qquad (4.27)$$

$$\Pi_{R}^{\text{fan}} = \exp \left( \log \Pi_{R} - \frac{v_{R}}{c_{s}} - \frac{v_{R}c_{s}}{c^{2}} + \frac{v_{R}^{\text{fan}}}{c_{s}} + \frac{v_{R}^{\text{fan}}c_{s}}{c^{2}} \right) \qquad (4.28)$$

We compare the results of these approximations to the exact solution for the tests in table 1 in figures 7-9. We can see that the PRVS performs very well in the case of rarefactions, where the solution is indistinguishable from the exact one, and only slightly misestimates $\Pi^*$ in the presence of shock waves. On the other hand, the LPVRS gives a larger error on the value of $\Pi^*$ in all the tests and on the value of $v^*$ in Test 1, especially for the larger value of $c_s$. While the PVRS approximation gives better results with respect to the LPVRS in all the tests, it is however susceptible to numerical instabilities, given its non-linear nature.

### 4.2.2 Harten, Lax and van Leer Riemann solver

The Harten, Lax and van Leer (HLL) Riemann solver approximates the solution as a system of two waves that separate three constant states. The intermediate state is found by averaging the exact solution between the fastest waves that travel to the left and right of the discontinuity [for a more detailed explanation see 43]. Given an estimate of these velocities $s_L$ and $s_R$, the HLL solver approximates the flux directly as:

$$\mathbf{F} = \begin{cases} \mathbf{F}_L & \text{if } x - x_0 \le s_L \tau \\ \mathbf{F}_{\text{HLL}} & \text{if } s_L \tau \le x - x_0 \le s_R \tau \\ \mathbf{F}_R & \text{if } x - x_0 \ge s_R \tau \end{cases} \qquad (4.29)$$

where

$$\mathbf{F}_{\text{HLL}} = \frac{s_R \mathbf{F}_L - s_L \mathbf{F}_R + s_L s_R (\mathbf{U}_R - \mathbf{U}_L)}{s_R - s_L}. \qquad (4.30)$$

The main advantage of this solver is its efficiency since it gives an estimate of the flux that involves only a few basic operations.

To test this solver we employ the simple estimates:

$$s_L = \min \left\{ v_L - c_s \sqrt{1 - \frac{v_L^2}{c^2}}, v_R - c_s \sqrt{1 - \frac{v_R^2}{c^2}} \right\}, \qquad (4.31)$$

$$s_R = \max \left\{ v_L + c_s \sqrt{1 - \frac{v_L^2}{c^2}}, v_R + c_s \sqrt{1 - \frac{v_R^2}{c^2}} \right\}. \qquad (4.32)$$

We verified that using more complex estimates for the speed does not significantly improve the accuracy of the solution, so that it is not worth spoiling the efficiency of this approximation. Since this solver gives an approximation for the flux, its performance is best shown using a Godunov solver. We compare the results of the tests in table 1 using the MUSCL-Hancock method described in section 5.1 with PVRS and HLL Riemann solvers in figures 10-12. We can see that the HLL solver performs almost identically to PVRS in these tests.
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Solution of Test 1 for Π (left panel) and υ (right panel) at time τ = 600 for $c_s = 10^{-3}$ (blue) and $10^{-4}$ (red) respectively for the exact solver (continuous line), LPVRS (dashed line) and PVRS (dotted line).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8.** As in figure 7 for Test 2.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 9.** As in figure 7 for Test 3.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Solution of Test 1 for Π (left panel) and $v$ (right panel) at time $\tau = 600$ for $c_s = 10^{-3}$ (blue) and $10^{-4}$ (red) for the exact Riemann solver (continuous line) and the MUSCL-Hancock method with PVRS (dotted line) and HLL Riemann solver (dashed line) respectively. Bottom panels: the absolute difference between the numerical and the exact solution with opposite signs for the two Riemann solvers to better see the distinction.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 11.** As in figure 10 for Test 2.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 12.** As in figure 10 for Test 3.</figcaption></figure>
## 5 Finite volume conservative methods

Let us consider the problem of solving eq. (3.1) numerically, provided with a set of initial and boundary conditions. The first step consists in discretizing the time interval into a finite number of time-steps $\tau^n$ of size $\Delta\tau$ and the spatial domain in a finite number of cells $I_i = [x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}]$ of size $\Delta x = L/M$, where $i = 1, \dots, M$ and $L$ is the size of the spatial domain of integration. This choice of spatial discretization is what characterizes finite volume methods, where the state $\mathbf{U}_i$ is regarded as a volume average of the variable $\mathbf{U}$ over the cell $I_i$. In the following we will denote the numerical approximation of $\mathbf{U}$ in the cell $I_i$ at time $\tau^n$ as $\mathbf{U}_i^n$.

A numerical solution of eq. (3.1) can be obtained by using the operator splitting approach, where the system is split into a purely advection problem with $\mathbf{U}_\tau + \mathbf{F}(\mathbf{U})_x = 0$ and an ordinary differential equation (ODE) $\mathbf{U}_\tau = \mathbf{S}(\mathbf{U})$. In such a case, one solves the PDE associated to the advection equations and then uses the solution to solve the ODE with a method that matches the order of accuracy of the solution of the advection problem such as to retain the same level of accuracy globally. In the following we focus on the solution of the advection equations, while the solution to the ODE will be implemented into the full cosmological simulation code that we will present in future work.

Following standard textbook presentations [see e.g. 43], the numerical solution to the advection problem can be obtained using a conservative scheme, such as that introduced by Godunov [37]:

$$ \mathbf{U}_i^{n+1} = \mathbf{U}_i^n + \frac{\Delta\tau}{\Delta x} [\mathbf{F}_{i-\frac{1}{2}} - \mathbf{F}_{i+\frac{1}{2}}], \quad (5.1) $$

where the $\mathbf{F}_{i\pm\frac{1}{2}}$ are inter-cell fluxes at the boundaries of the $i$-th cell given by

$$ \mathbf{F}_{i\pm\frac{1}{2}} = \mathbf{F}_{i\pm\frac{1}{2}}[\mathbf{U}_{i\pm\frac{1}{2}}(0)], \quad (5.2) $$

where $\mathbf{U}_{i+1/2}(0)$ is the solution of the Riemann problem $RP(\mathbf{U}_i^n, \mathbf{U}_{i+1}^n)$ at $x/t = 0$ and $\mathbf{U}_{i-1/2}(0)$ is the solution of the Riemann problem $RP(\mathbf{U}_{i-1}^n, \mathbf{U}_i^n)$ at $x/t = 0$. The size of the time-step $\Delta\tau$ satisfies the condition $\Delta\tau \le \Delta x/s_{\max}^n$, where $s_{\max}^n$ is the maximum wave velocity at time $\tau^n$ in the spatial domain of integration. This condition ensures that no wave in the solution of the Riemann problem travels more than the size of the cell $\Delta x$ in the time interval $\Delta\tau$. A time step satisfying the above condition can be set by introducing the Courant-Friedrichs-Lewy (CFL) coefficient $C_{cfl}$, such that

$$ \Delta\tau = C_{cfl} \frac{\Delta x}{s_{\max}^n}, \quad (5.3) $$

with $0 < C_{cfl} < 1$.

Notice that using piecewise constant data as initial condition for the Riemann problem, this method gives solutions that are only first order accurate in space and time. However, higher-order accuracy in space and/or time can be achieved by reconstructing the state variables using a truncated Taylor expansion and then use the reconstructed state as input for the Riemann solver. Here, we will derive working examples of some of these higher-order schemes such as the Monotonic Upstream-Centered Scheme for Conservation Laws (MUSCL-Hancock), the Piecewise Linear Method (PLM) and the Piecewise Parabolic Methods (PPM).

We test these schemes against the exact solutions of the Sod, the two rarefaction, and the two shock waves tests, which we have presented in section 4. Moreover, as we expect some of these schemes to improve the calculation of the flow in smooth regions rather than at discontinuities, we perform a pure advection test, which we discuss in section 5.4.
Hereafter, we use the Primitive Variables Riemann Solver (PVRS) described in section 4.2.1. We set $C_{cfl} = 0.9$ and estimate the maximum wave velocity as

$$s_{\max}^{n} = \max_{i}\{|v_{i}^{n}| + c_{s}\}, \qquad (5.4)$$

which provides a very good approximation in the case of rarefaction waves, though it may underestimate the speed of propagation in the case of shocks.

We assume periodic boundary conditions for the advection test:

$$\Pi_{1}^{n} = \Pi_{M}^{n}, \quad v_{1}^{n} = v_{M}^{n}, \qquad (5.5)$$

while for the other test we assume transmissive boundary conditions:

$$\Pi_{M+1}^{n} = \Pi_{M}^{n}, \quad v_{M+1}^{n} = v_{M}^{n} \qquad (5.6)$$

so that the boundaries do not affect the propagation of waves.

## 5.1 MUSCL-Hancock method

As already mentioned, Godunov's original scheme, which is based on a piecewise constant distribution of data, is first-order accurate in both space and time. In a series of seminal papers [46–48, 48–50] Van Leer presented a modification of Godunov's method to achieve higher-order accuracy. The idea is to use the initial piecewise constant data to reconstruct the distribution of state variables inside cells and extrapolate their values at the cell interfaces. Then, using the conservation equations these boundary values are evolved by half time-step to achieve second order accuracy in time. It is these reconstructed values that provide the initial data for the Riemann problem at the cell interface. While first-order schemes are guaranteed to preserve the monotonicity of the solution by Godunov's theorem, this is not the case for higher-order methods. In fact, the data reconstruction inside cells may generate local extrema that lead to spurious numerical oscillations in the solution. We thus need to impose monotonicity constraints on the reconstructed states, as we will see below.

The MUSCL-Hancock method introduced by Van Leer<fnref n="1" /> [51] is based on a linear reconstruction of the data distribution inside cells. The algorithm consists of the following step:

* Starting with the initial piecewise constant data $U_i^n$, the left and right cell boundary value are linearly extrapolated as

$$U_{i,L}^{n} = U_{i}^{n} - \frac{1}{2} \partial U_{i}, \quad U_{i,R}^{n} = U_{i}^{n} + \frac{1}{2} \partial U_{i}, \qquad (5.7)$$

where $\partial U_i$ denotes the slope of the linear interpolation in the i-th cell, whose determination will be described below.

* Evolve $U_{i,L}^n$ and $U_{i,R}^n$ by a time $\Delta\tau/2$ using the conservative update

$$\bar{U}_{i,L}^{n+1/2} = U_{i,L}^{n} + \frac{1}{2} \frac{\Delta\tau}{\Delta x} [\mathbf{F}(U_{i,L}^{n}) - \mathbf{F}(U_{i,R}^{n})], \qquad (5.8)$$

$$\bar{U}_{i,R}^{n+1/2} = U_{i,R}^{n} + \frac{1}{2} \frac{\Delta\tau}{\Delta x} [\mathbf{F}(U_{i,L}^{n}) - \mathbf{F}(U_{i,R}^{n})]. \qquad (5.9)$$

<footnote n="1">¹In [51] the idea behind this scheme is attributed to S. Hancock, see also [52], from which the acronym of MUSCL-Hancock method.</footnote>
<figure><img src="image_16.png" /><figcaption>Figure 13. MUSCL-Hancock (circles) vs. Exact solution (lines) of Test 1 for Π (left panel) and v (right panel) at time τ = 600 units for c_s = 10⁻³ (blue solid line) and 10⁻⁴ (red dashed line) respectively. In the bottom panel the absolute difference between the numerical and exact solution.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 14. As in figure 13 for Test 2.</figcaption></figure>

* Solve the Riemann Problem at inter-cell boundary with piecewise constant data $\bar{U}_{i,L}^{n+1/2}$ and $\bar{U}_{i,R}^{n+1/2}$:

$$ \mathbf{U}_{\tau} + \mathbf{F}(\mathbf{U})_x = 0, \qquad (5.10) $$

with initial conditions

$$ \mathbf{U}(x, 0) = \begin{cases} \bar{U}_{i,L}^{n+1/2} & \text{if } x < 0 \\ \bar{U}_{i+1,R}^{n+1/2} & \text{if } x > 0 \end{cases} \qquad (5.11) $$

A key point concerns the evaluation of the slopes $\partial \mathbf{U}_i$, which must satisfy monotonicity constraints. We implement the MUSCL-Hancock scheme with the MINMAX slope limiter which reads as:

$$ \partial \mathbf{U}_i = \frac{1}{2} [\mathrm{sign}(\mathbf{U}_i^n - \mathbf{U}_{i-1}^n) + \mathrm{sign}(\mathbf{U}_{i+1}^n - \mathbf{U}_i^n)] \times \min (|\mathbf{U}_i^n - \mathbf{U}_{i-1}^n|, |\mathbf{U}_{i+1}^n - \mathbf{U}_i^n|). \quad (5.12) $$

In figures 13, 14 and 15 we plot the numerical solutions (dotted lines) to test cases given in table 1 obtained with the MUSCL-Hancock method against the exact solutions at $\tau = 600$
<figure><img src="image_18.png" /><figcaption>Figure 15. As in figure 13 for Test 3.</figcaption></figure>

units for $c_s = 10^{-3}$ (solid blue lines) and $c_s = 10^{-4}$ (dash red lines) respectively. In the bottom panels we plot the absolute value of the difference between the numerical and exact solutions.

Overall the scheme performs rather well. A zoom on the value of $\Pi$ near the shocks in figures 13 and 15 reveals that shock waves are spread over ~3 cells, rather than having the zero-width of the exact solution. On the other hand, the speed is correctly estimated as can be seen on the right panels, thus indicating that the average shock position is well determined by the numerical scheme. Another characteristic feature of the solutions of Test 1 and 3 near the shocks is the absence of spurious oscillations. In the case of rarefaction waves, larger errors occur near the head and the tail as can be seen from figures 13 and 14, discontinuities are spread over ~10 cells. The numerical solution converges to the exact by increasing the spatial resolution of the grid.

## 5.2 Piecewise Linear Method (PLM)

This scheme, introduced by Colella [53], builds upon the MUSCL method and uses a linear reconstruction of the data to achieve second order accuracy in both space and time. We will follow the implementation presented in the lecture notes by M. Zingale [54] to which we refer the reader for a clear summary of higher-order reconstruction methods of non-linear advection equations.

The basic idea is to obtain left and right inter-cell states as a Taylor expansion of the primitive variables to first-order in $\Delta x/2$ to displace from the cell-centered value to the interface, and to first-order in $\Delta \tau/2$ to evolve the states to mid-point in time. A crucial aspect of this scheme is the so called *characteristic tracing step*, where only the waves that travel towards the interface are selected during the reconstruction process. This ensures that the correct flow of information is preserved.

These states are constructed as:

$$ \bar{\Pi}_{i,L}^{n+1/2} = \Pi_i^n + \frac{1}{2} [\tilde{\Pi}_i^{-} + \tilde{\Pi}_i^{+}] , \quad \bar{\Pi}_{i+1,R}^{n+1/2} = \Pi_{i+1}^n - \frac{1}{2} [\tilde{\Pi}_{i+1}^{-} + \tilde{\Pi}_{i+1}^{+}] , \quad (5.13) $$

with

$$ \tilde{\Pi}_i^{\pm} = \begin{cases} (1 - \frac{\Delta \tau}{\Delta x} \lambda_{\pm}(v_i)) (\tilde{L}_{1,i}^{\pm} \tilde{R}_{1,i}^{\pm} \partial \Pi_i + \tilde{L}_{2,i}^{\pm} \tilde{R}_{1,i}^{\pm} \partial v_i) & \text{if } \lambda_{\pm}(v_i) \ge 0 \\ 0 & \text{otherwise} \end{cases} \quad (5.14) $$
<figure><img src="image_19.png" /><figcaption>**Figure 16.** PLM (circles) vs. Exact solution (lines) of Test 1 for Π (left panel) and v (right panel) at time τ = 600 units for c_s = 10⁻³ (blue solid line) and 10⁻⁴ (red dashed line) respectively. In the bottom panel the absolute difference between the numerical and exact solution.</figcaption></figure>

and for the velocity:

$$ \tilde{v}_{i}^{\pm} = \begin{cases} \left(1 - \frac{\Delta\tau}{\Delta x}\lambda_{\pm}(v_i)\right) \left(\tilde{L}_{1,i}^{\pm} \tilde{R}_{2,i}^{\pm} \partial\Pi_i + \tilde{L}_{2,i}^{\pm} \tilde{R}_{2,i}^{\pm} \partial v_i\right) & \text{if } \lambda_{\pm}(v_i) \ge 0 \\ 0 & \text{otherwise} \end{cases} \quad (5.15) $$

where $\tilde{L}_1^\pm$ and $\tilde{L}_2^\pm$ are the elements of the left eigenvectors eq. (3.18) and $\tilde{R}_1^\pm$ and $\tilde{R}_2^\pm$ are the elements of the right eigenvectors eq. (3.17).

The characteristic tracing is evident because we put to 0 the contributions to $\tilde{\Pi}_i^\pm$ and $\tilde{v}_i^\pm$ when the wave $\lambda_\pm$ is moving away from the interface that we are considering.

In this scheme we also need to enforce monotonicity of the solution. To achieve this we compute the slope $\partial W_i = \{\partial \Pi_i, \partial v_i\}$ using the SUPERBEE slope limiter:

$$ \partial W_i = [\operatorname{sign}(\mathbf{W}_i^n - \mathbf{W}_{i-1}^n) + \operatorname{sign}(\mathbf{W}_{i+1}^n - \mathbf{W}_i^n)] \times \min \left[ |\mathbf{W}_i^n - \mathbf{W}_{i-1}^n|, |\mathbf{W}_{i+1}^n - \mathbf{W}_i^n|, \frac{1}{2}\max(|\mathbf{W}_i^n - \mathbf{W}_{i-1}^n|, |\mathbf{W}_{i+1}^n - \mathbf{W}_i^n|) \right]. \quad (5.16) $$

In figures 16, 17 and 18 we plot the PLM numerical solutions (dotted lines) to test cases given in table 1 against the exact solutions at $\tau = 600$ units for $c_s = 10^{-3}$ (solid blue lines) and $c_s = 10^{-4}$ (dash red lines) respectively. In the bottom panels we plot the absolute value of the difference between the numerical and exact solutions. We can see only minor differences with the respect to the results obtained with the MUSCL-Hancock method shown in figures 13, 14 and 15. In particular we may notice that the former is less accurate than MUSCL-Hancock in resolving rarefaction waves, while it performs better in the case of shocks.

## 5.3 Piecewise Parabolic Method (PPM)

Originally introduced by Colella & Woodward [55], this scheme achieves second-order accuracy using a parabolic reconstruction of the data with characteristic tracing. We refer again to Zingale's notes [54] for a clear summary of PPM in the case of non-linear advection equations. The expressions of the parabolic reconstruction of primitive variables with characteristics tracing step in the case of dark fluids are quite lengthy and here we only present the results of the numerical tests of table 1 shown in figures 19, 20 and 21 (dotted lines) against the
<figure><img src="image_20.png" /><figcaption>Figure 17. As in figure 16 for Test 2.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 18. As in figure 16 for Test 3.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Figure 19. PPM (circles) vs. Exact solution (lines) of Test 1 for Π (left panel) and v (right panel) at time τ = 600 units for $c_s = 10^{-3}$ (blue solid line) and $10^{-4}$ (red dashed line) respectively. In the bottom panel the absolute difference between the numerical and exact solution.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>Figure 20. As in figure 19 for Test 2.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Figure 21. As in figure 19 for Test 3.</figcaption></figure>

exact solutions at $\tau = 600$ units for $c_s = 10^{-3}$ (solid blue lines) and $c_s = 10^{-4}$ (dash red lines) respectively. In the bottom panels we plot the absolute value of difference between the exact and the numerical solutions.

## 5.4 Advection test for higher-order schemes

We consider a pure 1D advection problem and advect an initial Gaussian density profile

$$ \Pi(x, 0) \equiv \Pi_0(x) = 1 + e^{-\frac{(x-x_0)^2}{\sigma^2}}, \qquad (5.17) $$

located at $x_0 = 0.5$ with $\sigma = 0.1$ that moves with constant speed $v = 0.1$ toward the positive values of the x-axis. The exact solution to the advection equation is simply $\Pi(x, \tau) = \Pi_0(x - v \cdot \tau)$.

We solve numerically the 1D advection equation for the initial density profile $\Pi_0(x)$ in the case of $c_s = 10^{-3}$ with periodic boundary condition using the MUSCL-Hancock method, the PLM and PPM schemes previously discussed. The results are shown in figure 22, where we plot the exact solution at time $\tau = 5$ (blue solid line) and that obtained using MUSCL-Hancock (blue points), PLM (red points) and PPM (green points). The bottom panel shows
<figure><img src="image_25.png" /><figcaption>**Figure 22.** Advection of a Gaussian density profile (located at $x_0 = 0.5$ at the initial time) moving at constant speed $v = 0.1$ at time $\tau = 5$ units for $c_s = 10^{-3}$. Shown in the plot are the exact solution (blue solid line) and the numerical solutions obtained using the MUSCL-Hancock (blue points), PLM (red points) and PPM (green points) schemes respectively. In the bottom panel we plot the absolute difference between the exact and the numerical solutions.</figcaption></figure>

the absolute difference with respect to the exact solution. We can see the improvement in accuracy of the PPM scheme compared to the PLM and MUSCL-Hancock methods.

# 6 Conclusions

A dynamical form of dark energy is expected to lead to dark energy inhomogeneities whose properties are characterized by an effective speed of sound parameter. Cosmological observations of the cosmic microwave background as well as large scale structure have so far left the dark energy speed of sound unconstrained. However, ongoing and future galaxy surveys may have the potential to constrain this parameter, especially in the case of dark energy models with small sound speed values for which dark energy perturbations can become significant at small scales and late time. The detection of the effects of such inhomogeneities can provide smoking gun evidence against the cosmological constant paradigm, thus unveiling the dynamical nature of this exotic component. Nevertheless, this requires the availability of theoretical predictions capable of following the gravitational collapse of the coupled system of dark matter and dark energy throughout the non-linear regime of cosmic structure formation.

The gravitational collapse of dark matter density fluctuations can be followed using the standard N-body method. However, this is not the case for a dark energy fluid characterized by a negative equation of state and a non-vanishing sound speed. This is because in a given cosmological volume the mass of the particles tracing the clustering dark energy field is not conserved, rather it varies in space and time in a way that cannot be determined by the N-body equations of motion in a closed form. In contrast, it is possible to directly solve the Euler equations describing the evolution of the dark energy fluid using hydrodynamical
methods. Matter and dark energy perturbations are mutually coupled gravitationally via the Poisson equation for the gravitational potential.

Finite volume conservative methods have been developed in a vast literature to numerically solve hydrodynamic equations. These are implemented in established cosmological simulation codes that simulate the gravitational collapse of baryonic fluids. Hence, the use of such methods to solve the dark energy fluid equations can facilitate their implementation in simulation codes, thus allowing to perform full cosmological simulations of clustering dark energy models. One of the most widely adopted conservative numerical schemes is the Godunov method. This relies on the solution of the Riemann Problem of the advection part of the hydrodynamic equations considered. As such, the use of such a method to solve the dynamics of dark energy fluctuation in a cosmological simulation code requires the derivation of exact and approximate solutions that are specific to the form of the dark energy fluid equations.

In this work we have presented a reformulation of the Euler equations describing a dark energy fluid in a quasi-conservative form of hyperbolic partial differential equations with source terms. We have derived the exact solution of the Riemann Problem associated with the advection part of these equations in the non-relativistic limit with a small non-vanishing speed of sound. Physically, this corresponds to finding the wave structure of dark fluid perturbations propagating in a non-expanding background in the absence of gravitational interactions.

We have constructed a number of approximate solvers of the RP specific to the dark energy fluid equations and tested their accuracy against exact solutions of standard test cases, such as the Sod test, the propagation of two rarefaction waves and two shock waves. The use of such approximate solvers is necessary to significantly speed up the computation of the numerical solution of the Euler equations with the Godunov method. We have presented the equations describing the original formulation of the Godunov method that is first-order accurate in space and time. We have also derived the formula specific to the Euler equations of the dark fluid considered here for higher-order schemes based on the Godunov method, such as the MUSCL-Hancock method, the Piecewise Linear Method and the Piecewise Parabolic Method. We have tested the accuracy of these methods against exact solutions of standard test cases considered for the RP as well as an advection test.

In a follow-up work we will present the full cosmological simulation code, where we account for the cosmological source terms as well as the gravitational interaction and test the numerical solution in 3D.

To conclude, this work lays down the first stepping stone to constrain the speed of sound of Dark Energy perturbations using large-scale-structure data and opens the way to learning more about the Dark Energy phenomenon.

## Acknowledgments

We thank Eiichiro Komatsu, Emiliano Sefusatti, Volker Springel and Romain Teyssier for useful discussions, and the anonymous referee for helpful comments which lead to a clarified presentation. The research leading to these results has received funding from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013 Grant Agreement no. 279954). LB and FS acknowledge support from the Starting Grant (ERC-2015-STG 678652) "GrInflaGal" of the European Research Council.
# A Elementary wave solutions in 1D

## A.1 Rarefaction waves

We can use the Riemann invariants eq. (3.6) to relate a known state to the left (right) of a rarefaction wave to that on the right (left). These relations can be written in terms of primitive variables as:

$$ \frac{d\Pi}{\Pi} = \frac{\left(1 + \frac{c_s^2}{c^2}\right) dv}{c_s \sqrt{1 - \frac{v^2}{c^2} - \frac{c_s^2}{c^2} v}} \quad (\lambda = \lambda_+) \quad (A.1) $$

$$ \frac{d\Pi}{\Pi} = \frac{\left(1 + \frac{c_s^2}{c^2}\right) dv}{-c_s \sqrt{1 - \frac{v^2}{c^2} - \frac{c_s^2}{c^2} v}} \quad (\lambda = \lambda_-) \quad (A.2) $$

and can be integrated analytically to obtain:

$$ \ln \Pi - \frac{c}{c_s} \arcsin \left( \frac{v}{c} \right) + \ln \left( \sqrt{1 - \frac{v^2}{c^2} - \frac{vc_s}{c^2}} \right) = \text{const.} \quad (\lambda = \lambda_+) \quad (A.3) $$

$$ \ln \Pi + \frac{c}{c_s} \arcsin \left( \frac{v}{c} \right) + \ln \left( \sqrt{1 - \frac{v^2}{c^2} + \frac{vc_s}{c^2}} \right) = \text{const.} \quad (\lambda = \lambda_-) \quad (A.4) $$

which express the conservation of these quantities across the waves associated to the characteristic fields $\lambda_+$ and $\lambda_-$ respectively.

### A.1.1 Left rarefaction

Let us consider a left rarefaction wave associated with $\lambda_-$-field. Using the Riemann invariant eq. (A.4) we can relate the known state $U_L$ to the unknown state in the star region $U^*$ shown in figure 2. It is useful to introduce the auxiliary function

$$ g_{\text{Lr}}(v) \equiv \frac{c}{c_s} \arcsin \left( \frac{v}{c} \right) + \ln \left( \sqrt{1 - \frac{v^2}{c^2} + \frac{vc_s}{c^2}} \right), \quad (A.5) $$

hence from eq. (A.4) we obtain

$$ g_{\text{Lr}}(v^*) = \ln \frac{\Pi_L}{\Pi^*} + g_{\text{Lr}}(v_L). \quad (A.6) $$

This implicit relation can be inverted to obtain

$$ v^* = f_{\text{Lr}}(\Pi^*, \mathbf{W}_L), \quad (A.7) $$

where we have denoted with $f_{\text{Lr}}$ the inverse function $g_{\text{Lr}}^{-1}$.

The rarefaction wave is enclosed in a region delimited by the head of the wave moving with speed $s_{\text{HLr}} = \lambda_-(v_L)$ and the tail with speed $s_{\text{TLr}} = \lambda_+(v^*)$. The evolution of this region, also dubbed rarefaction fan, can be obtained by integrating the equation $dx/d\tau = \lambda_-(v)$
along the characteristic with initial condition $x(0) = x_0$ (the location of the discontinuity at initial time). This gives

$$ \frac{x - x_0}{\tau} = \lambda_-(v_{\text{Lfan}}), \qquad (A.8) $$

which can be inverted to obtain $v_{\text{Lfan}}(x)$ at any given time. $\Pi_{\text{Lfan}}(x)$ inside the fan is then obtained using the generalized Riemann invariant eq. (A.4):

$$ \ln \Pi_{\text{Lfan}}(x) = \ln \Pi_{\text{L}} + g_{\text{Lr}}(v_{\text{L}}) - g_{\text{Lr}}(v_{\text{Lfan}}(x)). \qquad (A.9) $$

### A.1.2 Right rarefaction

Let us now consider a right rarefaction wave associated with the characteristic field $\lambda_+$. Using the Riemann invariant eq. (A.3) we can relate the known state $U_R$ to the right of the wave to the unknown state in the star region $U^*$ shown in figure 2. Let us introduce the auxiliary function

$$ g_{\text{Rr}}(v) = -\frac{c}{c_s} \arcsin\left(\frac{v}{c}\right) + \ln\left(\sqrt{1 - \frac{v^2}{c^2}} - \frac{vc_s}{c^2}\right), \qquad (A.10) $$

from eq. (A.4) we obtain

$$ g_{\text{Rr}}(v^*) = \ln \frac{\Pi_{\text{R}}}{\Pi^*} + g_{\text{Rr}}(v_{\text{R}}). \qquad (A.11) $$

This implicit relation can be inverted to obtain

$$ v^* = f_{\text{Rr}}(\Pi^*, \mathbf{W}_R), \qquad (A.12) $$

where we have denoted with $f_{\text{Rr}}$ the inverse function $g_{\text{Rr}}^{-1}$.

The right rarefaction fan is delimited by the head of the wave moving with speed $s_{\text{HRr}} = \lambda_+(v_R)$ and the tail with speed $s_{\text{TRr}} = \lambda_+(v^*)$. The evolution of the fan region is given by solving $dx/d\tau = \lambda_+(v)$ with $x(0) = x_0$ (the location of the discontinuity at initial time):

$$ \frac{x - x_0}{\tau} = \lambda_+(v_{\text{Rfan}}), \qquad (A.13) $$

this can be inverted to infer $v_{\text{Rfan}}(x)$. Then, using the Riemann invariant eq. (A.3) we obtain $\Pi_{\text{Rfan}}(x)$ inside the fan:

$$ \ln \Pi_{\text{Rfan}}(x) = \ln \Pi_{\text{R}} + g_{\text{Rr}}(v_{\text{R}}) - g_{\text{Lr}}(v_{\text{Rfan}}(x)). \qquad (A.14) $$

## A.2 Shock waves

We can use the RH conditions in eq. (3.7) to relate a known state to the left (right) of the shock to that unknown to the right (left).

### A.2.1 Left shock

Let us consider a shock moving from right to left separating the known initial state $U_L$ from the unknown state $U^*$ in the star region represented in figure 2. Using eq. (3.7) we have a system of algebraic equations which reads as

$$ \left(1 + \frac{c_s^2}{c^2}\right) (\Pi_L v_L - \Pi^* v^*) = s_{\text{Ls}}(\Pi_L - \Pi^*) \qquad (A.15) $$

$$ \frac{c_s^2}{1 + \frac{c_s^2}{c^2}} (\Pi_L - \Pi^*) + (\Pi_L v_L^2 - \Pi^* v^2) = s_{\text{Ls}}(\Pi_L v_L - \Pi^* v^*). \qquad (A.16) $$
This is an algebraic system of two equations and three unknowns ($\Pi^*, v^*, s_{Ls}$). We can solve for $s_{Ls}$ and $v^*$ to find

$$ \frac{s_{Ls}}{1 + \frac{c_s^2}{c^2}} = \frac{\Pi_L v_L}{\Pi_L + \frac{c_s^2}{c^2} \Pi^*} - \sqrt{\frac{\Pi_L^2 v_L^2}{\left(\Pi_L + \frac{c_s^2}{c^2} \Pi^*\right)^2} - \frac{\Pi_L v_L^2}{\Pi_L + \frac{c_s^2}{c^2} \Pi^*} + \frac{\Pi^* c_s^2}{\left(\Pi_L + \frac{c_s^2}{c^2} \Pi^*\right) \left(1 + \frac{c_s^2}{c^2}\right)}} \quad (A.17) $$

$$ v^* \equiv f_{Ls}(\Pi^*, \mathbf{W}_L) = \frac{\Pi_L}{\Pi^*} v_L - \frac{s_{Ls}}{1 + \frac{c_s^2}{c^2}} \frac{\Pi_L - \Pi^*}{\Pi^*}, \quad (A.18) $$

where $s_{Ls}$ is the negative root of the quadratic equation obtained from eqs. (A.15) and (A.16), consistently with the fact that the left shock moves from the right to the left.

### A.2.2 Right shock

Let us consider a shock moving from the left to the right separating the known initial state $\mathbf{U}_R$ from the unknown state $\mathbf{U}^*$ in the star region represented in figure 2. Using the Rankine-Hugoniot conditions we obtain the system of algebraic equations

$$ \left(1 + \frac{c_s^2}{c^2}\right) (\Pi^* v^* - \Pi_R v_R) = s_{Rs}(\Pi^* - \Pi_R) \quad (A.19) $$

$$ \frac{c_s^2}{1 + \frac{c_s^2}{c^2}} (\Pi^* - \Pi_R) + (\Pi^* v^* - \Pi_R v_R^2) = s_{Rs}(\Pi^* v^* - \Pi_R v_R). \quad (A.20) $$

This is an algebraic system of two equations and three unknowns ($\Pi^*, v^*, s_{Rs}$), solving for $s_{Rs}$ and $v^*$ we find

$$ \frac{s_{Rs}}{1 + \frac{c_s^2}{c^2}} = \frac{\Pi_R v_R}{\Pi_R + \frac{c_s^2}{c^2} \Pi^*} + \sqrt{\frac{\Pi_R^2 v_R^2}{\left(\Pi_R + \frac{c_s^2}{c^2} \Pi^*\right)^2} - \frac{\Pi_R v_R^2}{\Pi_R + \frac{c_s^2}{c^2} \Pi^*} + \frac{\Pi^* c_s^2}{\left(\Pi_R + \frac{c_s^2}{c^2} \Pi^*\right) \left(1 + \frac{c_s^2}{c^2}\right)}} \quad (A.21) $$

$$ v^* \equiv f_{Rs}(\Pi^*, \mathbf{W}_R) = \frac{\Pi_R}{\Pi^*} v_R + \frac{s_{Rs}}{1 + \frac{c_s^2}{c^2}} \frac{\Pi^* - \Pi_R}{\Pi^*}, \quad (A.22) $$

where $s_{Rs}$ is the positive root of the quadratic equation obtained from eqs. (A.19) and (A.20), consistently with the fact that the right shock moves from the left to the right.

## B Elementary wave solutions in 3D

### B.1 Rarefaction waves

Defining:

$$ f_{1,4}(v_x) = 1 \pm \frac{v_x c_s}{c^2 \sqrt{1 - \frac{v_x^2}{c^2}}} \quad (B.1) $$

the Riemann Invariants for the $\lambda_{1,4}$ eigenvalue are given by:

$$ \frac{d\Pi}{1 + \frac{c_s^2}{c^2}} = \frac{d(\Pi v_x)}{\lambda_1} = \frac{d(\Pi v_y)}{v_y f_{1,4}(v_x)} = \frac{d(\Pi v_z)}{v_z f_{1,4}(v_x)}. \quad (B.2) $$
From the first equality we see that the Riemann invariant for $\Pi$ and $v_x$ remain the same as the 1D case. To find the values of $v_y$ across the wave we can solve:

$$ \frac{d\Pi}{1 + \frac{c_s^2}{c^2}} f_{1,4}(v_x) = \frac{d(\Pi v_y)}{v_y} \qquad (B.3) $$

By eliminating $\Pi$ using the first equality we can get to an equation that relates $v_y$ to $v_x$:

$$ \frac{-v_x \pm c_s \sqrt{1 - \frac{v_x^2}{c^2}}}{c^2 - v_x^2 \pm c_s v_x \sqrt{1 - \frac{v_x^2}{c^2}}} dv_x = \frac{dv_y}{v_y} \qquad (B.4) $$

which can be integrated to give:

$$ \frac{\frac{v_y}{c}}{\pm \frac{c_s v_x}{c^2} + \sqrt{1 - \frac{v_x^2}{c^2}}} = \text{const.} \qquad (B.5) $$

Similarly, for the third component of the velocity $v_z$:

$$ \frac{\frac{v_z}{c}}{\pm \frac{c_s v_x}{c^2} + \sqrt{1 - \frac{v_x^2}{c^2}}} = \text{const.} \qquad (B.6) $$

We can then find the values of the transverse velocities in the star region as:

$$ v_{y;L,R}^* = v_{y;L,R} \frac{\pm \frac{c_s u^*}{c^2} + \sqrt{1 - \frac{v_x^*}{c^2}}}{\pm \frac{c_s v_{x;L,R}}{c^2} + \sqrt{1 - \frac{v_{x;L,R}^2}{c^2}}} \qquad (B.7) $$

$$ v_{z;L,R}^* = v_{z;L,R} \frac{\pm \frac{c_s v_x^*}{c^2} + \sqrt{1 - \frac{v_x^*}{c^2}}}{\pm \frac{c_s v_{x;L,R}}{c^2} + \sqrt{1 - \frac{v_x^2}{c^2}}} \qquad (B.8) $$

Similarly, to find the solution inside the rarefaction fan we make use of the Riemann Invariants and compute:

$$ v_{y;L,R}^{\text{fan}} = v_{y;L,R} \frac{\pm \frac{c_s v_{x;L,R}^{\text{fan}}}{c^2} + \sqrt{1 - \frac{v_{x;L,R}^{\text{fan}}}{c^2}}}{\pm \frac{c_s v_{x;L,R}}{c^2} + \sqrt{1 - \frac{v_{x;L,R}^2}{c^2}}} \qquad (B.9) $$

$$ v_{z;L,R}^{\text{fan}} = v_{z;L,R} \frac{\pm \frac{c_s v_{x;L,R}^{\text{fan}}}{c^2} + \sqrt{1 - \frac{v_{x;L,R}^{\text{fan}}}{c^2}}}{\pm \frac{c_s v_{x;L,R}}{c^2} + \sqrt{1 - \frac{v_{x;L,R}^2}{c^2}}} \qquad (B.10) $$

## B.2 Shock waves

The Rankine-Hugoniot conditions in this case read:

$$ \left(1 + \frac{c_s^2}{c^2}\right) \left(\Pi_{L,R} v_{x;L,R} - \Pi^* v_x^*\right) = s_{L,R;s} \left(\Pi_{L,R} - \Pi^*\right) \qquad (B.11) $$

$$ \Pi_{L,R} v_x^2_{x;L,R} - \Pi^* v_x^{*2} \pm \frac{c_s^2}{1 + \frac{c_s^2}{c^2}} \left(\Pi_{L,R} - \Pi^*\right) = s_{L,R;s} \left(\Pi_{L,R} v_x_{x;L,R} - \Pi^* v_x^*\right) \qquad (B.12) $$

$$ \Pi_{L,R} v_x_{x;L,R} v_y_{y;L,R} - \Pi^* v_x^* v_y^* = s_{L,R;s} \left(\Pi_{L,R} v_y_{y;L,R} - \Pi^* v_y^*\right) \qquad (B.13) $$

$$ \Pi_{L,R} v_x_{x;L,R} v_z_{z;L,R} - \Pi^* v_x^* v_z^* = s_{L,R;s} \left(\Pi_{L,R} v_z_{z;L,R} - \Pi^* v_z^*\right) \qquad (B.14) $$
Again we see that the solution for $\Pi$ and $v_x$ remains identical to the 1D case, while we can easily solve equations (B.13) and (B.14) to find:

$$v_{y;L,R}^* = \Pi_{L,R} v_{y;L,R} \frac{s_{L,R;s} - v_{x;L,R}}{\Pi^*(s_{L,R;s} - v_x^*)} \quad (B.15)$$

$$v_{z;L,R}^* = \Pi_{L,R} v_{z;L,R} \frac{s_{L,R;s} - v_{x;L,R}}{\Pi^*(s_{L,R;s} - v_x^*)} \quad (B.16)$$

## References

[1] PLANCK collaboration, Planck 2018 results. VI. Cosmological parameters, Astron. Astrophys. 641 (2020) A6 [Erratum ibid. 652 (2021) C4] [arXiv:1807.06209] [INSPIRE].

[2] BOSS collaboration, The clustering of galaxies in the completed SDSS-III Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample, Mon. Not. Roy. Astron. Soc. 470 (2017) 2617 [arXiv:1607.03155] [INSPIRE].

[3] DES collaboration, Dark Energy Survey Year 3 results: Cosmological constraints from galaxy clustering and weak lensing, Phys. Rev. D 105 (2022) 023520 [arXiv:2105.13549] [INSPIRE].

[4] T. Hamana et al., Cosmological constraints from cosmic shear two-point correlation functions with HSC survey first-year data, Publ. Astron. Soc. Jap. 72 (2020) 16 [arXiv:1906.06041] [INSPIRE].

[5] S.M. Carroll, The Cosmological constant, Living Rev. Rel. 4 (2001) 1 [astro-ph/0004075] [INSPIRE].

[6] L. Verde, T. Treu and A.G. Riess, Tensions between the Early and the Late Universe, Nature Astron. 3 (2019) 891 [arXiv:1907.10625] [INSPIRE].

[7] C. Heymans et al., KiDS-1000 Cosmology: Multi-probe weak gravitational lensing and spectroscopic galaxy clustering constraints, Astron. Astrophys. 646 (2021) A140 [arXiv:2007.15632] [INSPIRE].

[8] S. Tsujikawa, Modified gravity models of dark energy, Lect. Notes Phys. 800 (2010) 99 [arXiv:1101.0191] [INSPIRE].

[9] T. Clifton, P.G. Ferreira, A. Padilla and C. Skordis, Modified Gravity and Cosmology, Phys. Rept. 513 (2012) 1 [arXiv:1106.2476] [INSPIRE].

[10] E.J. Copeland, M. Sami and S. Tsujikawa, Dynamics of dark energy, Int. J. Mod. Phys. D 15 (2006) 1753 [hep-th/0603057] [INSPIRE].

[11] P. Creminelli, G. D'Amico, J. Noreña and F. Vernizzi, The Effective Theory of Quintessence: the $w < -1$ Side Unveiled, JCAP 02 (2009) 018 [arXiv:0811.0827] [INSPIRE].

[12] J. Gleyzes, D. Langlois and F. Vernizzi, A unifying description of dark energy, Int. J. Mod. Phys. D 23 (2015) 1443010 [arXiv:1411.3712] [INSPIRE].

[13] R.R. Caldwell, R. Dave and P.J. Steinhardt, Cosmological imprint of an energy component with general equation of state, Phys. Rev. Lett. 80 (1998) 1582 [astro-ph/9708069] [INSPIRE].

[14] C. Armendariz-Picon, V.F. Mukhanov and P.J. Steinhardt, Essentials of k essence, Phys. Rev. D 63 (2001) 103510 [astro-ph/0006373] [INSPIRE].

[15] D. Bertacca, S. Matarrese and M. Pietroni, Unified Dark Matter in Scalar Field Cosmologies, Mod. Phys. Lett. A 22 (2007) 2893 [astro-ph/0703259] [INSPIRE].

[16] J. Khoury, Alternative to particle dark matter, Phys. Rev. D 91 (2015) 024022 [arXiv:1409.0012] [INSPIRE].
[17] J.K. Erickson et al., *Measuring the speed of sound of quintessence*, Phys. Rev. Lett. **88** (2002) 121301 [astro-ph/0112438] [INSPIRE].

[18] R. Bean and O. Doré, *Probing dark energy perturbations: The Dark energy equation of state and speed of sound as measured by WMAP*, Phys. Rev. D **69** (2004) 083503 [astro-ph/0307100] [INSPIRE].

[19] P.-S. Corasaniti, T. Giannantonio and A. Melchiorri, *Constraining dark energy with cross-correlated CMB and large scale structure data*, Phys. Rev. D **71** (2005) 123521 [astro-ph/0504115] [INSPIRE].

[20] R. de Putter, D. Huterer and E.V. Linder, *Measuring the Speed of Dark: Detecting Dark Energy Perturbations*, Phys. Rev. D **81** (2010) 103513 [arXiv:1002.1311] [INSPIRE].

[21] J.-M. Alimi et al., *Imprints of dark energy on cosmic structure formation — I. Realistic quintessence models and the non-linear matter power spectrum*, Mon. Not. Roy. Astron. Soc. **401** (2010) 775 [arXiv:0903.5490] [INSPIRE].

[22] E. Jennings, C.M. Baugh, R.E. Angulo and S. Pascoli, *Simulations of Quintessential Cold Dark Matter: beyond the cosmological constant*, Mon. Not. Roy. Astron. Soc. **401** (2010) 2181 [arXiv:0908.1394] [INSPIRE].

[23] J. Courtin et al., *Imprints of dark energy on cosmic structure formation: II) Non-Universality of the halo mass function*, Mon. Not. Roy. Astron. Soc. **410** (2011) 1911 [arXiv:1001.3425] [INSPIRE].

[24] D.F. Mota and C. van de Bruck, *On the Spherical collapse model in dark energy cosmologies*, Astron. Astrophys. **421** (2004) 71 [astro-ph/0401504] [INSPIRE].

[25] P. Creminelli et al., *Spherical collapse in quintessence models with zero speed of sound*, JCAP **03** (2010) 027 [arXiv:0911.2701] [INSPIRE].

[26] T. Basse, O.E. Bjaelde and Y.Y.Y. Wong, *Spherical collapse of dark energy with an arbitrary sound speed*, JCAP **10** (2011) 038 [arXiv:1009.0010] [INSPIRE].

[27] E. Sefusatti and F. Vernizzi, *Cosmological structure formation with clustering quintessence*, JCAP **03** (2011) 047 [arXiv:1101.1026] [INSPIRE].

[28] S. Anselmi, G. Ballesteros and M. Pietroni, *Non-linear dark energy clustering*, JCAP **11** (2011) 014 [arXiv:1106.0834] [INSPIRE].

[29] S. Anselmi, D. López Nacir and E. Sefusatti, *Nonlinear effects of dark energy clustering beyond the acoustic scales*, JCAP **07** (2014) 013 [arXiv:1402.4269] [INSPIRE].

[30] F. Hassani, J. Adamek, M. Kunz and F. Vernizzi, *k-evolution: a relativistic N-body code for clustering dark energy*, JCAP **12** (2019) 011 [arXiv:1910.01104] [INSPIRE].

[31] F. Hassani, J. Adamek and M. Kunz, *Clustering dark energy imprints on cosmological observables of the gravitational field*, Mon. Not. Roy. Astron. Soc. **500** (2020) 4514 [arXiv:2007.04968] [INSPIRE].

[32] F. Hassani et al., *New nonlinear instability for scalar fields*, Phys. Rev. D **105** (2022) L021304 [arXiv:2107.14215] [INSPIRE].

[33] F. Hassani et al., *A new instability in clustering dark energy?*, arXiv:2204.13098 [INSPIRE].

[34] G.L. Bryan and M.L. Norman, *Simulating x-ray clusters with adaptive mesh refinement*, ASP Conf. Ser. **123** (1997) 363 [astro-ph/9710186] [INSPIRE].

[35] B. Fryxell et al., *FLASH: An Adaptive Mesh Hydrodynamics Code for Modeling Astrophysical Thermonuclear Flashes*, Astrophys. J. Suppl. **131** (2000) 273 [INSPIRE].

[36] R. Teyssier, *Cosmological hydrodynamics with adaptive mesh refinement: a new high resolution code called ramses*, Astron. Astrophys. **385** (2002) 337 [astro-ph/0111367] [INSPIRE].
[37] S. Godunov, *A finite difference method for the computation of discontinuous solutions of the equations of fluid dynamics*, Sbornik Math. **47** (1959) 357.

[38] P. Colella and H.M. Glaz, *Efficient solution algorithms for the Riemann problem for real gases*, J. Comput. Phys. **59** (1985) 264.

[39] J.J. Gottlieb and C.P.T. Groth, *Assessment of Riemann solvers for unsteady one-dimensional inviscid flows of perfect gases*, J. Comput. Phys. **78** (1988) 437.

[40] R. Menikoff and B.J. Plohr, *The Riemann problem for fluid flow of real materials*, Rev. Mod. Phys. **61** (1989) 75 [INSPIRE].

[41] M.J. Ivings, D.M. Causon and E.F. Toro, *On Riemann solvers for compressible liquids*, Int. J. Numer. Meth. Fluids **28** (1998) 395.

[42] R.J. LeVeque, *Finite volume methods for hyperbolic problems*, Cambridge University Press (2004) [DOI:10.1017/CBO9780511791253].

[43] E.F. Toro, *Riemann solvers and numerical methods for fluid dynamics: a practical introduction*, Springer-Verlag (2009) [DOI:10.1007/b79761].

[44] A. Jeffrey, *Quasilinear hyperbolic systems and waves*, Pitman Publishing (1976).

[45] G.A. Sod, *A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws*, J. Comput. Phys. **27** (1978) 1.

[46] B. van Leer, *Towards the ultimate conservative difference scheme I. The quest of monotonicity*, in *Proceedings of the Third International Conference on Numerical Methods in Fluid Mechanics*, H. Cabannes and R. Temam eds., Lect. Notes Phys. **18** (1973) 163.

[47] B. van Leer, *Towards the ultimate conservative difference scheme. II. Monotonicity and conservation combined in a second-order scheme*, J. Comput. Phys. **14** (1974) 361.

[48] B.V. Leer, *Towards the ultimate conservative difference scheme III. Upstream-centered finite-difference schemes for ideal compressible flow*, J. Comput. Phys. **23** (1977) 263.

[49] B.V. Leer, *Towards the ultimate conservative difference scheme. IV. A new approach to numerical convection*, J. Comput. Phys. **23** (1977) 276.

[50] B. van Leer, *Towards the ultimate conservative difference scheme. V. A second-order sequel to Godunov's method*, J. Comput. Phys. **32** (1979) 101.

[51] B. Van Leer, *On the relation between the upwind-differencing schemes of Godunov, Engquist-Osher and Roe*, SIAM J. Sci. Stat. Comput. **5** (1984) 1.

[52] B. Van Leer, *Upwind and high-resolution methods for compressible flow: From donor cell to residual-distribution schemes*, in the proceedings of the 16th AIAA Computational Fluid Dynamics Conference, (2006), p. 3559 [DOI:10.2514/6.2003-3559].

[53] P. Colella, *A direct Eulerian MUSCL scheme for gas dynamics*, SIAM J. Sci. Stat. Comput. **6** (1985) 104.

[54] M. Zingale, *Introduction to Computational Astrophysical Hydrodynamics*, http://bender.astro.sunysb.edu/hydro_by_example/CompHydroTutorial.pdf.

[55] P. Colella and P.R. Woodward, *The Piecewise Parabolic Method (PPM) for Gas Dynamical Simulations*, J. Comput. Phys. **54** (1984) 174 [INSPIRE].