# Nonideal self-gravity and cosmology: Importance of correlations in the dynamics of the large-scale structures of the Universe

P. Tremblin, G. Chabrier, T. Padioleau, S. Daley-Yates

► To cite this version:

P. Tremblin, G. Chabrier, T. Padioleau, S. Daley-Yates. Nonideal self-gravity and cosmology: Importance of correlations in the dynamics of the large-scale structures of the Universe. Astronomy and Astrophysics - A&A, 2022, 659, pp.A108. 10.1051/0004-6361/202142103. cea-03609788

HAL Id: cea-03609788

https://cea.hal.science/cea-03609788

Submitted on 15 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# Nonideal self-gravity and cosmology: Importance of correlations in the dynamics of the large-scale structures of the Universe

P. Tremblin<fnref n="1" />, G. Chabrier<fnref n="2" />,<fnref n="3" />, T. Padioleau<fnref n="1" />, and S. Daley-Yates<fnref n="1" />

Received 28 August 2021 / Accepted 9 December 2021

## ABSTRACT

**Aims.** Inspired by the statistical mechanics of an ensemble of interacting particles (BBGKY hierarchy), we propose to account for small-scale inhomogeneities in self-gravitating astrophysical fluids by deriving a nonideal virial theorem and nonideal Navier-Stokes equations. These equations involve the pair radial distribution function (similar to the two-point correlation function used to characterize the large-scale structures of the Universe), similarly to the interaction energy and equation of state in liquids. Within this framework, small-scale correlations lead to a nonideal amplification of the gravitational interaction energy, whose omission leads to a missing mass problem, for instance, in galaxies and galaxy clusters.

**Methods.** We propose to use a decomposition of the gravitational potential into a near- and far-field component in order to account for the gravitational force and correlations in the thermodynamics properties of the fluid. Based on the nonideal virial theorem, we also propose an extension of the Friedmann equations in the nonideal regime and use numerical simulations to constrain the contribution of these correlations to the expansion and acceleration of the Universe.

**Results.** We estimate that the nonideal amplification factor of the gravitational interaction energy of the baryons lies between 5 and 20, potentially explaining the observed value of the Hubble parameter (since the uncorrelated energy accounts for ~5%). Within this framework, the acceleration of the expansion emerges naturally because the number of substructures induced by gravitational collapse increases, which in turn increases their contribution to the total gravitational energy. A simple estimate predicts a nonideal deceleration parameter $q_{ni} \simeq -1$; this is potentially the first determination of the observed value based on an intuitively physical argument. We also suggest that small-scale gravitational interactions in bound structures (spiral arms or local clustering) could yield a transition to a viscous regime that can lead to flat rotation curves. This transition could also explain the dichotomy between (Keplerian) low surface brightness elliptical galaxy and (nonkeplerian) spiral galaxy rotation profiles. Overall, our results demonstrate that nonideal effects induced by inhomogeneities must be taken into account, potentially with our formalism, in order to properly determine the gravitational dynamics of galaxies and the large-scale Universe.

**Key words.** equation of state – gravitation – cosmology: theory

## 1. Introduction

Astrophysical flows are by nature multiscale systems exhibiting a wide range of dynamical regimes that is remains quite challenging to understand. Since the 1970s and 1980s, hydro- or N-body high- performance numerical simulations have been key tools for understanding self-gravitating systems. However, a proper understanding of collective effects in these simulations remains challenging, notably because of numerical artifacts. Therefore, theory remains vital to assess the genuine validity of the numerical simulations.

The theoretical description of an ensemble of particles is a challenging task that has motivated the development of a major domain in physics: statistical mechanics. The statistical descriptions of a gas and a solid are relatively easy in two limiting cases: perfect disorder, often referred to as molecular chaos for ideal gases, or perfect order for ordered solids. The description of liquids, or ill-condensed systems, however, is much more complex and has only been possible after the pioneering work of Bogolyubov, Born, Green, Kirkwood, and Yvon, referred to as the BBGKY hierarchy (e.g., Yvon 1935; Born & Green 1946; Bogoliubov 1946; Kirkwood 1946). Within this framework, the

ensemble of particles is described by the one-particle probability density function, the two-particle probability density function, and so on. Truncated at the second order, the system is described by the one-particle density function and the pair correlation function, generally referred to as the radial distribution function in liquid statistical mechanics. In the case of perfect chaos for an ideal gas, the radial distribution function is uniformly equal to 1, which means that there is no correlation in the fluid: the probability of finding two particles at two given points is equal to the product of the probabilities of finding one particle at these given points.

Liquids differ from an ideal gas because of the presence of interparticle forces that lead to correlations at short distances, and they differ from solids because of the lack of (periodic) long-range order (i.e., the pair correlation function tends to 1 at long distance). For water, the dominant effects responsible for these intermolecular forces are dipolar interactions and hydrogen bonding. Many other processes can be at play and make the picture more complex. However, from a statistical point of view, the knowledge of the radial distribution function is sufficient to derive the thermodynamic properties of the fluid, such as the equation of state (see Hansen et al. 2006; Aslangul 2006).

<footnote n="1">¹ Maison de la Simulation, CEA, CNRS, Univ. Paris-Sud, UVSQ, Université Paris-Saclay, 91191 Gif-sur-Yvette, France
e-mail: pascal.tremblin@cea.fr</footnote>

<footnote n="2">² Ecole Normale Supérieure de Lyon, CRAL, UMR CNRS 5574, 69364 Lyon Cedex 07, France</footnote>

<footnote n="3">³ Astrophysics Group, University of Exeter, EX4 4QL Exeter, UK</footnote>
Most of the developments of liquid statistical mechanics deeply rely on the fact that the interactions at play are short range: even though the interaction energy is proportional to the number of particle pairs in the fluid, only the close neighbors actually matter. This allows determining an extensive definition of the interaction energy that does not diverge in the thermodynamic limit. A noticeable exception is the case of Coulombic fluids, which interact throughout the Coulomb potential ($\propto 1/r$); large-scale convergence in this case is ensured by the neutralizing electron background. Unfortunately, this (screening) property is not met by a self-gravitating fluid because the gravitational force is long-range and always attractive: the interaction energy is a priori nonextensive and might require a new extension of statistical mechanics to include nonextensive hierarchical systems (see, e.g., de Vega & Sánchez 2002; Pfenniger 2006). This prevents astrophysicists from benefiting from the insight found through developments of statistical mechanics in the presence of interactions. For this reason, at large scales such as galaxies or the large-scale structures of the Universe, the gravitational interactions of unresolved structures in the simulations are usually ignored or neglected because the thermodynamic properties of the fluid are calculated assuming ideal gas approximation. This approximation, however, can lead to inconsistencies in our description of the Universe.

The present paper proposes properly applying statistical mechanics to gravitational systems by solving the problem of nonextensivity. This leads to the derivation of a nonideal virial theorem and nonideal Navier-Stokes equations to be used for self-gravitating hydrodynamic fluids when, as is the case within the Universe, correlations due to self-gravitating substructures are present. In Sect. 2 we recall classical results of statistical mechanics in order to define a nonideal equation of state (EOS) in the presence of interactions. In Sect. 3 we propose a way to adapt these tools to the long-range gravitational force by using a near-field and far-field decomposition of the potential. In Sect. 4 we present a semianalytical example using polytropic stellar structures to illustrate the importance of inhomogeneities and correlations when computing the correct interaction energy in the virial theorem. Then, we explore the role of the gravitational force and the correlations at small scales in galaxies to explain the observed flat rotation curves. Finally, using the Newtonian limit for now, we define nonideal Friedmann equations and explore the possible consequences of these nonideal effects in the context of the expanding and accelerating Universe. In Sect. 5 we summarize our conclusions, discuss the limitations of our work, and propose ways to proceed in exploring this potentially interesting idea.

## 2. Statistical mechanics of an ensemble of interacting particles

### 2.1. BBGKY hierarchy

Following the BBGKY hierarchy, we introduce some notations of statistical mechanics. We consider N particles of mass m within a volume V. We can define $P_s(r_1, r_2, ..., r_s)$, the probability of finding s particles at the positions $(r_1, r_2, ..., r_s)$. These probability density functions are normalized by the total number of N-uplets $N^s$. Two probability functions play an important role when there is no interaction involving more than two particles: $P_1(r)$, the probability of finding one particle at r and $P_2(r, r')$, and the probability of finding one particle at r and one particle at r'. These two probability functions can be used to exactly

define fluid quantities involving one and two particles from their N-body counterparts.

For instance, using $P_1(r)$, we can equivalently define the total acceleration $\langle a \rangle$ of the fluid volume V and the sum of all the N-body particle accelerations,

$$ \langle a \rangle = \sum_i a_i \\ = N \int_V a(r) P_1(r) dV. \qquad (1) $$

The use of $P_1(r)$ here involves counting the number of particles in the volume V since by definition, the density inside the fluid volume is

$$ \rho(r) = N P_1(r). \qquad (2) $$

In the case of a homogeneous and isotropic fluid, $P_1(r)$ is simply equal to $1/V$ for all r, and $\rho$ is equal to $N/V$, the number density.

Similarly, the interaction energy of the ensemble of particles $\langle H_{\text{int}} \rangle$ can be defined by summing all the interaction energies between distinct pairs of particles in the N-body description or by integrating over $P_2(r, r')$ for the fluid volume,

$$ \langle H_{\text{int}} \rangle = \sum_{i<j} \phi(|r_i - r_j|) \\ = \frac{N(N-1)}{2} \iint_{V,V} \phi(|r - r'|) P_2(r, r') dV dV', \qquad (3) $$

where $N(N-1)/2$ is the number of independent pairs in the fluid. In the homogeneous and isotropic case, $P_2(r, r')$ only depends on the distance modulus $||r - r'||$. It is then convenient to define the radial distribution function

$$ g(\|r - r'\|) = V^2 P_2(r, r'), \qquad (4) $$

or equivalently, the correlation function, $\xi(\|r - r'\|)$, for a homogeneous isotropic fluid,

$$ \begin{aligned} \xi(\|r - r'\|) &= V^2(P_2(r, r') - P_1(r)P_1(r')) \\ &= V^2 P_2(r, r') - 1 \\ &= g(\|r - r'\|) - 1. \end{aligned} \qquad (5) $$

For an isotropic fluid, the correlation function $\xi(r)$ is the angle-averaged value of $\xi(r) = \langle \delta(r')\delta(r + r') \rangle$, with $\delta(r)$ the density perturbation defined as $\delta(r) = (\rho(r) - \langle \rho \rangle) / \langle \rho \rangle$ (see Springel et al. 2017). Here, $\langle \cdot \rangle$ denotes a volume average. The power spectrum can then be defined as the Fourier transform of $\xi(r)$,

$$ \xi(r) = \frac{V}{(2\pi)^3} \int |\delta_k|^2 e^{-ik \cdot r} d^3 k, \qquad (6) $$

which is commonly used in cosmological studies to characterize the large-scale structures of the Universe (Peebles 1980; Peacock 1999). We show here that the correlation function is not a simple diagnostic of structures, but plays an active role in determining the dynamics of the fluid. For the sake of similarity with the tools that are commonly used in statistical physics, we use the radial distribution function rather than the correlation function in our calculations.

It is easy to see that in the absence of correlation in the fluid, $P_2(r, r') = P_1(r)P_1(r')$, that is, $\xi(\|r - r'\|) = 0$ and $g(\|r - r'\|) = 1$ for all r and r' in the case of a homogeneous fluid. This situation is the limit case of an ideal fluid, often referred to as “perfect
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Schematic illustration of situations in which the ideal hypothesis is valid and in which it breaks down. Left: fluid is uncorrelated at all scales. Right: on the scale of control volume $V_1$, the fluid can be considered locally as inhomogeneous and uncorrelated, while on the scale of control volume $V_2$, the fluid can be considered as globally homogeneous and correlated.</figcaption></figure>

molecular chaos”, for which the radial distribution function is equal to 1 everywhere. A schematic representation of such an ideal and homogeneous fluid is given in the left panel of Fig. 1, where the blue dots represent point-like interactionless particles. In this example, the fluid is homogeneous at all scales: the density inside a small ($V_1$) or large ($V_2$) control volume is independent of the location of the control volume.

In the presence of interactions, the density distribution may no longer be homogeneous at all scales, for instance, in the presence of clustering in the particle distribution. This situation is shown in the right panel of Fig. 1. For a small control volume $V_1$, the distribution of particles within the volume is ideal and there is no correlation: $P_2(\mathbf{r}, \mathbf{r}') \neq P_1(\mathbf{r})P_1(\mathbf{r}')$. However, the density depends on the location of the control volume, for example, whether it lies within a cluster or outside. This means that inhomogeneities are resolved and the fluid is ideal (i.e., uncorrelated), but inhomogeneous. For a large control volume $V_2$, the distribution of particles within the volume is not ideal and correlations are induced by the clustering inside the volume: $P_2(\mathbf{r}, \mathbf{r}') \neq P_1(\mathbf{r})P_1(\mathbf{r}')$. However, the density does not depend on the location of the control volume provided that the control volume is sufficiently large. In that case, the inhomogeneities are not resolved and the fluid is nonideal (correlated) and homogeneous.

This distinction is very important for the correct treatment of inhomogeneities (correlations) in a fluid. Because it depends on the scale, this is of prime importance for astrophysical flows, which can only be observed or simulated at large scales. The main consequence is that we cannot simply average over large volumes as this masks the complexity of the density inhomogeneities: these inhomogeneities then result in correlations in the fluid, which requires a change in dynamical equations. This can be seen with the computation of the interaction energy, which intrinsically depends on the distance between particles (see Eq. (3)). Within a statistical (fluid) description, the interaction energy fundamentally involves the pair correlation function, $P_2(\mathbf{r}, \mathbf{r}') \neq P_1(\mathbf{r})P_1(\mathbf{r}')$. In a simulation assuming an ideal fluid, the correlations will be missed, that is, the fluid at small scales may erroneously be considered as ideal if correlations are present. This yields an erroneous estimate of the interaction

energy. The proper calculation of the interaction energy must thus be done accurately and ensure that the small-scale interaction energy is properly accounted for. This may be done through Eq. (3) and the use of a nonideal framework.

## 2.2. Virial theorem for correlated fluids

For a system of interacting particles of mass $m$, the N-body form of the virial theorem can be stated as

$$ -\sum_{i} m v_{i}^{2} + \frac{1}{2} \sum_{i} \frac{d^{2} I_{i}}{dt^{2}} = \sum_{i, j, i \neq j} \mathbf{F}_{j \to i} \cdot \mathbf{r}_{i}, \qquad (7) $$

where $I_i = mr_i^2$ is the moment of inertia of the particle i. In a fluid (statistical) approach, the equivalent form reads

$$ -N \int_V m v^2(\mathbf{r}) P_1(\mathbf{r}) dV + \frac{1}{2} N \int_V \frac{d^2 I}{dt^2}(\mathbf{r}) P_1(\mathbf{r}) dV \\ = N(N-1) \iint_{V,V} \mathbf{F}_{\mathbf{r}' \to \mathbf{r}} \cdot \mathbf{r} P_2(\mathbf{r}, \mathbf{r}') dV dV'. \qquad (8) $$

Assuming a force deriving from a potential $\phi \propto |\mathbf{r} - \mathbf{r}'|^{\alpha}$, the last integral can be symmetrized by using $P_2(\mathbf{r}, \mathbf{r}') = P_2(\mathbf{r}', \mathbf{r})$ to obtain

$$ N(N-1) \iint_{V,V} \mathbf{F}_{\mathbf{r}' \to \mathbf{r}} \cdot \mathbf{r} P_2(\mathbf{r}, \mathbf{r}') dV dV' \\ = \frac{N(N-1)}{2} \iint_{V,V} \mathbf{F}_{\mathbf{r}' \to \mathbf{r}} \cdot (\mathbf{r} - \mathbf{r}') P_2(\mathbf{r}, \mathbf{r}') dV dV' \\ = -\alpha \frac{N(N-1)}{2} \iint_{V,V} \phi(|\mathbf{r} - \mathbf{r}'|) P_2(\mathbf{r}, \mathbf{r}') dV dV', \qquad (9) $$

in which we recognize the interaction energy defined in Eq. (3). The virial theorem now takes the classic form

$$ \langle mv^2 \rangle - \alpha \langle H_{\text{int}} \rangle = \frac{\langle d^2 I / dt^2 \rangle}{2}, \qquad (10) $$

in which we note the presence of the pair correlation function in the interaction energy.

Similarly, we can define the total mechanical energy $E_m$ of the N-body ensemble of particles,

$$ E_m = \langle mv^2/2 \rangle + \langle H_{\text{int}} \rangle \\ = \frac{1}{2} \sum_i m v_i^2 + \sum_{i<j} \phi(|\mathbf{r}_i - \mathbf{r}_j|), \qquad (11) $$

whose fluid counterpart is given by

$$ E_m = \frac{N}{2} \int_V m v^2(\mathbf{r}) P_1(\mathbf{r}) dV + \frac{N(N-1)}{2} \iint_{V,V} \phi(|\mathbf{r} - \mathbf{r}'|) P_2(\mathbf{r}, \mathbf{r}') dV dV'. \qquad (12) $$

If we neglect correlations in the interaction energy, then $P_2(\mathbf{r}, \mathbf{r}') = P_1(\mathbf{r})P_1(\mathbf{r}')$, which in the limit $N \gg 1$ yields

$$ \begin{aligned} \langle H_{\text{int}} \rangle &\approx \frac{N^2}{2} \iint_{V,V} \phi(|\mathbf{r} - \mathbf{r}'|) P_1(\mathbf{r}) P_1(\mathbf{r}') dV dV' \\ &\approx \frac{1}{2} \iint_{V,V} \phi(|\mathbf{r} - \mathbf{r}'|) \rho(\mathbf{r}) \rho(\mathbf{r}') dV dV' \\ &\approx \frac{1}{2} \int_V \rho(\mathbf{r}) \Phi(\mathbf{r}) dV, \end{aligned} \qquad (13) $$
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Schematic illustrating the interactions between two fluid volumes used to derive Newton's laws for pressureless correlated fluids.</figcaption></figure>

where we have defined the mean-field potential $\Phi(r) = \int_V \phi(|r - r'|)\rho(r')dV'$, which can be computed with the Poisson equation for the gravitational or the electrostatic force, $\nabla^2\Phi(r) \propto \rho(r)$. However, we stress that the possibility of defining this mean-field potential is only possible in the absence of correlations: as soon as $P_2(r, r') \neq P_1(r)P_1(r')$, the use of a mean-field and the Poisson equation is not correct as it will ignore the correlations present in the system.

### 2.3. Newton's laws of motion for pressureless correlated fluids

At a more fundamental level, we can also show that correlations should be taken into account in Newton's laws for a (pressureless) correlated fluid. We consider now an ensemble of $N_1$ particles of mass $m_1$ in a volume $V_1$ interacting with $N_2$ particles of mass $m_2$ in a volume $V_2$ (see Fig. 2). In an N-body approach, the second law can be written

$$ \sum_i m_1 a_i = \sum_{j,i} F_{j \to i}, \quad (14) $$

with $i$ in $[0, N_1]$ and $j$ in $[0, N_2]$. Newton's second law for a pressureless correlated fluid is then given by

$$ N_1 \int_{V_1} m_1 a(r) P_1(r) dV_1 = N_1 N_2 \iint_{V_1, V_2} F_{r_2 \to r_1} P_2(r_1, r_2) dV_1 dV_2, \quad (15) $$

which we write in a more compact way as $\langle m_1 a \rangle_1 = \langle F \rangle_{2 \to 1}$. We have generalized the two-point probability density function to different volumes $V_1$ and $V_2$, and it is normalized by $N_1 N_2$.

The first law for an isolated fluid of volume $V_1$ is then $\langle m_1 a \rangle_1 = \langle F \rangle_{1 \to 1} = 0$ and can be deduced from the second law by taking $V_2 = V_1$ and using $P_2(r_1, r_1') = P_2(r_1', r_1)$ and $F_{r_1 \to r_1'} = -F_{r_1' \to r_1}$. Similarly, the third law $\langle F \rangle_{2 \to 1} = -\langle F \rangle_{1 \to 2}$ can be obtained using the same arguments.

We now assume that there is no correlation, that is, $P_2(r_1, r_2) = P_1(r_1)P_1(r_2)$, and that the fluid is isotropic and homogeneous: $P_1(r_1) = 1/V_1$. Choosing the origin of positions, $\mathbf{x}$, at the center of $V_1$, we have $N_1/V_1 = \rho_1(\mathbf{x})$ (see Sect. 2.1), and we define

$$ \begin{aligned} a(\mathbf{x}) &= \frac{1}{V_1} \int_{V_1} a(\mathbf{r}) dV_1, \\ \Phi(\mathbf{x}) &= \frac{1}{V_1} \iint_{V_1, V_2} \phi(|\mathbf{r}_2 - \mathbf{r}_1|) \rho_2(\mathbf{r}_2) dV_1 dV_2. \end{aligned} \quad (16) $$

The second law then becomes

$$ \begin{aligned} \frac{\rho_1(\mathbf{x})m_1}{V_1} \int_{V_1} a(\mathbf{r}) dV_1 &= \frac{\rho_1(\mathbf{x})m_1}{V_1} \iint_{V_1, V_2} \frac{F_{r_2 \to r_1}}{m} \rho_2(\mathbf{r}_2) dV_1 dV_2 \\ \rho_1(\mathbf{x})m_1 a(\mathbf{x}) &\approx \rho_1(\mathbf{x})m_1 \frac{-\nabla_x \Phi(\mathbf{x})}{m_1}, \end{aligned} \quad (17) $$

where we recognize the usual form of an external mean force deriving from the mean-field potential $\Phi(\mathbf{x})$. In the case of the gravitational force, $\Phi(\mathbf{x})/m_1$ is independent of $m_1$, and the equivalence principle for a pressureless ideal fluid simply consists of simplifying Eq. (17) by $\rho_1(\mathbf{x})m_1$ on each side,

$$ a(\mathbf{x}) \approx \frac{-\nabla_x \Phi(\mathbf{x})}{m_1}. \quad (18) $$

However, this simplification is only possible in the absence of correlations. That is, as soon as $P_2(r, r') \neq P_1(r)P_1(r')$, the use of the equivalence principle is not correct as it will ignore the correlations present in the system.

### 2.4. Nonideal equation of state

We recall that for a homogeneous and isotropic fluid, $P_2(r, r') = g(|r - r'|)/V^2$. In this case (for $N \gg 1$), the interaction energy becomes

$$ \langle H_{\text{int}} \rangle = \frac{2\pi N^2}{V} \int_0^R g(r) \phi(r) r^2 dr, \quad (19) $$

with $r = |r - r'|$ and $V = 4\pi R^3/3$. Statistical mechanics in the canonical ensemble in thermal equilibrium with a heat bath at a fixed temperature $T$ can be used to define the EOS. We refer to classical textbooks such as Hansen et al. (2006) for details. The energy of the system in the absence of external forces is then given by

$$ E = \frac{N k_B T}{\gamma - 1} + \frac{2\pi N^2}{V} \int_0^R g(r) \phi(r) r^2 dr, \quad (20) $$

with $\gamma = C_p/C_V$ the adiabatic index of the fluid, defined as the ratio of the specific heats at constant pressure and volume, respectively ($\gamma = 5/3$ for mono-atomic particles in their ground state). The first term in Eq. (20) is the perfect gas contribution, and the second term stems from the interactions between particles. The nonideal pressure is given by

$$ PV = N k_B T - \frac{2\pi N^2}{3V} \int_0^R g(r) \frac{d\phi(r)}{d} r r^3 dr. \quad (21) $$

When the temperature is identified with the velocity dispersion of particles within the fluid ($k_B T = 1/2\langle mv_i^2 \rangle$), Eq. (21) can be interpreted as the isotropic form of the virial theorem. $P = 0$ corresponds to virial equilibrium, for which the ideal pressure is balanced by the interaction contribution. $P > 0$ corresponds to a fluid that expands because either the ideal (kinetic) pressure dominates the interaction term or the interactions are repulsive, $d\phi(r)/dr < 0$. In contrast, $P < 0$ corresponds to a fluid that collapses because the contribution due to attractive interactions $d\phi(r)/dr > 0$ dominates the ideal pressure. This case corresponds to a phase transition from gas to liquid, for instance, with the attractive interactions being the intermolecular dipole forces, or to gravitational collapse in the case of the gravitational force (see, e.g., de Vega & Sánchez 2002).

Strictly speaking, a self-interacting fluid should always be described with a nonideal EOS. Nevertheless, an ideal approximation is possible if the interaction term is negligible in Eq. (21).
In this case, $g(r) \approx 1$, so that with $\rho = N/V$, the condition becomes

$$ \rho k_B T \gg \frac{2\pi\rho^2}{3} \int_0^R \frac{d\phi(r)}{d} rr^3 dr. \qquad (22) $$

For a potential of the form $\phi = -k/r$ (with $k > 0$ for the gravitational force or attractive Coulomb interactions),

$$ \rho k_B T \gg \pi \rho^2 k R^2 / 3 \Rightarrow R \ll \sqrt{\frac{k_B T}{k \rho}}. \qquad (23) $$

In this case, we recognize the characteristic Debye screening length in the case of Coulomb interactions or the Jeans length in the case of the gravitational force. The condition in Eq. (23) states that an ideal approximation in the treatment of the fluid under consideration is acceptable if (and only if) the volumes considered are sufficiently small so that the Debye length or the Jeans length is well resolved. Whereas such a condition is usually fulfilled in numerical simulations of star formation and is known as the Truelove condition (see Truelove et al. 1997), it is not easily done for large-scale numerical simulations of galaxies or larger structures. In these cases, it is mandatory to use a non-ideal EOS and to include the terms involving the pair-correlation function in the calculations in this way.

## 2.5. Nonideal stress tensor

More generally, we express the above theory in the form of a stress tensor, following the pioneering work of Kirkwood et al. (1949). The Newtonian expression for the stress tensor of a system of N particles of mass m is given by

$$ \sigma_{ij} = - \left( P + \left( \frac{2}{3}\eta - \chi \right) \sum_k \partial_k u_k \right) \delta_{ij} + \eta \left( \partial_i u_j + \partial_j u_i \right), \qquad (24) $$

with P the nonideal pressure given in Sect. 2.4, and $\eta$ and $\chi$ the coefficients of shear and bulk viscosity, which can be computed as

$$ \eta = \frac{1}{2} \frac{N}{V} m D + \frac{\pi}{15 D} \frac{N^2}{V^2} \int_0^R \psi_2(r) g(r) \frac{d\phi(r)}{d} rr^3 dr, $$

$$ \chi = \frac{1}{3} \frac{N}{V} m D + \frac{\pi}{9 D} \frac{N^2}{V^2} \int_0^R \psi_0(r) g(r) \frac{d\phi(r)}{d} rr^3 dr, \qquad (25) $$

with D the self-diffusion coefficient in the fluid (from Fick's law). As defined in Kirkwood et al. (1949), $\psi_0(r)$ is the surface harmonic of zeroth order arising from the dilatation component of the rate of strain, and $\psi_2(r)$ is the surface harmonic of order 2 arising from the shear component. They obey the following differential equations:

$$ \begin{aligned} \frac{d}{dr} \left( r^2 g(r) \frac{d\psi_2(r)}{dr} r \right) - 6\psi_2(r) g(r) &= r^3 \frac{dg(r)}{dr} \\ \frac{d}{dr} r \left( r^2 g(r) \frac{d\psi_0(r)}{dr} \right) &= r^3 \frac{dg(r)}{dr}. \end{aligned} \qquad (26) $$

As for the nonideal pressure, the nonideal coefficients of shear and bulk viscosity contain two contributions: an ideal (Brownian) contribution arising from the momentum transfer between colliding particles, and a nonideal part arising from the direct transfer by interaction forces. The second part strongly dominates when interactions are present in the fluid. As for the nonideal pressure, it is also crucial to take interactions in the stress tensor into account when the Jeans length is not resolved.

## 2.6. Thermodynamic limit

We can take the thermodynamic limit by imposing $(N, V) \to \infty$ while keeping $N/V = \rho$ constant. Defining $\rho_m e$ the internal energy per unit volume (with $\rho_m = \rho m$),

$$ P = \rho k_B T - \frac{2\pi\rho^2}{3} \int_0^\infty g(r, \rho, T) \frac{d\phi(r)}{dr} r^3 dr, \qquad (27) $$

$$ \rho_m e = \frac{\rho k_B T}{\gamma - 1} + 2\pi\rho^2 \int_0^\infty g(r, \rho, T) \phi(r) r^2 dr, \qquad (28) $$

where we have made the dependence of the radial distribution function on the external variables $\rho$ and $T$ explicit. Determining the EOS now reduces to the determination of $g(r, \rho, T)$, which can be done by three different technics in the case of liquids: calculation of the Ornstein-Zernike equation, N-body Monte Carlo or molecular dynamic simulations, or experiments by X-ray or neutron diffraction (see Aslangul 2006). We emphasize that we implicitly assume here that an equilibrium can be reached. This allows us to define a radial distribution function that does not depend on time and can be determined from external variables ($\rho$ and $T$). Strictly speaking, this might not be the case for self-gravitating systems that are collapsing, hence are not in equilibrium. This assumption then assumes that some sort of feedback stabilizes small scales, allowing for the definition of an EOS for the scales that are not resolved and described by statistical mechanics. This assumption can be relaxed when the dynamical equations in the BBGKY hierarchy are used that describe the evolution of the correlations. However, as discussed in Davis & Peebles (1977), we need to find a closure relation for the hierarchy by expressing the last N-point correlation functions as a function of the others. This closure is likely possible only at very high N in order to accurately describe the structures that are observed in the Universe (see Balian & Schaeffer 1989a,b), which would result in a complex system of dynamical equations describing the evolution of correlations. We therefore use here a pragmatic approach (similarly to the approach used for liquids in physical chemistry) and assume that an equilibrium, or more exactly, a stationary state, exists at small scales. We also assume that we can determine the two-point correlation function either by explicitly resolving inhomogeneities in small-scale numerical simulations (e.g., with molecular dynamics for liquids or N-body simulations for self-gravitating systems) or by using astrophysical observations.

An important point is that the thermodynamic limit in Eq. (28) is well defined only if we obtain an extensive definition of the total internal energy $E$ (or intensive for the internal energy density per unit volume $\rho_m e$), that is, the integral in the expression must converge when $N, V \to \infty$. Because $\lim_{r \to \infty} g(r) = 1$, this means that the potential from which the force derives must decrease faster than $r^{-4}$ at infinity. This is the case for the Lennard-Jones potential, which is classically used to model intermolecular forces in liquid water or molecular liquids. An interpretation of this requirement is that the interaction energy in Eq. (20) is not proportional to $N^2$ but to $N \times N_{\text{neighbors}}$, with $N_{\text{neighbors}}$ the number of particles involved in the short-range interactions. Therefore, the interaction energy is proportional to $N$ and the definition of $E$ is extensive in the thermodynamic limit. If the interactions are long range, however, the interaction energy is really proportional to $N^2$. In this case, the total energy $E$ is not extensive, and we cannot use the thermodynamic limit. This problem is usually thought to prevent the use of statistical mechanics to describe systems characterized by a long-range force, such as gravity. We show in Sect. 3 a correct treatment in the case of the gravitational force.
## 3. Nonideal self-gravity

### 3.1. Hierarchical virial theorem

As we showed in Sect. 2.4, characteristic length scales such as the Debye length $\lambda_D$ or Jeans length $\lambda_J$ can be defined to determine whether nonideal (interaction) contributions should be taken into account. In the case of the electrostatic force, the Coulomb potential is split into a short-range and a long-range component, $\phi = \exp(-r/\lambda_D)\phi + (1 - \exp(-r/\lambda_D))\phi$. The first part corresponds to the screened Coulomb potential, due to the polarizable electron background, and the residual corresponds to the departure between the Coulomb and the screened potential (e.g., Chabrier 1990). We propose to decompose the gravitational force in a similar manner.

Hydrogen atoms become gravitationally bounded in the interior of a star, for example. Using a classical polytropic stellar model, we can define a virial theorem that balances the internal energy in the star and the gravitational interaction energy (see, e.g., Sect. 4.1). If now we consider a cluster of interacting polytropic stars, we can consider each star as isolated in a first approximation and use the previous virial theorem to obtain their internal structure and then describe their collective dynamics assuming they are point masses (i.e., we neglect tidal interactions). However, we can easily see that we should not integrate the gravitational force at infinity when using the virial theorem in a single star, otherwise we would include some gravitational interactions in the internal stellar structures that contribute to the cluster dynamics. A solution for this issue is to decompose the gravitational potential into a near-field and far-field component using a length scale $\lambda_0$ larger than the size of a star,

$$ \phi = e^{-r/\lambda_0}\phi + (1 - e^{-r/\lambda_0})\phi, \quad (29) $$

and we apply the virial theorem inside a star using the near-field $\phi_0 = \exp(-r/\lambda_0)\phi$. We emphasize that the use of exponential damping is arbitrary. Another possibility is to use an erf function, as is often done with the Coulomb interaction in physical chemistry. We can then proceed recursively in the residual part. We now use a second length scale $\lambda_1$ in order to define a near field $\phi_1$ for the virial theorem describing the dynamics inside the stellar cluster and a far field for the dynamics of a collection of stellar clusters inside a galaxy, for example. $\phi_1$ is then given by $\exp(-r/\lambda_1)(1 - \exp(-r/\lambda_0))\phi$. We can continue to define the gravitational energy contributing to the dynamics of the stellar clusters inside the galaxy and then the dynamics of the galaxies inside a galaxy cluster, and so on. The hierarchical decomposition of the gravitational potential with a collection of length scales $\lambda_0, \dots, \lambda_i, \dots$ is then given by

$$ \phi = \sum_{i=0}^{\infty} \phi_i \\ \phi_i = e^{-r/\lambda_i} \prod_{j=0}^{i-1} (1 - e^{-r/\lambda_j}) \phi. \quad (30) $$

This expansion is similar to the BBGKY hierarchy for the N-body problem in statistical physics, where the potential is split into a (short-range) interparticle pair potential and a (long-range) external-field potential (see, e.g., Chavanis 2013). The series of length scales is a-priori arbitrary, and for a purely self-gravitating collapsing fluid, there is no particular length scale that would physically make sense. The characteristic length scales that should be used are likely linked to feedback and support processes that define the possible virial equilibria at different scales:

active galactic nucleus feedback for galaxy clusters, rotational support for galaxies, stellar feedback for molecular clouds and star clusters, and pressure support for stars. When this stationary state is reached, for instance, inside a star, we can account for the velocity dispersion and gravitational energy at this scale as internal degrees of freedom (similarly to an EOS) when going to a larger scale. When tidal effects are neglected, this allows decoupling the different scales based on feedback processes.

We can now use the virial theorem defined in Sect. 2.2 to link the dynamic properties $v_i, I_i$ of objects of mass $m_i$ at each scale $\lambda_i$ of the hierarchy to the corresponding gravitational interaction potential $\phi_i$,

$$ -\langle m_i v_i^2 \rangle + \langle d^2 I_i / dt^2 \rangle / 2 = \langle F_i \cdot r \rangle, \quad (31) $$

in which we can replace kinetic energy by thermal energy for $i=0$. As mentioned above, we emphasize that this hierarchy neglects tidal interactions. A possible extension would be to also decompose the energy going into tides, if they are not negligible at a given scale. This is likely to be the case at the scale of galaxies (see, e.g., Barnes & Hernquist 1992). For the time being, however, we assume that these contributions represent a small correction to the total gravitational energy.

Because each term $\langle F_i \cdot r \rangle$ is computed with the short-range potential $\phi_i$, damped by a factor $e^{-r/\lambda_i}$, we can now safely take the thermodynamic limit in Eq. (31), that is, $N, V \to \infty$, because the integral converges. For a homogeneous and isotropic correlated fluid, this term is given by

$$ \langle F_i \cdot r \rangle = 2\pi\rho^2 \int_0^\infty g(r) \frac{d\phi_i(r)}{dr} r^3 dr. \quad (32) $$

Similarly, the total mechanical energy of the system at a scale $\lambda_i$ is given by

$$ E_{m,i} = \langle m_i v_i^2 / 2 \rangle + \langle H_{\text{int},i} \rangle, \quad (33) $$

where we can take the thermodynamic limit for the interaction energy,

$$ \langle H_{\text{int},i} \rangle = 2\pi\rho^2 \int_0^\infty g(r) \phi_i(r) r^2 dr. \quad (34) $$

We show in the next section that a similar decomposition can be used to define nonideal self-gravitating hydrodynamics in numerical simulations.

### 3.2. Nonideal self-gravitating hydrodynamics

A standard strategy to obtain a large-scale hydrodynamic model of complex multiphasic systems relies on the homogenization procedure by averaging volumes of characteristic size $\lambda$ (Whitaker 1998). We can then decompose the interaction potential into a near-field and a far-field contribution using this characteristic scale. However, in numerical simulations, as soon as the grid size (or the kernel extension in SPH methods) is such that $\Delta x < \lambda$, inhomogeneities at scales smaller than $\lambda$ have no physical meaning in such a homogenized model, which is only valid at scales larger than $\lambda$.

In numerical simulations, it is then natural to use the grid resolution $\Delta x$ as the characteristic length to split the gravitational interaction

$$ \phi(r) = e^{-r/\Delta x}\phi(r) + (1 - e^{-r/\Delta x})\phi(r). \quad (35) $$

We can then define a near field $\phi_{\text{int}}(r) = \exp(-r/\Delta x)\phi(r)$ to handle the nonideal effects in the EOS inside the simulation control
volumes, and a far field $\phi_{\text{ext}}(r) = (1 - \exp(-r/\Delta x))\phi(r)$ to define the external force between the simulation control volumes. Following Irving & Kirkwood (1950), the Navier-Stokes equations for a nonideal fluid can be derived from the Liouville equation in the BBGKY hierarchy, and in the case of a self-gravitating fluid, they are given by

$$ 
\begin{aligned}
& \frac{\partial \rho_m}{\partial t} + \nabla (\rho_m \mathbf{u}) = 0, \\
& \frac{\partial \rho_m \mathbf{u}}{\partial t} + \nabla (\rho_m \mathbf{u} \otimes \mathbf{u} + \sigma) = \mathbf{F}_{\text{ext}}, \\
& \frac{\partial \rho_m E}{\partial t} + \nabla (\rho_m \mathbf{u} E + \sigma \cdot \mathbf{u}) = \mathbf{F}_{\text{ext}} \cdot \mathbf{u},
\end{aligned}
\quad (36) $$

using the following relations for $E$ and $\sigma$

$$ 
\begin{aligned}
E &= e + \frac{1}{2} \mathbf{u}^2, \\
\sigma_{ij} &= - \left( P + \left(\frac{2}{3}\eta - \chi\right) \sum_k \partial_k u_k \right) \delta_{ij} + \eta \left( \partial_i u_j + \partial_j u_i \right),
\end{aligned}
\quad (37) $$

in which the nonideal (internal) quantities are given by

$$ 
\begin{aligned}
P &= \rho k_B T - \frac{2\pi\rho^2}{3} \int_0^\infty g(r) \frac{d\phi_{\text{int}}(r)}{dr} r^3 dr \\
\rho_m e &= \frac{\rho k_B T}{\gamma - 1} + 2\pi\rho^2 \int_0^\infty g(r) \phi_{\text{int}}(r) r^2 dr \\
\eta &= \frac{1}{2} \rho_m D + \frac{\pi}{15D} \rho^2 \int_0^\infty \psi_2(r) g(r) \frac{d\phi_{\text{int}}(r)}{dr} r^3 dr \\
\chi &= \frac{1}{3} \rho_m D + \frac{\pi}{9D} \rho^2 \int_0^\infty \psi_0(r) g(r) \frac{d\phi_{\text{int}}(r)}{dr} r^3 dr \\
\text{with } \phi_{\text{int}}(r) &= -\frac{Gm^2}{r} e^{-r/\Delta x},
\end{aligned}
\quad (38) $$

with $\psi_2(r)$ and $\psi_0(r)$ given by Eq. (26). The external force is defined by an integral form in the whole simulation domain $V_{\text{sim}}$

$$ 
\begin{aligned}
F_{\text{ext}}(\mathbf{x}) &= -\rho(\mathbf{x})\nabla_{\mathbf{x}}(\Phi_{\text{ext}}(\mathbf{x})), \\
\Phi_{\text{ext}}(\mathbf{x}) &= \int_{V_{\text{sim}}} \phi_{\text{ext}}(|\mathbf{x} - \mathbf{x}'|)\rho(\mathbf{x}')dV_{\text{sim}}, \\
\text{with } \phi_{\text{ext}}(r) &= -\frac{Gm^2}{r}(1 - e^{-r/\Delta x}).
\end{aligned}
\quad (39) $$

As for the virial theorem in the previous section, all the integrals in the EOS quantities are defined with the damped potential $\phi_{\text{int}}(r) = e^{-r/\Delta x}\phi(r)$. Hence, all these integrals are well defined and converge in the thermodynamic limit. When $\Delta x \ll \lambda_J$ and notably in the limit $\Delta x \to 0$, we have $\phi_{\text{int}}(r) \to 0$ and $\phi_{\text{ext}}(r) \to \phi(r)$, which means that all the gravitational force is entirely resolved externally and the nonideal contributions to the fluid EOS given by Eq. (38) are negligible: We recover then the classical (ideal) form of the Navier-Stokes equations for a self-gravitating fluid. As soon as $\Delta x > \lambda_J$, however, the nonideal contributions to the EOS due to interactions within the control volumes must be included for a correct calculation of the gravitational energy of the whole system.

# 4. Applications

## 4.1. Impact of inhomogeneities on the virial theorem

The virial theorem is usually applied in astrophysics in a simplified way. We define the total mass $M = Nm$, and the average total energy $\langle E/M \rangle$, kinetic energy $\langle E_{\text{kin}}/M \rangle$, and potential

energy $\langle E_p/M \rangle$ per unit mass. The virial theorem, assuming that the particles are gravitationally bound in a sphere of radius $r_{\text{bound}}$, is usually evaluated with these quantities per unit mass

$$ 
\begin{aligned}
\langle E/M \rangle &= \langle E_{\text{kin}}/M \rangle + \langle E_p/M \rangle, \\
-2\langle E_{\text{kin}}/M \rangle &= \langle E_p/M \rangle \approx -\frac{GM}{r_{\text{bound}}}. 
\end{aligned}
\quad (40) $$

The use of a radius $r_{\text{bound}}$ is similar to the decomposition into near and far fields that we introduced in Sect. 3.1. It defines a scale up to which the gravitational energy is considered to contribute to the virial balance, while the rest is implicitly assumed to contribute to larger-scale dynamics. However, this form of the virial theorem contains no information about the substructures present within the bound volume: only the total mass and size of the system enters the virial equation. This latter thus ignores nonideal effects due to correlations induced by these substructures.

In order to show the importance of inhomogeneities in the virial theorem, we apply it as a simple semianalytical example to the case of a polytropic stellar structure. In this case, we can explicitly compute $P_2(r, r')$ by semianalytically computing an inhomogeneous density field and compare the result with a homogeneous assumption on the virial theorem. A polytropic stellar structure can be calculated by the Lane-Emden equation,

$$ \frac{1}{z^2} \frac{d}{dz} \left( z^2 \frac{dw}{dz} \right) + w^n = 0, \quad (41) $$

with $\rho = \rho_c u^n$, $P = K\rho^\gamma$, $n = 1/(\gamma - 1)$ the polytropic index, $z = Ar$, and $A^2 = (4\pi G/(n+1))\rho_c^{1-1/n}/K$. We choose a dimensionless unit system in which $\rho_c = K = G = 1$. By solving the Lane-Emden equation for different polytropic indexes $n < 5$, we can compute inhomogeneous density profiles that have a finite radius $R$ (see the top panel of Fig. 3). For a polytrope, the following virial theorem (per unit mass) reads

$$ 
\begin{aligned}
H/M &= 2H_{\text{int}}/M, \\
H &= \int_V \frac{\gamma}{\gamma - 1} P dV, \\
M &= \int_V \rho dV,
\end{aligned}
\quad (42) $$

where $H$ is the total enthalpy, $M$ is the total mass, and $H_{\text{int}}$ is given by Eq. (13) with the zero of potential energy defined at the surface of the star. For a homogeneous density profile in Eq. (13), with $\rho = M/(4\pi R^3/3)$, we obtain

$$ \frac{H}{M} = \frac{GM}{5R}. \quad (43) $$

For the inhomogeneous case, as illustrated in Fig. 1, there are two possibilities for the description of the inhomogeneities. In the first case (with $V_1$), they are resolved (e.g., in high-resolution simulations or observations) so that $\rho(r)$ is not constant. In the opposite case (with $V_2$), $\rho(r)$ is constant and the (unresolved) inhomogeneities are included in the pair correlation function $g(r)$ in our formalism. Determining the impact of ignoring the inhomogeneities thus leads to two different possibilities. For the resolved case, it means comparing the computation with $\rho(r)$ not constant with the computation in which $\rho(r)$ is replaced by its mean value. For the unresolved case, it means comparing the computation with $g(r) \neq 1$ with the copmutation assuming an uncorrelated fluid, $g(r) = 1$. The semianalytical polytropic model allows us to easily perform the first (resolved)
<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Polytropic solutions. Top: density profiles solutions to the Lane-Emden equation with different polytropic indexes n. Bottom: total enthalpy per unit mass compared to gravitational interaction energy per unit mass assuming either a homogeneous density profile or accounting for inhomogeneities. Blue and green crosses are indistinguishable.</figcaption></figure>

test. Because we can explicitly compute the inhomogeneous profile $\rho(r)$ and all the necessary integral quantities, we can compare this exact solution with the solution in which $\rho(r)$ is replaced by its mean value.

We now assume that an observer has only access to the total enthalpy, mass, and radius of different objects and does not know whether these objects have substructures. The naive assumption, as explained above, is to ignore the possible inhomogeneities and assume that the density is constant, $\rho = M/(4\pi R^3/3)$, hence using Eq. (43) for the virial equation. This assumption, however, is only correct for a polytropic index $n = 0$. As shown in the bottom panel of Fig. 3, for profiles $n \neq 0$, the error on the total energy is significant and can be up to a factor 25 for $n = 4$. Therefore, the assumption of a homogeneous profile within the structures will lead to a missing-mass problem. For a polytropic structure, however, we can calculate the inhomogeneous density profile semianalytically and thus calculate exactly the contribution of inhomogeneities to the total energy. We characterize this contribution by an inhomogeneous amplification factor $\alpha_{ih}$, which can be calculated exactly,

$$ \alpha_{ih} = \frac{N^2 \iint_{V,V} P_2(\mathbf{r}, \mathbf{r}') \phi(|\mathbf{r} - \mathbf{r}'|) dV dV'}{N^2 \iint_{V,V} 1/(4\pi R^3/3)^2 \phi(|\mathbf{r} - \mathbf{r}'|) dV dV'} \\ = \frac{\iint_{V,V} \rho(r) \rho(r') (1/|\mathbf{r} - \mathbf{r}'| - 1/|\mathbf{R} - \mathbf{r}'|) dV dV'}{\iint_{V,V} (M/(4\pi R^3/3))^2 (1/|\mathbf{r} - \mathbf{r}'| - 1/|\mathbf{R} - \mathbf{r}'|) dV dV'} \quad (44) $$

with $\mathbf{R}$ a vector on a sphere of radius $R$ (the radius of the polytropic structure), used to define the zero of the gravitational potential at the stellar surface, and $V = 4\pi R^3/3$. In our calculation, $N^2 P_2(\mathbf{r}, \mathbf{r}') = N^2 P_1(\mathbf{r}) P_1(\mathbf{r}') = \rho(\mathbf{r})\rho(\mathbf{r}')$ for the inhomogeneous structure, while the homogeneous calculation corresponds to $N^2 P_2(\mathbf{r}, \mathbf{r}') = N^2 P_1(\mathbf{r}) P_1(\mathbf{r}') = N^2/V^2$, with $P_1(\mathbf{r}) = 1/V$. We can now correct Eq. (43) to account for inhomogeneities,

$$ \frac{H}{M} = \frac{GM}{5R} \alpha_{ih}. \quad (45) $$

As shown in the bottom panel of Fig. 3, when this factor is taken into account, we exactly recover the correct enthalpy and the correct gravitational interaction energy for the different polytropic structures (the green and blue crosses in the figure are indistinguishable). This simple example demonstrates that it is mandatory to account for the contribution of substructures to correctly evaluate interaction energies in astrophysical structures.

## 4.2. Gravitational viscous regime and the rotation curve of galaxies

Another consequence of substructures and interactions is a possible transition to a viscous regime when interactions dominate at small scales. Heuristically, the rotation curve of a fluid is expected to be strongly impacted by viscous stresses. For high viscosity, the rotation tends to the rotation of a Taylor-Couette flow with $u_\theta \rightarrow \Omega r$ in cylindrical coordinates.

We can explore different regimes by performing a Hilbert expansion with a small parameter $\epsilon$ that characterizes the regime of the flow, that is, we can take $\epsilon$ equal to the inverse of the Reynolds number $Re$ such that $\epsilon \rightarrow 0$ in the inertial limit. We assume a stationary 2D stellar flow in cylindrical coordinates using either stellar hydrodynamics (see Binney & Tremaine 1987; Burkert & Hensler 1988) or gas hydrodynamics in the HI disk with the nonideal momentum evolution equation presented in Sect. 3.2. In the case of stellar hydrodynamics, we can replace $k_B T$ in Eq. (36) by the velocity dispersion tensor, as is done in the Jeans equation, while the interaction part in the nonideal term remains unchanged. We develop all the variables with respect to $\epsilon$,

$$ \begin{align} \rho_m(r, \theta) &= \rho_m^0(r) + \epsilon \rho_m^1(r, \theta) + \dots \notag \\ u_r(r, \theta) &= \epsilon u_r^1(r, \theta) + \dots \notag \\ u_\theta(r, \theta) &= u_\theta^0(r) + \epsilon u_\theta^1(r, \theta) + \dots \notag \\ P(r, \theta) &= P^0(r) + \epsilon P^1(r, \theta) + \dots \notag \\ \mathbf{F}_{\text{ext}}(r, \theta) &= \mathbf{F}_{\text{ext}}^0(r) + \epsilon \mathbf{F}_{\text{ext}}^1(r, \theta) + \dots \notag \\ \eta &= \eta^0 + \epsilon \eta^1 + \dots, \tag{46} \end{align} $$

where we have made the hypothesis that the zeroth order is axisymmetric with $u_\theta^0 = 0$ and verifies the condition $\nabla(u^0) = 0$.

We can now explore the inertial regime by assuming the following dependence of the shear viscosity and pressure, $\eta^0 \approx 0$ and $P^0(r) \approx 0$ in the ideal regime, thus with no contribution from interactions. The zeroth order of the nonideal Navier-Stokes equation (Eq. (36)) then gives

$$ \begin{align} 0^{\text{th}}, \text{r-component: } \frac{-(u_{\theta}^{0})^{2}}{r} &= \frac{1}{\rho_{m}^{0}} F_{\text{ext},r}^{0} \notag \\ 0^{\text{th}}, \theta\text{-component: } \theta = 0, \tag{47} \end{align} $$

which gives the classical velocity profile $u_\theta^0 = \sqrt{GM/r}$, decreasing as $r^{-1/2}$ at large distance.
In the viscous regime dominated by the nonideal contribution, we now have $\eta^0 \neq 0$ and $P^0(r) \neq 0$ because of the presence of interactions and bound structures. The zeroth order of the non-ideal Navier-Stokes equation then gives

$$ 
\begin{aligned}
0^{\text{th}}, r\text{-component}: \frac{-(u_{\theta}^0)^2}{r} &= -\frac{1}{\rho_m^0} \frac{\partial P^0}{\partial r} + \frac{1}{\rho_m^0} F_{\text{ext},r}^0 \\
0^{\text{th}}, \theta\text{-component}: 0 &= \eta^0 \left( \frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial u_{\theta}^0}{\partial r} \right) - \frac{u_{\theta}^0}{r^2} \right).
\end{aligned}
\quad (48) $$

We can obtain the velocity profile from the $\theta$ component of the Navier-Stokes equation. Under appropriate boundary conditions, the fluid can exhibit a Taylor-Couette flow in the viscous limit: $u_{\theta}^0 = \Omega r$. The associated centrifugal force in the $r$ component is then compensated for by pressure gradients. We recall that the nonideal pressure contains attractive gravitational interactions at small scale here, hence it can compensate for the centrifugal outward acceleration.

In principle, it should be possible to verify whether this interpretation is correct by exploring such regimes with N-body numerical simulations. The result of this study is presented in Fig. 4 from the N-body simulations of Voglis et al. (2006), who increased the total angular momentum at fixed total mass $M$ and size $R$, which means that the gravitational energy $GM/R$ is constant for all the simulations. The spiral-arm structures emerge on a scale $\lambda_{\text{arm}}$ that is smaller than $R$. When the full N-body dynamics cannot be explored and the simulation is approximated by a fluid, the structure bounded by the spiral arms indicates that the Jeans length is smaller than the size of the system. This implies that a nonideal description is needed when the spiral-arm structures have formed. In this case, the gravitational interactions at small scales inside the spiral arms have to be taken into account in the viscous stress, indicating a transition to the viscous regime. Figure 5 displays the rotation curves associated with the simulations shown in Fig. 4, assuming a half-mass radius of 3 kpc and a total mass of $10^{11} M_{\odot}$ for the galaxy. A deviation from the inertial regime as soon as spiral-arm structures form is evident. This corroborates the assumption that a transition to a viscous regime occurs.

These results also suggest that the change in the rotation curve between spiral and dwarf galaxies and some low surface brightness (LSB) elliptical galaxies exhibiting Keplerian profiles can be interpreted as a phase transition from a nonviscous to a viscous regime that depends on the presence of bound substructures in the galaxy (e.g., local clustering or spiral arms). An example of such a transition is illustrated in Fig. 6 with the observation data of Paolo et al. (2019).

### 4.3. Nonideal cosmology

As discussed in Sects. 2.2 and 2.3, the Poisson equation and the equivalence principle can be safely used only for uncorrelated fluids. The Poisson and Einstein equations are not designed to account for (short-scale) correlations that can be present in a fluid. The existence of large-scale structures that have formed during the evolution of the Universe means, however, that we are in the presence of a correlated fluid, even though at very large scales the Universe can be considered as homogeneous and isotropic. It might even be argued that at scales characteristic of these structures, the Universe is the most correlated fluid in nature. Figure 7 shows that the two-point correlation function and radial distribution function can reach values $g(r) \approx 10^8$, to be compared with $g(r) = 1$ for an ideal (uncorrelated) fluid (see

<figure><img src="image_5.png" /><figcaption>Fig. 4. Snapshots of the N-body simulations from Fig. 1 in Voglis et al. (2006). The arm or spiral structures are obtained by increasing the total angular momentum from models QR1 to QR4 for the same mass and size for all the simulations.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Rotation curves of the different models from QR1 to QR4 from Fig. 2 in Harsoula & Kalapotharakos (2009).</figcaption></figure>

Fig. 7). We show below how the effect of these correlations can be incorporated into the Friedmann equations. For simplicity, we restrict ourselves to the Newtonian limit and do not intend to explore an extension of general relativity to the nonideal regime at the moment. It can therefore be considered as a preliminary approach at this stage. As shown in many studies and textbooks (see, e.g., Peacock 1999), however, Newtonian cosmology has proved to be a very useful framework for exploring various physical effects in cosmology while keeping the calculations relatively simple.

Assuming a homogeneous and isotropic universe, we can obtain the well-known Friedmann equations from Einstein's equation,

$$ 
\begin{aligned}
H^2 + \frac{K}{a^2} &= \frac{8\pi G}{3c^2} \rho_e, \\
\dot{H} + \frac{3}{2}H^2 + \frac{K}{2a^2} &= -\frac{4\pi G}{c^2} P,
\end{aligned}
\quad (49) $$
<figure><img src="image_7.png" /><figcaption>Fig. 6. Relation between total acceleration and baryonic component for LSB and dwarf galaxies (Paolo et al. 2019). Acceleration deduced from the rotation curves of models QR1 and QR4 of Voglis et al. (2006) and Harsoula & Kalapotharakos (2009) assuming a half-mass radius of 3 kpc and a total mass of $10^{11} M_{\odot}$ for the galaxy. QR1 is taken as the reference for Keplerian rotation. The bifurcation between Keplerian rotation (QR1) and anomalous rotation (QR4) can be interpreted as a phase transition induced by the presence of substructures, i.e., spiral arms that have developed in model QR4.</figcaption></figure>

with $\rho_e$ and $P$ the energy density and pressure, $H = \dot{a}(t)/a(t)$ the Hubble parameter, $K$ the Gaussian curvature, and $a(t)$ the scale factor defining the geometry of the Universe, with $ds^2 = a(t)^2d\Sigma^2 - c^2dt^2$ and $d\Sigma^2 = dr^2/(1 - Kr^2) + r^2d\Omega^2$. We consider flat geometries, with $K = 0$, and we define the deceleration parameter $q = -1 - \dot{H}/H^2$. The Friedmann equations are then given by

$$H^2 = \frac{8\pi G}{3c^2} \rho_b, \\ q = -1 - \frac{\dot{H}}{H^2} \approx \frac{1}{2}, \qquad (50)$$

with $\rho_b$ the energy density of baryons. Based on the observed density of baryons, the expansion from the Hubble parameter should be about $15 \text{ km}^{-1} \text{ s}^{-1} \text{ Mpc}^{-1}$ and the deceleration parameter close to $1/2$. As a consequence, this model cannot account for the observed expansion of the Universe close to 67 km/s/Mpc from the latest data of the Planck mission (Planck Collaboration VI 2020), and the acceleration of the expansion measured by type Ia supernovae (Riess et al. 1998), $q = -1.0 \pm 0.4$.

To account for these discrepancies, the standard cosmological model must invoke the presence of cold dark matter (CDM) and dark energy in the form of a cosmological constant $\Lambda$, leading to the equations

$$H^2 = \frac{8\pi G}{3c^2}(\rho_b + \rho_{\text{CDM}}) + \frac{\Lambda c^2}{3}, \\ q = -1 - \frac{\dot{H}}{H^2} \approx \frac{1}{2} - \frac{\Lambda c^2}{2H}. \qquad (51)$$

This adjustment leads to the classical energy budget of the Universe in which baryonic matter accounts for ~5%, cold dark matter for ~20%, and dark energy ($\Lambda$) for ~75%.

We now show how correlations can be included in the Friedmann equations. As mentioned above, we present the approach within the Newtonian limit to illustrate our point here. Following Peacock (1999), we assume an ensemble of particles of mass $m$ in a volume $V$ of radius $R$, whose fluid density is

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

<figure><img src="image_9.png" /><figcaption>Fig. 7. Radial distribution function of stars and gas corresponding to the parameters given in Table 1. Case 1 corresponds to a fit of the correlation functions of the large-scale structures of the Universe at redshift 0 in the IllustrisTNG simulation from Springel et al. (2017). Case 2 uses the same stellar component, but more compact substructures in the gas component.</figcaption></figure>

assumed to be homogeneous and isotropic at very large scales, but can be correlated at smaller scales, that is $P_1(\mathbf{r}) = 1/V$ but $P_2(\mathbf{r}, \mathbf{r}') \neq P_1(\mathbf{r})P_1(\mathbf{r}')$. The total mechanical energy of this system under its $N$-body or fluid form is given by (see Sect. 2.2)

$$E_m = \langle mv^2/2 \rangle + \langle H_{\text{int}} \rangle. \qquad (52)$$

Similarly to $\alpha_{ih}$ in Sect. 4.1, we define $\alpha_{ni}$ by

$$\alpha_{ni} = \frac{\iint_{V,V} \phi(|\mathbf{r} - \mathbf{r}'|)P_2(\mathbf{r}, \mathbf{r}')dVdV'}{\iint_{V,V} \phi(|\mathbf{r} - \mathbf{r}'|)P_1(\mathbf{r})P_1(\mathbf{r}')dVdV'}. \qquad (53)$$

The ideal (uncorrelated) fluid hypothesis corresponds to $\alpha_{ni} = 1$. Equation (52) can then be rewritten as

$$E_m = \langle mv^2/2 \rangle + \alpha_{ni} \langle H_{\text{int,ideal}} \rangle. \qquad (54)$$

Following Peacock (1999), we take $\langle mv^2/2 \rangle = M\dot{R}^2/2$ and $\langle H_{\text{int,ideal}} \rangle = -GM^2/R$ and rewrite this equation as

$$E_m/M = \frac{\dot{R}^2}{2} - \frac{4\pi}{3}\alpha_{ni}G\rho_m R^2. \qquad (55)$$

Rearranging the different terms, we obtain

$$\frac{\dot{R}^2}{R^2} - \frac{2E_m/M}{R^2} = \frac{8\pi G}{3c^2} \alpha_{ni} \rho_m c^2. \qquad (56)$$
In the uncorrelated case, $\alpha_{ni} = 1$, and with the substitutions $R \to a$ and $-E_m/M \to K$, we recognize the first Friedmann equation in Eq. (50). Therefore, Eq. (56) indicates that correlations should be accounted for through a multiplicative factor $\alpha_{ni}$ of the energy density for a proper nonideal generalization of the Friedmann equation (this could also be seen as a multiplicative factor of the gravitational constant). The last step in obtaining a nonideal first Friedmann equation is to take the thermodynamic limit in $\alpha_{ni}$, $(N, V) \to \infty$ keeping $N/V = \rho$ constant. This can be done by introducing a near-field approximation on a scale $\lambda_H$, that is, by replacing $\phi$ by $\exp(-r/\lambda_H)\phi$ in Eq. (53), or equivalently, introducing a cutoff of the integrals at a radius $r_{bound}$ of approximately $\lambda_H$. In the latter case, and assuming a large-scale flat, homogeneous, and isotropic Universe, we obtain the following first nonideal Friedmann equation:

$$H_{\text{ni}}^2 = \frac{8\pi G}{3c^2} \alpha_{\text{ni}} \rho_b \qquad (57)$$

with

$$\alpha_{\text{ni}} = \frac{\int_{0}^{r_{\text{bound}}} g(r) r dr}{\int_{0}^{r_{\text{bound}}} r dr}. \qquad (58)$$

We can now study the acceleration of the expansion. Because we can link the value of the Hubble parameter to the nonideal amplification induced by bound substructures, $\alpha_{ni}$, it is natural to expect an acceleration of the expansion linked to an increase in the densities of these structures, due to the ongoing gravitational collapse. In order to derive the second nonideal Friedmann equation, we use the first law of thermodynamics in an expanding universe,

$$\frac{d(\rho_e a^3)}{dt} = -P_{\text{ni}} \frac{da^3}{dt}, \qquad (59)$$

where P is a nonideal pressure. This yields the relation $\dot{\rho}_e + 3a/\rho_e = -3\dot{P}/a$. Using this relation and the derivation of the first nonideal Friedmann relation, we obtain the nonideal second Friedmann equation,

$$\dot{H}_{\text{ni}} + \frac{3}{2} H_{\text{ni}}^2 + \frac{K}{2a^2} = \frac{-4\pi G}{c^2} P_{\text{ni}} \alpha_{\text{ni}} + \frac{4\pi G}{3c^2 H_{\text{ni}}} \rho_e \dot{\alpha}_{\text{ni}}. \qquad (60)$$

Assuming a flat geometry, $K = 0$, we obtain the deceleration parameter as

$$q_{\text{ni}} = \frac{1}{2} + \frac{4\pi G}{c^2 H_{\text{ni}}} P_{\text{ni}} \alpha_{\text{ni}} - \frac{\dot{\alpha}_{\text{ni}}}{2H_{\text{ni}}\alpha_{\text{ni}}}. \qquad (61)$$

The nonideal pressure term is negligible as in the ideal case. However, we can obtain a negative deceleration parameter with the third term because of $\dot{\alpha}_{ni}$, which is positive because the densities of substructures increase with time due to the ongoing gravitational collapse. This term can then replace the contribution of the cosmological constant that needs to be introduced in the Friedmann equations to explain the acceleration of the expansion in the standard ΛCDM approach.

To summarize, the two nonideal Friedmann equations (for $K = 0$) in our formalism are given by

$$H_{\text{ni}}^2 = \frac{8\pi G}{3c^2} \alpha_{\text{ni}} \rho_b, \\ q_{\text{ni}} = -1 - \frac{\dot{H}_{\text{ni}}}{H_{\text{ni}}^2} \approx \frac{1}{2} - \frac{\dot{\alpha}_{\text{ni}}}{2H_{\text{ni}}\alpha_{\text{ni}}}, \qquad (62)$$

with

$$\alpha_{\text{ni}} = \frac{\int_{0}^{r_{\text{bound}}} g(r) r dr}{\int_{0}^{r_{\text{bound}}} r dr}. \qquad (63)$$

We can now estimate $g(r)$ from available numerical simulations. We first decompose $\alpha_{ni}$ into the dominant components of visible matter, namely gas and stars,

$$\alpha_{\text{ni}} = X_{\text{star}}^2 \alpha_{\text{ss}} + (1 - X_{\text{star}})^2 \alpha_{\text{gg}} + 2X_{\text{star}}(1 - X_{\text{star}})\alpha_{\text{gs}}, \qquad (64)$$

where ss, gg, and sg denote the star-star, gas-gas, and star-gas interaction contributions, respectively, and $X_{\text{star}}$ denotes the mass fraction in stars relative to the total visible baryonic matter in gas and stars ($X_{\text{star}} \approx 5\%$). We can parameterize the radial distribution functions for gas and stars from the auto- and cross-correlation functions of the large-scale structures of the Universe (Peebles 1980) by using simulations at redshift $z = 0$ from (Springel et al. 2017)<fnref n="1" />, as illustrated in Fig. 7. It is important to stress that the decomposition $\alpha_{ni}$ must be made with squared mass fractions and account for cross-correlations because the interaction energy depends formally on $P_2(r, r')$, which can be seen as the square of the density (similarly to Eq. (15) in Springel et al. 2017). We recall the link between the correlation function and the radial distribution function $\xi(r) = g(r) - 1$. The radial distribution functions can be parameterized with the simple forms

$$g_{ij}(r) = \max \left( g_{1\text{Mpc}} \times \left( \frac{r}{1\text{Mpc}} \right)^{-\beta_0}, 1 \right), \quad r < 1\text{Mpc} \\ = \max \left( g_{1\text{Mpc}} \times \left( \frac{r}{1\text{Mpc}} \right)^{-\beta_1}, 1 \right), \quad r \ge 1\text{Mpc} \qquad (65)$$

We assume for simplicity that the radial distribution function is equal to 1 at large distances, while strictly speaking, it should be slightly smaller than 1 and tend to 1 at infinity. The parameters $g_{1\text{Mpc}}$, $\beta_0$, and $\beta_1$ are given in Table 1 for two cases. Case 1 corresponds to the fit of the correlation functions presented in Springel et al. (2017), assuming $\xi(r) \approx g(r)$ for $\xi(r) \gg 1$. Case 2 has the same parameterization of the stellar component as case 1, is well constrained by observations (e.g., Li & White 2009), but has more compact substructures in the gas component. This modification is ad hoc for the moment and further investigations are needed to confirm if its amplitude is compatible with the observational constraints that are available for the gas component. The main point of case 2, however, is to explore the sensitivity of the nonideal amplification to the gas substructures.

We then define the bound radius $r_{bound}$ as the minimum radius, such that $g_{ij}(r) < g_{th}$ for $r > r_{bound}$, and we explore two values for the threshold value: $g_{th} = 1$ and $g_{th} = 4$.

Table 2 gives the obtained nonideal amplification factors for the different components for the two threshold values, as well as the total nonideal value $\alpha_{ni}$. We also give the corresponding nonideal Hubble parameter for star mass fractions $X_{\text{star}} = 5\%$ and $X_{\text{star}} = 10\%$, respectively. In the most conservative case (case 1 with $g_{th} = 1$), we obtain an overall amplification factor between 5 and 6. With a modification of the gas substructures (case 2 with $g_{th} = 1$) or a different threshold (case 1 with $g_{th} = 4$),

<footnote n="1">¹ As the structures obtained in ΛCDM cosmological simulations represent the observations well, we use them as the reference in our calculations. In contrast to the conventional ΛCDM model, however, dark matter and dark energy in our formalism are in reality proxies for non-ideal effects induced by substructures.</footnote>
<table><caption>**Table 1.** Parameters used to parameterize the radial distribution function of gas and stars for case 1 and case 2.</caption><thead><tr><th colspan="4">Case 1</th></tr><tr><th></th><th>g<sub>1Mpc</sub></th><th>β<sub>0</sub></th><th>β<sub>1</sub></th></tr></thead><tbody><tr><th>g<sub>gg</sub></th><td>50</td><td>0.75</td><td>1.9</td></tr><tr><th>g<sub>gs</sub></th><td>120</td><td>1.15</td><td>1.9</td></tr><tr><th>g<sub>ss</sub></th><td>200</td><td>2.40</td><td>1.9</td></tr></tbody><thead><tr><th colspan="4">Case 2</th></tr><tr><th></th><th>g<sub>1Mpc</sub></th><th>β<sub>0</sub></th><th>β<sub>1</sub></th></tr></thead><tbody><tr><th>g<sub>gg</sub></th><td>50</td><td>1.2</td><td>3.2</td></tr><tr><th>g<sub>gs</sub></th><td>120</td><td>1.5</td><td>3.2</td></tr><tr><th>g<sub>ss</sub></th><td>200</td><td>2.4</td><td>1.9</td></tr></tbody></table>

**Notes.** The corresponding radial distribution functions are plotted in Fig. 7.

<table><caption>**Table 2.** Nonideal amplification of the different component (gas-gas, star-star, and gas-star) for two different values of the threshold $g_{th}$ used to define the bound radius.</caption><thead><tr><th colspan="5">Case 1</th></tr></thead><tbody><tr><th>g<sub>th</sub></th><td>1</td><td>4</td><td></td><td></td></tr><tr><th>α<sub>gg</sub></th><td>5.02</td><td>15.56</td><td></td><td></td></tr><tr><th>α<sub>gs</sub></th><td>6.28</td><td>21.0</td><td></td><td></td></tr><tr><th>α<sub>ss</sub></th><td>61.0</td><td>257</td><td></td><td></td></tr><tr><th>X<sub>star</sub></th><td>0.05</td><td>0.1</td><td>0.05</td><td>0.1</td></tr><tr><th>α<sub>ni</sub></th><td>5.28</td><td>5.80</td><td>16.7</td><td>18.9</td></tr><tr><th>H<sub>ni</sub></th><td>34.2</td><td>35.8</td><td>60.7</td><td>64.7</td></tr><tr><th>α&#x0307;<sub>ni</sub></th><td>0.12</td><td>0.13</td><td>0.20</td><td>0.21</td></tr><tr><th>q<sub>ni</sub></th><td>-1.20</td><td>-1.21</td><td>-1.12</td><td>-1.09</td></tr></tbody><thead><tr><th colspan="5">Case 2</th></tr></thead><tbody><tr><th>g<sub>th</sub></th><td>1</td><td>4</td><td></td><td></td></tr><tr><th>α<sub>gg</sub></th><td>16.4</td><td>36.2</td><td></td><td></td></tr><tr><th>α<sub>gs</sub></th><td>31.7</td><td>72.7</td><td></td><td></td></tr><tr><th>α<sub>ss</sub></th><td>61.0</td><td>257</td><td></td><td></td></tr><tr><th>X<sub>star</sub></th><td>0.05</td><td>0.1</td><td>0.05</td><td>0.1</td></tr><tr><th>α<sub>ni</sub></th><td>17.9</td><td>19.6</td><td>40.2</td><td>45.0</td></tr><tr><th>H<sub>ni</sub></th><td>62.9</td><td>65.8</td><td>94.3</td><td>94.3</td></tr><tr><th>α&#x0307;<sub>ni</sub></th><td>0.21</td><td>0.21</td><td>0.26</td><td>0.27</td></tr><tr><th>q<sub>ni</sub></th><td>-1.09</td><td>-1.09</td><td>-0.87</td><td>-0.83</td></tr></tbody></table>

**Notes.** Total nonideal amplification $\alpha_{ni}$ and the corresponding nonideal Hubble parameter ($H_{ni}$ in units of km s$^{-1}$ Mpc$^{-1}$) for two different values of the star mass fraction relative to the total baryon mass in the Universe (we assume a present-day baryon density of 0.25 particle/m$^3$). Estimation of $\dot{\alpha}_{ni}$ assuming a 14 Gyr old Universe and the corresponding nonideal deceleration parameter $q_{ni}$.

we obtain a total amplification of about 20, which yields the observed value of the Hubble parameter (about 67 km/s/Mpc, Planck Collaboration VI 2020). The last case is extreme (case 2 with $g_{th} = 4$) and shows that the nonideal model can give a very rapid expansion, depending on the contribution of substructures to the gravitational energy.

It is also interesting to consider the nonideal amplification per component: because the stars are far more substructured (correlated) than the gas, the nonideal amplification for stars is between 5 and 20 times higher than the amplification for the gas. This effect could thus explain observations like the bullet cluster (Markevitch et al. 2004): when stars and gas are spatially separated, the nonideal amplification in the stellar component

is significantly higher than the amplification in the gas. If the difference is about 20 or more, the nonideal amplification can compensate the difference in mass fractions, potentially explaining the observed high gravitational weak lensing related to the stellar component.

We also give in Table 2 estimates for $\dot{\alpha}_{ni}$ and the corresponding deceleration parameter $q_{ni}$, assuming $\dot{\alpha}_{ni}/\alpha_{ni} \approx \ln(\alpha_{ni})/t_u$, with $t_u \approx 14$ Gyr the age of the Universe. The deceleration parameter is always close to the observed value $q = -1.0 \pm 0.4$ (Riess et al. 1998). An interesting feature of this nonideal model is that in contrast to a $\Lambda$CDM approach that needs to introduce two parameters $\rho_{CDM}$ and $\Lambda$ to explain the expansion and its acceleration, we only introduce one parameter $\alpha_{ni}$. We therefore link the expansion to its acceleration and provide an analytical expression for $q_{ni}$, using the critical density $\rho_c = 3H^2/8\pi G$. Assuming the nonideal amplification factor $\alpha_{ni}$ entirely explains the observed present-day expansion, that is, it varied from $\alpha_{ni} = 1$ to $\alpha_{ni} = \rho_c/\rho_b \approx 20$ over the age of the universe, $t_u \sim 14$ Gyr, we can obtain an order-of-magnitude analytical expression,

$$q_{ni} \approx \frac{1}{2} - \frac{\ln(\rho_c/\rho_b)}{2t_u H}, \quad (66)$$

which gives a deceleration parameter $q_{ni} = -1.06$ that is compatible with the value based on type Ia supernovae ($q = -1.0 \pm 0.4$, Riess et al. 1998).

This value contradicts current estimates using all cosmological constraints ($q = -0.52 \pm 0.07$, Planck Collaboration VI 2020). It must be stressed, however, that the aforementioned estimates of the deceleration parameter based on the Planck results rely on the $\Lambda$CDM model. The current nonideal cosmological formalism might relax this tension: an older universe with an age of about 20 Gyr would give a nonideal deceleration parameter of about -0.5, for example. Another possibility is that the formation of structures has slowed down in the recent universe, for instance, because of feedback processes. Observations show that the cosmic star formation rate has decreased by almost a factor 10 since redshift $z \sim 1-2$, that is, within the past ~8 - 10 Gyr (Madau & Dickinson 2014). An increase in nonideal amplification factor of about 2 within the last 7 Gyr would agree with current estimates of $q$. This possibility can be explored in detail with numerical simulations. Nonetheless, because the physics responsible for the acceleration of the Universe remains mysterious so far, it is quite remarkable that we can derive a good order-of-magnitude estimate of the acceleration of the expansion based on a physically intuitive and appealing argument, namely the increasing density induced by the gravitational collapse of the large-scale structures.

# 5. Discussion and conclusion

## 5.1. Nonideal Einstein equation?

As discussed in Sects. 2.2 and 2.3, we may have to step away from the traditional forms of the Poisson and Einstein equations to obtain the nonideal Friedmann equations. For a nonrelativistic fluid, the mass energy density, $E_m = \rho c^2$, is dominant in the zeroth component of the stress energy tensor, $T_{00} \approx \rho c^2$, and it would also dominate all other forms of energy, kinetic and potential: $E_m \gg E_{\text{kin}}$, $E_p$, and so on. At first sight, it seems impossible for an interaction potential energy $E_p$ to dominate the mass energy, implying that correlations in the interaction energy should not modify the Einstein equations. It must be stressed, however, that in the nonrelativistic limit of the Einstein equations, the gravitational interaction energy has a special status and should not
be identified with $E_p$, but rather directly with the integral of $T_{00}$. This can be seen by deriving the Poisson equation from the Einstein equations: $\Delta\phi \approx -4\pi GT_{00}/c^2$. The gravitational interaction energy in the nonrelativistic limit is thus given by

$$ H_{\text{int}} = \frac{1}{2} \int_V \rho(\mathbf{r})\phi(\mathbf{r})dV \\ \approx -\frac{1}{2} \iint_{V,V'} G \frac{\rho(\mathbf{r})T_{00}(\mathbf{r}')/c^2}{\|\mathbf{r} - \mathbf{r}'\|} dVdV', \quad (67) $$

which corresponds to the uncorrelated (ideal) version of the interaction energy. It ignores the distance between pairs of particles in the fluid because it does not depend on the correlation function. Heuristically, this explains why the contribution to the gravitational interaction energy arising from correlations should be introduced as a multiplicative factor of the mass energy (or of the gravitational constant) and not as an additional term.

Another way to show the limits of using the Einstein equation for self-gravitating correlated fluids is to consider the relativistic version of the nonideal Navier-Stokes equations introduced in Sect. 3.2: $\nabla_\mu T^{\mu\nu} = 0$. As explained in Sect. 3.2, as soon as the Jeans length is not resolved (hence when the cosmological principle is applied at large scale), part of the gravitational interactions should be taken into account in the EOS of the fluid, hence in the definition of $T^{\mu\nu}$. The rest can be accounted for as external interactions with a mean field that can be included in the covariant derivative $\nabla_\mu$ in the context of General Relativity. The gravitational interactions that contribute to the EOS cannot in general be included in the covariant derivative, which is another way to show that the equivalence principle for all the gravitational interactions is broken by the presence of correlations in a self-gravitating nonideal fluid.

We stress that the concept of nonideal self-gravity developed in the present paper does not modify the fundamental law of gravity between two point particles, which relies only on Newton's law. Furthermore, it does not contradict Poisson and Einstein theories: they always remain valid if one can have access to or compute exactly the full inhomogeneous density distributions down to the scale below which we can assume the gas to be uncorrelated (see the example of polytropic stars in Sect. 4.1). However, for a large ensemble of particles with unresolved small-scale inhomogeneities, a nonideal modification of the Poisson and Einstein equations is required to account for the correlations between these substructures. It is worth pointing out that a modification of the Poisson equation because of correlations is a procedure already used in the case of the electrostatic force in the context of physical chemistry (see, e.g., Storey & Bazant 2012). In this field, Poisson and Einstein theories would be qualified as mean field theories, that is, they would only be valid in the uncorrelated case. A more complete theoretical approach to develop a nonideal formalism for the Einstein equations could come from two routes. Either from a homogenization procedure of the inhomogeneous Einstein equations, a path that is already currently explored (see, e.g., Buchert 2000, 2008), or from an extension of the concepts used in statistical mechanics, as developed in this paper, within a general relativistic framework.

## 5.2. Conclusion

Starting from the well-established BBGKY formalism, we have shown that correlations induced by bound substructures should be accounted for when describing the energy-mass density of the Universe. This should be done via the use of a nonideal virial

theorem formalism and nonideal Navier-Stokes equations. When these nonideal effects are taken into account, an amplification of the gravitational interaction energy results that could account, at least partly, for the missing mass problem in galaxies, clusters of galaxies, and large-scale structures in general. The strength of this model is that the radial distribution function (or the two-point correlation function) can be well determined in galaxies and the large-scale structures either by observations or by simulations.

By examining the viscous limit of the nonideal Navier-Stokes equations induced by the presence of small-scale interactions, we showed that the presence of substructures (e.g., spiral arms and local clustering) can produce nonkeplerian rotation profiles that could explain the observed flat rotation curves. A consequence of this approach is to show that the observed bifurcation between spiral and dwarf galaxies and some LSB galaxies with Keplerian rotation profiles could be explained by the presence or absence of substructures within these galaxies. In the context of statistical mechanics, this bifurcation could be interpreted as a phase transition between two different equilibrium (viscous and nonviscous) states.

Using this nonideal virial theorem, we have derived nonideal Friedmann equations within the Newtonian limit. Using the radial distribution function of visible gas and stars obtained in large-scale simulations, we can compute the nonideal amplification factor of gravitational mass energy and show that this factor can easily be about 5 to 20 and thus could explain the observed value of the Hubble parameter. Furthermore, we show that the amplification is much stronger for the stellar component than for the gas component, which could also explain the bullet cluster observation with a high gravitational weak lensing in the stellar component compared to the gas component.

Furthermore, because most of the contribution to the value of the Hubble parameter stems from the nonideal amplification caused by the interactions between substructures, this naturally explains the acceleration of the expansion of the Universe as a consequence of the ongoing collapse, thus the increasing densities of these substructures. An amplification factor $\alpha_{ni} \simeq 20$ during most of the lifetime of the Universe yields a deceleration parameter $q_{ni} \simeq -1$ (see Eq. (66)), close to the estimates based on type Ia supernovae. To the best of our knowledge, this model is the first that can predict a coherent value for the acceleration of the expansion of the Universe with first-principle physical arguments.

Observations of the *Euclid* satellite will be crucial to probe and further constrain the radial distribution function and the nonideal effects of self-gravitating matter at large scales. A precise characterization of these effects will be essential to build a full nonideal cosmological model addressing all cosmological constraints (e.g., cosmic microwave background data) and to eventually reveal what is really dark in our Universe.

*Acknowledgements.* We thank the anonymous referee for his/her constructive comments. We thank S. Kokh and E. Audit for co-supervising the PhD of T. Padioleau. This PhD subject is at the frontier between astrophysics and two-phase flow simulations for the cooling system of nuclear power plants and the link between liquid statistical mechanics and gravity at large scale has emerged from this trans-disciplinary approach. They also thank CEA in general (DRF and DES in particular) for the funding of the PhD and for providing a favorable environment for trans-disciplinarity at Maison de la Simulation. P.T. also thanks the Astrosim 2019 school and in particular J. Rosdhal and O. Hahn for their lectures: part of the ideas presented in this paper has emerged after the school. The authors also thank F. Sainsbury-Martinez, J. Faure, T. Buchert, G. Laibe, Q. Vigneron, P. Hennebelle, F. Bournaud, D. Elbaz, D. Borgis, R. Balian, and L. Saint-Raymond for interesting discussions and useful comments on this paper. P.T. acknowledges supports by the European Research Council under Grant Agreement ATMO 757858.
References

Aslangul, C. 2006, *Physique statistique des Fluides Classiques*

Balian, R., & Schaeffer, R. 1989a, A&A, 220, 1

Balian, R., & Schaeffer, R. 1989b, A&A, 226, 373

Barnes, J. E., & Hernquist, L. 1992, ARA&A, 30, 705

Binney, J., & Tremaine, S. 1987, Galactic Dynamics

Bogoliubov, N. N. 1946, J. Phys.-USSR, 10, 265

Born, M., & Green, H. S. 1946, Proc. Roy. Soc. London Ser. A Math. Phys. Sci., 188, 10

Buchert, T. 2000, in *Proceedings of 9th Workshop on General Relativity and Gravitation conference*, 306

Buchert, T. 2008, Gen. Relativ. Gravit., 40, 467

Burkert, A., & Hensler, G. 1988, A&A, 199, 131

Chabrier, G. 1990, J. Phys. France, 51, 1607

Chavanis, P. H. 2013, A&A, 556, A93

Davis, M., & Peebles, P. J. E. 1977, ApJS, 34, 425

de Vega, H., & Sánchez, N. 2002, Nucl. Phys. B, 625, 409

Hansen, J.-P., & McDonald, I. R. 2006, Theory of Simple Liquids, 3rd edn., eds. J.-P. Hansen, & I. R. McDonald (Burlington: Academic Press), 11

Harsoula, M., & Kalapotharakos, C. 2009, MNRAS, 394, 1605

Irving, J. H., & Kirkwood, J. G. 1950, J. Chem. Phys., 18, 817

Kirkwood, J. G. 1946, J. Chem. Phys., 14, 180

Kirkwood, J. G., Buff, F. P., & Green, M. S. 1949, J. Chem. Phys., 17, 988

Li, C., & White, S. D. M. 2009, MNRAS, 398, 2177

Madau, P., & Dickinson, M. 2014, ARA&A, 52, 415

Markevitch, M., Gonzalez, A. H., Clowe, D., et al. 2004, ApJ, 606, 819

Paolo, C. D., Salucci, P., & Fontaine, J. P. 2019, ApJ, 873, 106

Peacock, J. 1999, *Cosmological Physics*, Cambridge Astrophysics (Cambridge University Press)

Peebles, P. J. E. 1980, *The Large-Scale Structure of the Universe*

Pfenniger, D. 2006, Comp. Ren. Phys., 7, 360

Planck Collaboration VI. 2020, A&A, 641, A6

Riess, A. G., Filippenko, A. V., Challis, P., et al. 1998, AJ, 116, 1009

Springel, V., Pakmor, R., Pillepich, A., et al. 2017, MNRAS, 475, 676

Storey, B. D., & Bazant, M. Z. 2012, Phys. Rev. E, 86, 056303

Truelove, J. K., Klein, R. I., McKee, C. F., et al. 1997, ApJ, 489, L179

Voglis, N., Stavropoulos, I., & Kalapotharakos, C. 2006, MNRAS, 372, 901

Whitaker, S. 1998, *The Method of Volume Averaging, Theory and Applications of Transport in Porous Media* (Netherlands: Springer)

Yvon, J. 1935, *La théorie statistique des fluides et l'équation d'état, Actualités Scientifiques et Industrielles: Théories Mécaniques* (Hermann & cie)