text:
# A Model For Infall Around Virialized Halos

## 1 Introduction

Numerical simulations of hierarchical halo formation indicate a roughly universal spherically-averaged density profile within virialized dark matter halos (Navarro, Frenk, & White, 1997). Such profiles are generated by complicated processes that involve violent relaxation and multiple shell crossing; the great interest in them results from the possibility of a direct comparison to the profiles inferred from galactic rotation curves (e.g., Salucci & Burkert, 2000; de Blok, McGaugh, & Rubin, 2001; van den Bosch & Swaters, 2001) and from lensing in clusters (e.g., Sand, Treu, & Ellis, 2002; , 2002). However, the analogous question of what is the typical density profile outside the virial radius has received relatively little attention, because of the difficulty of directly observing density profiles even out to the halo virial radius; this difficulty results from the rapid decrease in the gas and dark matter densities with distance from the halo center. However, the theoretical description is simpler for infall outside the virialization radius, since in cold dark matter (hereafter CDM) models there is no shell crossing until the final non-linear collapse.

The gas around halos can only be detected using physical effects that are very sensitive to low-density gas. One such phenomenon is resonant Ly$`\alpha `$ absorption; Gunn & Peterson (1965) first noted that even an H I neutral fraction $`10^5`$ suffices to make gas at the cosmic mean density strongly absorb photons at the Ly$`\alpha `$ resonance. Thus, the pattern of gas infall around halos may affect observable properties of the Ly$`\alpha `$ absorption in these regions. At moderate redshifts, these regions near halos play a prominent role in measurements of the proximity effect of quasars (e.g., Scott et al., 2000; Bajtlik, Duncan, & Ostriker, 1988). In these measurements, the effect of the high quasar intensity on the statistics of Ly$`\alpha `$ absorption is analyzed and used to infer the mean ionizing background, since the quasar intensity dominates only as long as it is significantly stronger than the background. A different application may be possible at very high redshifts, corresponding to the time before the universe was fully reionized; the presence of the fully neutral IGM could be detected through absorption of a bright source by the Ly$`\alpha `$ damping wing of the IGM (Miralda-Escudé, 1998). However, every bright source creates a surrounding H II region which complicated the interpretation (Cen & Haiman, 2000; Madau & Rees, 2000). The properties of the gas in the regions that lie just outside the host halo of the source galaxy or quasar should play a major role in these various applications of the observed absorption.

The various calculations of absorption cited above neglected the enhanced densities and the negative (i.e., infall) velocities expected for the gas surrounding a massive halo at high redshift. This infall, however, may be an essential component of any correct interpretation of observations; Loeb & Eisenstein (1995) applied it to the proximity effect and showed that neglecting infall can lead to an overestimate of the ionizing background flux by up to a factor of three. Applying the same infall profile, Barkana & Loeb (2002) have recently revealed the signature of infalling gas around high redshift quasars. The spectral pattern due to Ly$`\alpha `$ absorption by the infalling gas can be used to estimate the quasar’s ionizing intensity and the total mass of its host halo. While the spectra of individual objects are marked by the large density fluctuations expected in the intergalactic medium (hereafter IGM), if this pattern is averaged over many objects then the resulting mean profile could be compared directly with an accurate theoretical prediction for the density and velocity profile of the infalling gas.

Loeb & Eisenstein (1995) defined a halo in the same way as halos were defined in the classic model of Press & Schechter (1974). In this approach, halos on a given scale are defined as forming when the initial density averaged on this scale is higher than a threshold that depends on the collapse redshift and is fixed using spherical collapse. The halo mass function derived from this assumption, when multiplied by an ad-hoc correction factor of two, provides a reasonable match to the results of numerical simulations (e.g., Katz, Quinn, & Gelb, 1993).

The same mass function was later derived by Bond et al. (1991) using a self-consistent approach, with no need for external factors. In their derivation, the factor of two has a more satisfactory origin, namely the so-called “cloud-in-cloud” problem: if the average initial density on a certain scale is above the collapse threshold, this region may be contained within a larger region that is itself also denser than the collapse threshold; in this case the original region should be counted as belonging to the halo with mass corresponding to the larger collapsed region. Bond et al. (1991) and Press & Schechter (1974) both give the same halo mass function and thus appear to be at least statistically valid given the agreement with the mass function seen in numerical simulations. However, the Bond et al. (1991) model is more satisfactory since it accounts for the cloud-in-cloud problem and, more importantly, allows for calculations of other halo properties in addition to simple mass functions, all in an approach that is more self-consistent.

The Bond et al. (1991) approach (also known as extended Press-Schechter) has been used to calculate halo merger histories (Lacey & Cole, 1993), halo mass functions in models with warm dark matter (Barkana, Haiman, & Ostriker, 2001), halo mass functions in a model with ellipsoidal collapse (Sheth, Mo, & Tormen, 2001), and the nonlinear biasing of halos (Porciani et al., 1998; Scannapieco & Barkana, 2002). In this paper, we present a new application of this approach in which we obtain the expected profile of infalling matter around virialized halos. The rest of this paper is organized as follows. In § 2.1 we establish our notation and review the Bond et al. (1991) derivation of the Press & Schechter (1974) halo mass function. In § 2.2 we review the infall profile used by Loeb & Eisenstein (1995) and also derive an alternative model based on similar assumptions but carried out in Fourier space. In § 2.3 we derive our model based on the extended Press-Schechter formalism. In § 2.4 we then discuss the modified picture of halo virialization once infall is included. In § 3 we illustrate the predictions of our model for the initial and final density profiles surrounding halos. Finally in § 4 we give our conclusions.

## 2 A Model for Infall: Derivation

### 2.1 Review of Halo Collapse

Before addressing the problem of the initial density profile around a virialized halo, we briefly review in this section the approach of Bond et al. (1991) which leads to the halo mass function. We work with the linear overdensity field $`\delta (𝐱,z)\rho (𝐱,z)/\overline{\rho }(z)1`$, where $`𝐱`$ is a comoving position in space, $`z`$ is the cosmological redshift and $`\rho `$ is the mass density, with $`\overline{\rho }`$ being the cosmic mean density. In the linear regime, the density field maintains its shape in comoving coordinates and the overdensity simply grows as $`\delta =\delta _iD(z)/D(z_i)`$, where $`z_i`$ and $`\delta _i`$ are the initial redshift and overdensity, and $`D(z)`$ is the linear growth factor (Peebles, 1980). When the overdensity in a given region becomes non-linear, the expansion halts and the region turns around and collapses to form a virialized halo.

The time at which the region virializes can be estimated based on the initial linear overdensity, using as a guide the collapse of a spherical tophat perturbation. At the moment at which a tophat collapses to a point, the overdensity predicted by linear theory is $`\delta _c=1.686`$ (Peebles, 1980) in the Einstein-de Sitter model. This value depends slightly on cosmological parameters, and even in the Einstein-de Sitter case it is modified by infall, as we discuss in §2.4; our analytical expressions in this section are valid for any value of $`\delta _c`$.

A useful alternative way to view the evolution of density is to consider the linear density field extrapolated to the present time, i.e., the initial density field at high redshift extrapolated to the present by multiplication by the relative growth factor. In this case, the critical threshold for collapse at redshift $`z`$ becomes redshift dependent,

$$\delta _c(z)=\delta _c/D(z).$$
(1)
We adopt this view, and throughout this paper the power spectrum $`P(k)`$ refers to the initial power spectrum, linearly extrapolated to the present (i.e., not including non-linear evolution).

At a given $`z`$, we consider the smoothed density in a region around a fixed point $`A`$ in space. We begin by averaging over a large mass scale $`M`$, and then lower $`M`$ until we find the highest value for which the averaged overdensity is higher than $`\delta _c(z)`$; we assume that the point $`A`$ belongs to a halo with a mass $`M`$ corresponding to this filter scale.

In this picture we can derive the mass distribution of halos at a redshift $`z`$ by considering the statistics of the smoothed linear density field. If the initial density field is a Gaussian random field and the smoothing is done using sharp $`k`$-space filters, then the value of the smoothed $`\delta `$ undergoes a random walk as the cutoff value of $`k`$ is increased. If the random walk first hits the collapse threshold $`\delta _c(z)`$ at $`k`$, then at a redshift $`z`$ the point $`A`$ is assumed to belong to a halo with a mass corresponding to this value of $`k`$. Instead of using $`k`$ or the halo mass, we adopt as the independent variable the variance at a particular filter scale $`k`$,

$$S_k\frac{1}{2\pi ^2}_0^k𝑑k^{}k^2P(k^{}).$$
(2)
In order to construct the number density of halos in this approach, we need to solve for the evolution of the probability distribution $`Q(\delta ,S_k)`$, where $`Q(\delta ,S_k)d\delta `$ is the probability for a given random walk to be in the interval $`\delta `$ to $`\delta +d\delta `$ at $`S_k`$. Alternatively, $`Q(\delta ,S_k)d\delta `$ can also be viewed as the trajectory density, i.e., the fraction of the trajectories that are in the interval $`\delta `$ to $`\delta +d\delta `$ at $`S_k`$, assuming that we consider a large ensemble of random walks all of which begin with $`\delta =0`$ at $`S_k=0`$.

Bond et al. (1991) showed that $`Q`$ satisfies a diffusion equation,

$$\frac{Q}{S_k}=\frac{1}{2}\frac{^2Q}{\delta ^2},$$
(3)
which is solved by the Gaussian solution which we label $`Q_0`$:

$$Q_0(\delta ,S_k)=\frac{1}{\sqrt{2\pi S_k}}\mathrm{exp}\left[\frac{\delta ^2}{2S_k}\right].$$
(4)
To determine the probability of halo collapse at a redshift $`z`$, we consider the same situation but with an absorbing barrier at $`\delta =\nu `$, where we set $`\nu =\delta _c(z)`$. The fraction of trajectories absorbed by the barrier up to $`S_k`$ corresponds to the total mass fraction contained in halos with masses higher than the value associated with $`S_k`$. In this case, $`Q`$ again satisfies eq. (3) and the solution with the barrier in place is given by adding an extra image-solution:

$$Q(\nu ,\delta ,S_k)=Q_0(\delta ,S_k)Q_0(2\nu \delta ,S_k).$$
(5)
Using this expression, the fraction of all trajectories that have passed above the barrier $`\nu `$ by $`S_k`$ is

$$F(\nu ,S_k)=1_{\mathrm{}}^\nu 𝑑\delta Q(\nu ,\delta ,S_k),$$
(6)
and the differential mass function is found \[using eq. (3)\] to be

$$f(\nu ,S_k)=\frac{}{S_k}F(\nu ,S_k)=\frac{\nu }{\sqrt{2\pi }S_k^{3/2}}\mathrm{exp}\left[\frac{\nu ^2}{2S_k}\right].$$
(7)
As $`f(\nu ,S_k)dS_k`$ is the probability that point $`A`$ is in a halo with mass in the range corresponding to $`S_k`$ to $`S_k+dS_k`$, the halo abundance is then simply

$$\frac{dn}{dM}=\frac{\overline{\rho }}{M}\left|\frac{dS_k}{dM}\right|f(\nu ,S_k),$$
(8)
where $`dn`$ is the comoving number density of halos with masses in the range $`M`$ to $`M+dM`$. The cumulative mass fraction in halos above mass $`M`$ is similarly determined to be

$$F(>M|z)=\mathrm{erfc}\left(\frac{\nu }{\sqrt{2S_k}}\right).$$
(9)
While these expressions were derived in reference to density perturbations smoothed by a sharp $`k`$-space filter as given in eq. (2), $`S_k`$ is often replaced in the final results with the variance of the mass $`M`$ enclosed in a spatial sphere of comoving radius $`r`$:

$$\sigma ^2(M)=\sigma ^2(r)=\frac{1}{2\pi ^2}_0^{\mathrm{}}k^2𝑑kP(k)W^2(kr),$$
(10)
where $`W(x)`$ is the spherical tophat window function, defined in Fourier space as

$$W(x)3\left[\frac{\mathrm{sin}(x)}{x^3}\frac{\mathrm{cos}(x)}{x^2}\right],$$
(11)
and $`M=4\pi \overline{\rho }(0)r^3/3`$ in terms of the mean density of matter at $`z=0`$. With these replacements we recover the cumulative mass fraction that was originally derived (Press & Schechter, 1974) simply by considering the distribution of overdensities at a single point, smoothing with a tophat window function, and integrating from $`\delta _c`$ to $`\mathrm{}`$. In this derivation the authors were forced to multiply their result by an arbitrary factor of two, to account for the mass in underdense regions. The excursion-set derivation presented here, based on Bond et al. (1991), properly includes all the mass by accounting for small regions that lie within overdensities on larger scales. This approach also makes explicit the approximations involved in working with $`\sigma ^2(r).`$ Strictly speaking, dealing with a real-space filter requires a complete recalculation of $`f(\nu ,\sigma ^2)`$ which accounts for the correlations intrinsic to $`W(x)`$. However, simply replacing $`S_k`$ with $`\sigma ^2(r)`$ in eq. (8) has been shown to be in reasonable agreement with numerical simulations (e.g., Katz, Quinn, & Gelb, 1993), and is thus a standard approximation.

### 2.2 Infall based on the Press-Schechter Model

In the Press-Schechter model, a halo of mass $`M`$ is assumed to form at redshift $`z`$ if the corresponding comoving sphere (of radius denoted $`r_M`$) has a linear overdensity of $`\delta =\nu \delta _c(z)`$. Loeb & Eisenstein (1995) used this fact to derive a mean infall profile as follows. Consider spheres of various comoving radii $`r`$, all centered on the region in the initial density field that contains the mass that ends up in the final virialized halo. The mean enclosed overdensity $`\delta (r)`$ in each such sphere (linearly extrapolated to the present) is a Gaussian variable with variance $`\sigma ^2(r)`$ as given in eq. (10). The overdensities in two such spheres with radii $`r_1`$ and $`r_2`$ are correlated with a cross-correlation

$$\xi _r(r_1,r_2)\frac{1}{2\pi ^2}_0^{\mathrm{}}k^2𝑑kP(k)W(kr_1)W(kr_2).$$
(12)
Thus, the expected value of $`\delta (r)`$ given that $`\delta _M\delta (r_M)=\nu `$ is given by their joint Gaussian distribution as

$$\frac{\delta (r)_{\mathrm{PS}\mathrm{r}}}{\nu }=\frac{\xi _r(r_M,r)}{\sigma ^2(r_M)},$$
(13)
where the notation $`PSr`$ indicates a derivation based on the Press-Schechter model and carried out in $`r`$-space.

In addition to obtaining the mean profile, it is useful to calculate the scatter in the initial profiles around halos of a given mass in order to assess the applicability of the mean profile to individual objects. If in the future the model is to be compared to halos in numerical simulations or in real observations, a quantitative formula for the scatter helps us determine the number of objects necessary to average over before the mean profile emerges. In addition, the scatter itself is a prediction of the model that in principle can be compared to simulations or observations. In the present model, the ratio $`\delta (r)/\nu `$ follows a Gaussian distribution with the mean value given above and a standard deviation given by

$$\mathrm{\Sigma }(r)_{\mathrm{PS}\mathrm{r}}=\sqrt{\frac{\sigma ^2(r)}{\nu ^2}\frac{\xi _r^2(r_M,r)}{\nu ^2\sigma ^2(r_M)}}.$$
(14)
Alternatively, we can create an infall model based on similar assumptions, but (as in the previous subsection) calculated in $`k`$-space. Consider a halo that forms on a scale with a corresponding variance $`S_{k,M}`$, and consider a larger scale around it with variance $`S_k<S_{k,M}`$. We obtain here a Press-Schechter model in $`k`$-space by considering random trajectories but without including an absorbing barrier; the self-consistent inclusion of a barrier, which accounts for peaks-within-peaks and is the hallmark of the extended Press-Schechter model, is considered in the next subsection.

In the absence of a barrier, the probability distribution of the overdensity $`\delta `$ at any $`S_k`$ is Gaussian, and the values of $`\delta `$ at multiple points are jointly normal. Now, $`\delta `$ has variance $`S_k`$ and $`S_{k,M}`$ at the two scales of interest. Furthermore, since the random walk from $`S_k`$ to $`S_{k,M}`$ is independent of the random walk from 0 to $`S_k`$, the cross-correlation of the values of $`\delta `$ at $`S_k`$ and at $`S_{k,M}`$ is given by the overlapping portion of the two random walks: the cross-correlation is therefore simply $`S_k`$. Similarly to the $`r`$-space case above, the joint Gaussian distribution yields

$$\frac{\delta (r)_{\mathrm{PS}\mathrm{k}}}{\nu }=\frac{S_k}{S_{k,M}}.$$
(15)
As in the PS-r model, the distribution of $`\delta (r)/\nu `$ is Gaussian, where now the standard deviation is

$$\mathrm{\Sigma }(r)_{\mathrm{PS}\mathrm{k}}=\sqrt{\frac{S_k}{\nu ^2}\frac{S_k^2}{\nu ^2S_{k,M}}}.$$
(16)
Note that both of the models developed in this section satisfy the continuity condition $`\delta (r)/\nu 1`$ when $`rr_M`$ (or equivalently $`S_kS_{k,M}`$). However, we show in the next section that they fail to give a physically acceptable result in a different limit, when $`\nu \sqrt{S_{k,M}}`$. Moreover, we have ignored in this section the failure of the simple Press-Schechter approach to account for negative fluctuations, and have not included the factor-of-two normalization fudge factor; in the next section, in contrast, we develop a model that is self-consistent.

### 2.3 Infall based on the Extended Press-Schechter Model

In this section we derive a new model for the initial density profile around a virialized halo. As additional motivation, we first show that the two simple models considered in the previous section yield results that are physically unacceptable, once we consider the obvious possibility that a small halo will lie within a bigger one. Consider, for example, the PS-r model in a CDM cosmology (the cosmological parameters are listed in section 3). In order to explore the limit of $`\nu \sigma (r_M)`$ at $`z=0`$ (when $`\nu =1.68`$ in this model), we consider a halo of mass $`M=10^6M_{}`$ ($`r_M=18`$ kpc, $`\sigma (r_M)=8.5`$), and a surrounding region with $`r=2r_M`$ ($`M(r)=8\times 10^6M_{}`$, $`\sigma (r)=7.4`$). Consider randomly chosen regions of size $`r`$ in the universe. Half of them will have a negative overdensity and half will have a positive overdensity (where, as elsewhere in this paper, we consider the linearly-extrapolated overdensity). Those with a positive perturbation will almost certainly be contained within a halo of mass at least $`M(r)`$, since the typical $`\delta `$ on this scale is of order $`\sigma (r)`$, which in turn is much greater than the threshold $`\nu `$ for forming a halo. Even for those rare perturbations that are smaller than $`\nu `$, the fact that the perturbation is positive increases the (already very high) chance for the formation of some halo with radius between $`r_M`$ and $`r`$. Clearly, among regions with positive perturbations on the scale $`r`$, very few will contain an isolated halo of size $`r_M`$ that is not contained within a larger halo.

On the other hand, fully half of such large regions have a negative perturbation, again with magnitude of order $`\sigma (r)`$. Such a negative perturbation implies that certainly no halo has formed on the scale $`r`$; together with the correlation among spheres of different radii, it also implies a very strong reduction in the chance that any halo has formed on a scale between $`r_M`$ and $`r`$. This double effect implies that isolated halos of size $`r_M`$ will be found almost exclusively within voids on the scale $`r`$. Indeed, in the limit of $`\nu /\sigma (r_M)0`$, $`\delta `$ on the scale $`r`$ should be negative with a magnitude at least comparable to $`\sigma (r_M)`$, in order for $`\delta `$ to remain negative at all scales from $`r`$ down to $`r_M`$. Thus, in this limit we should find that the ratio $`\delta (r)/\nu \mathrm{}`$. This contrasts with the Press-Schechter models from the previous section, which predict that the typical such case should have a positive $`\delta `$ on the scale $`r`$ (in fact, a $`\delta `$ only slightly below $`\nu `$); this behavior arises from the general failure of these models to properly incorporate negative density fluctuations.

To derive our new, self-consistent model, we apply the extended Press-Schechter model in Fourier space and then use an ansatz to convert it to real space. We begin by following the derivation of $`\delta (r)_{\mathrm{PS}\mathrm{k}}`$ in the previous subsection but now including the absorption barrier at $`\delta =\nu `$. Thus, we consider once again the two scales with variances $`S_{k,M}`$ and $`S_k`$, but we now require $`S_{k,M}`$ to form a halo that is not contained within any larger virialized halo; in particular this implies that there is no halo on the scale corresponding to $`S_k`$ or on any larger scale. Then the probability distribution of the overdensity $`\delta `$ at $`S_k`$ is given by $`Q(\nu ,\delta ,S_k)`$ as in eq. (5). As before, the random walk from $`S_k`$ onwards is independent of the random walk leading up to $`S_k`$. In order to account for the extended Press-Schechter definition of halo formation, we use arguments similar to those leading to eq. (7) \[see Bond et al. (1991) and Lacey & Cole (1993) for other applications of similar arguments\]; the probability distribution of the initial density at $`S_k`$, given the formation of a halo at $`S_{k,M}`$, is found to be

$$P(\delta |\delta _M=\nu )_{\mathrm{ePS}\mathrm{k}}=Q(\nu ,\delta ,S_k)\frac{f(\nu \delta ,S_{k,M}S_k)}{f(\nu ,S_{k,M})}.$$
(17)
The mean expected $`\delta `$ is given by integrating $`\delta `$ over this distribution, with limits of $`\mathrm{}`$ to $`+\nu `$; we thus obtain the result of using the extended Press-Schechter model in $`k`$-space:

$`{\displaystyle \frac{\delta (r)_{\mathrm{ePS}\mathrm{k}}}{\nu }}`$ $`=`$ $`1\left(1\alpha +{\displaystyle \frac{\alpha }{\beta }}\right)\mathrm{erf}\left[\sqrt{{\displaystyle \frac{\beta (1\alpha )}{2\alpha }}}\right]`$ (18)
$`\sqrt{{\displaystyle \frac{2\alpha (1\alpha )}{\pi \beta }}}\mathrm{exp}\left[{\displaystyle \frac{\beta (1\alpha )}{2\alpha }}\right],`$
where we have defined

$$\alpha =\frac{S_k}{S_{k,M}};\beta =\frac{\nu ^2}{S_{k,M}}.$$
(19)
Unlike the models in the previous subsection, the present model does not yield a Gaussian probability distribution for the density. To find the scatter, we first integrate the probability distribution given by eq. (17) and find the cumulative probability distribution $`P_C`$. Using the variable $`\gamma =\delta /\nu `$ instead of $`\delta `$, the probability of obtaining a value of $`\delta (r)`$ that is between $`\mathrm{}`$ and $`\nu \gamma `$ equals

$`P_C(\gamma )_{\mathrm{ePS}\mathrm{k}}`$ $`=`$ $`1+{\displaystyle \frac{1}{2}}\left[\mathrm{erf}(A)+\mathrm{erf}(B)\right]`$ (20)
$`+\sqrt{{\displaystyle \frac{\alpha }{2\pi \beta (1\alpha )}}}\left[e^{A^2}e^{B^2}\right],`$
where we use the shorthand

$$A=\frac{(\gamma \alpha )\sqrt{\beta }}{\sqrt{2\alpha (1\alpha )}};B=\frac{(\gamma +\alpha 2)\sqrt{\beta }}{\sqrt{2\alpha (1\alpha )}}.$$
(21)
We can now numerically calculate the $`1\sigma `$ interval of $`\delta (r)`$ by bracketing the central $`68.3\%`$ of the probability. For a given $`M`$ and $`\nu `$, we thus define the minus $`1\sigma `$ density profile as given by values of $`\gamma `$ determined at each $`r`$ by $`P_C(\gamma )=15.9\%`$, while the plus $`1\sigma `$ profile is determined by $`P_C(\gamma )=84.1\%`$.

Having developed several possible infall models, we compare their relative merits by considering several limiting cases. First, when $`rr_M`$ (or equivalently $`S_kS_{k,M}`$), all the models satisfy the continuity condition $`\delta (r)/\nu 1`$ with zero scatter. In the limit of $`\beta 0`$, we require $`\delta (r)/\nu \mathrm{}`$ as explained earlier in this section. While $`\delta (r)_{\mathrm{ePS}\mathrm{k}}`$ behaves correctly in this limit (and approaches a constant that does not depend on $`\nu `$), $`\delta (r)_{\mathrm{PS}\mathrm{r}}`$ and $`\delta (r)_{\mathrm{PS}\mathrm{k}}`$ are proportional to $`\nu `$ and thus fail to correctly describe the $`\beta 0`$ limit. Thus, only $`\delta (r)_{\mathrm{ePS}\mathrm{k}}`$ behaves physically in both of the limits $`rr_M`$ and $`\beta 0`$.

In the limit of a halo that represents an extremely rare peak, i.e., $`\beta \mathrm{}`$, the Press-Schechter model becomes an accurate description that is equivalent to the extended Press-Schechter model; the reason is that if the central halo is extremely rare then it is exceedingly improbable to find a virialized halo on any scale larger than $`r_M`$, and thus it is acceptable to neglect the barrier at $`\delta =\nu `$ except right near the scale $`r_M`$. In this limit, $`\delta (r)_{\mathrm{ePS}\mathrm{k}}`$ and $`\delta (r)_{\mathrm{PS}\mathrm{k}}`$ both approach the value $`\nu \alpha `$. However, in this limit only, the PS-r model is the best available model since it directly considers spheres in real space, and these correspond most directly to the regions that actually collapse to form halos; this model also applies most directly to the threshold $`\nu `$ that is derived based on spherical collapse in real space. Indeed, the best possible model would be $`\delta (r)_{\mathrm{ePS}\mathrm{r}}`$, a quantity based on the extended Press-Schechter model applied directly in real space. However, such a model would be extremely difficult to develop, since excursion sets of $`\delta `$ as a function of $`\sigma ^2(r)`$ would no longer correspond to pure random walks \[see Bond et al. (1991) for related discussion\]. Thus, following the tradition of the extended Press-Schechter approach (see §1 and 2.1), we use $`\delta (r)_{\mathrm{ePS}\mathrm{k}}`$ along with an ansatz.

To arrive at the most logical ansatz, we simply replace $`\alpha `$ and $`\beta `$ in eq. (19) with their $`r`$-space equivalents. Based on the above arguments, we adopt $`\delta (r)_{\mathrm{PS}\mathrm{r}}`$ and $`\mathrm{\Sigma }(r)_{\mathrm{PS}\mathrm{r}}`$ as the correct $`\beta \mathrm{}`$ limits of the mean and the standard deviation, respectively. Indeed, $`\delta (r)`$ in the ePS-k model approaches a Gaussian distribution in the limit of $`\beta \mathrm{}`$, and therefore we can match this limit identically to the PS-r model, and this procedure fixes unambiguously the required $`r`$-space versions of both $`\alpha `$ and $`\beta `$. Thus, our final ansats for the mean initial profile of the density averaged in spheres of radius $`r`$, linearly extrapolated to $`z=0`$, around a virialized halo is

$$\frac{\delta (r)}{\nu }\frac{\delta (r)_{\mathrm{ePS}\mathrm{k}}}{\nu },$$
(22)
where we use eq. (18) except that we define

$$\alpha \frac{\xi _r(r_M,r)}{\sigma ^2(r_M)};\beta \frac{\nu ^2\alpha (1\alpha )}{\sigma ^2(r)\alpha \xi _r(r_M,r)}.$$
(23)
We apply these same definitions to eqs. (20) and (21) in order to estimate the scatter in the infall profiles as given by this final model.

### 2.4 Halo Virialization in the Presence of Infall

As noted in §2.1, the linearly extrapolated initial overdensity of a spherical tophat corresponding to the moment of collapse is (Peebles, 1980) $`\delta _c=3(12\pi )^{2/3}/20`$ if $`\mathrm{\Omega }_m=1`$. In reality, even a slight violation of the exact symmetry of the initial perturbation can prevent the tophat from collapsing to a point. Instead, the halo reaches a state of virial equilibrium by violent relaxation (phase mixing). Using the virial theorem $`U=2K`$ to relate the potential energy $`U`$ to the kinetic energy $`K`$ in the final state, the final overdensity relative to the cosmic mean density at the collapse redshift is (Peebles, 1980) $`\mathrm{\Delta }_c=18\pi ^2`$ in the Einstein-de Sitter model. Thus, if a tophat containing mass $`M`$ is predicted by spherical collapse to reach a physical radius $`R=0`$ at the collapse redshift $`z`$, this amount of mass is instead assumed to be distributed within a sphere of size equal to the (physical) virial radius $`R_{\mathrm{vir}}`$ determined by $`M=\frac{4}{3}\pi \overline{\rho }(z)\mathrm{\Delta }_cR_{\mathrm{vir}}^3`$.

In the presence of infall, however, the enclosed mass is modified. Indeed, at the moment when the tophat mass $`M`$ is predicted by spherical collapse to reach a radius $`R=0`$, the surrounding infall implies that additional mass has already accumulated within the sphere of radius $`R_{\mathrm{vir}}`$. Thus, the overdensity within the standard virial radius $`R_{\mathrm{vir}}`$ is predicted to be higher than the standard value of $`18\pi ^2`$. Therefore, once infall is included, even in the Einstein-de Sitter model the precise overdensity to be used to define a “virialized” halo becomes somewhat uncertain. In this paper we choose the overdensity of $`18\pi ^2`$ to define the virial radius in all cosmological models, since this overdensity has been frequently used to compare the analytical model predictions (e.g., for the halo mass function) with numerical simulations (e.g., Jenkins et al., 2001). Thus, we carry out spherical collapse calculations with the average infall profile derived in the previous section. Having fixed the virial overdensity we then find the enclosed virial mass $`M_{\mathrm{vir}}`$, which is greater than the original tophat mass $`M`$ by up to a factor of two (depending on $`M`$, $`z`$, and the cosmological parameters). We also find the linearly extrapolated initial overdensity of the virialized mass shell, which we denote $`\delta _v`$.

In the Einstein-de Sitter case we find $`\delta _v=1.5927`$ for the infall-based definition of virialization just noted, compared to the standard value of $`\delta _c=1.6865`$ for an isolated initial tophat perturbation. More generally, $`\delta _c`$ depends on the cosmological parameters. In a flat cosmology with $`\mathrm{\Omega }_m(z)+\mathrm{\Omega }_\mathrm{\Lambda }(z)=1`$ at a redshift $`z`$, we find the following fitting function accurate to $`5\times 10^5`$ in the range $`0.05<\mathrm{\Omega }_m(z)<1`$ (see also, e.g., Eke, Cole & Frenk, 1996):

$$\delta _c=1.6865\left[1+\frac{x}{1000}\left(4.86.51x.113x^2\right)\right],$$
(24)
where $`x=\mathrm{log}\left[\mathrm{\Omega }_m(z)\right]`$. Although $`\delta _v`$ does not strictly depend only on the cosmological parameters, we find numerically that it is essentially independent of halo mass. Indeed, over a wide parameter range of $`10^4M_{}<M<10^{16}M_{}`$ and $`0.05<\mathrm{\Omega }_m(z)<1`$, $`\delta _v`$ is given to within a relative accuracy of $`3\times 10^4`$ by

$$\delta _v=1.5927\left[1+\frac{x}{1000}\left(3.9.36x.09x^2\right)\right].$$
(25)
Thus, for example, $`\delta _v=1.573`$ for $`\mathrm{\Omega }_m(z)=0.05`$, compared to $`\delta _c=1.659`$ in this case.

We note that the Press-Schechter halo mass function is known to underestimate the abundance of rare halos relative to the results of numerical simulations (e.g., Jenkins et al., 2001). This can be fixed in the Einstein-de Sitter model if the collapse threshold is multiplied by a factor of 0.87, while we find that infall lowers the threshold, from $`\delta _c`$ to $`\delta _v`$, by a factor of only 0.94. Thus, infall by itself cannot fully explain the discrepancy in the Press-Schechter prediction for the number of high-mass halos.

## 3 Results

We illustrate results for the standard cosmological parameters $`\mathrm{\Omega }_m=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, $`h=0.7`$ and a CDM power spectrum (with initial slope $`n=1`$) normalized to a variance $`\sigma _8=0.8`$ in spheres of radius 8 $`h^1`$ Mpc. We consider a wide range of redshifts but focus on the higher values since at $`z<2`$ the infalling intergalactic gas is very tenuous and may thus be undetectable; however, note that at low redshift a significant fraction of the infalling baryons may already be in the form of virialized dwarf galaxies that are themselves detectable.

Figure 1 shows the predictions of our model \[eqs. (22) and (23)\] for the mean initial profiles of enclosed density around virialized halos. In most cases, the profiles lie significantly below the asymptotic case of a rare halo ($`\beta \mathrm{}`$) for which the profile is given by $`\delta (r)_{\mathrm{PS}\mathrm{r}}`$ as in eq. (13). In particular, when $`M_{\mathrm{vir}}=10^8M_{}`$, a halo forming at $`z2`$ is most likely to be found in a void on all scales larger than a few times $`r_M`$. Note that we restrict all calculations to radii where the initial profile of enclosed density is monotonically decreasing, since otherwise shells would cross during the subsequent collapse; this cutoff affects only the $`z=0`$ cases as well as the $`M_{\mathrm{vir}}=10^8M_{}`$ halo forming at $`z=2`$ or 3.

We calculate spherical collapse with the initial profiles shown in Figure 1. The resulting final density profiles are shown in Figure 2. In the limiting case of an extremely rare halo (i.e., very high formation redshift at a fixed halo mass), low mass halos have higher densities of dark matter surrounding their virial radius. This can be understood roughly as follows. For a power spectrum of density fluctuations approximated over some range of scales as a power law $`P(k)k^n`$, the correlation function over the same range of scales varies as $`\xi (r)r^{(n+3)}`$. The initial density profile around a rare halo follows, roughly, the dark matter correlation function. On small scales, the CDM power spectrum approaches $`n=3`$ and the initial density profile is rather flat, while on larger scales the power index increases toward $`n=2`$ for large galaxies, which are therefore surrounded by an initial profile that falls off more steeply (see Figure 1). These differences are magnified by the process of gravitational collapse, resulting in final values for $`\rho /\overline{\rho }`$ at the virial radius of 41, 34, and 27 for $`M_{\mathrm{vir}}=10^8`$, $`10^{10}`$, and $`10^{12}M_{}`$, respectively (see Figure 2). However, at moderate redshifts high mass halos represent rarer, more extreme peaks than do low mass halos. Thus, e.g., even at redshift 5 the profile surrounding $`M_{\mathrm{vir}}10^{10}M_{}`$ halos falls significantly below the profile of the extreme limiting case.

Next, we use eq. (20) to predict the expected scatter in the density profiles around virialized halos. Figure 3 shows the initial and final density profiles for $`M_{\mathrm{vir}}=10^{10}M_{}`$ at two different redshifts, for the mean expected profile that we have thus far considered as well as for the plus and for the minus $`1\sigma `$ profiles. The final density at $`R_{\mathrm{vir}}`$, for example, is higher than the mean density by a factor of $`12_4^{+8}`$ at $`z=2`$ and $`28_7^{+12}`$ at $`z=10`$. We thus find a significant scatter, although it is smaller at the higher redshift. If the scatter as calculated within our spherical model is representative of the true scatter among halos, its magnitude suggests that our predicted mean profile can be compared accurately to an observed sample if the sample contains at least a few dozen halos of comparable masses and redshifts.

Figure 4 shows the final profiles of peculiar velocity $`V_{\mathrm{pec}}`$ of the infalling dark matter corresponding to the halos considered in Figures 1 and 2. The same trends are seen in this figure as in the previous ones; in the extreme case of a rare halo, infall is stronger for low halo masses, but at moderate redshifts, the higher-mass halos actually produce the stronger infall. Figure 4 also shows the peculiar velocity that would correspond to stationary matter, i.e., with zero total velocity (long-dashed curve). This curve is given by minus the Hubble expansion velocity $`V_{\mathrm{Hub}}`$, which in these units is

$$\frac{V_{\mathrm{Hub}}(R)}{V_c}=\frac{1}{3\pi }\frac{R}{R_{\mathrm{vir}}}.$$
(26)
Values of $`V_{\mathrm{pec}}`$ that lie below $`V_{\mathrm{Hub}}`$ give a total velocity corresponding to infall.

When gas falls into a dark matter halo, incoming streams from all directions strike each other at supersonic speeds, creating a strong shock wave. Three-dimensional hydrodynamic simulations show that the most massive halos at any time in the universe are indeed surrounded by strong, quasi-spherical accretion shocks (Keshet et al., 2002; Abel, Bryan, & Norman, 2002). The radius of this shock depends in general on the history of gas infall into the halo, and will thus vary somewhat among halos depending on the history of their assembly through mergers. In particular, in the spherically-symmetric secondary infall solution, when a trace amount of adiabatic gas is allowed to flow in a gravitational potential determined by collisionless dark matter with density $`\mathrm{\Omega }_m=1`$, the accretion shock forms at the radius $`R_{90.8}`$ which encloses a mean dark matter density of 90.8 times the cosmic mean (Bertschinger, 1985). However, the simulation of Abel, Bryan, & Norman (2002) shows that along most lines of sight the shock occurs at a smaller distance, close to the virial radius.

The position of the accretion shock is an important theoretical prediction that can be probed by observations, since the infalling gas that is about to cross the shock front produces a sharp absorption feature in high-redshift quasar spectra (Barkana & Loeb, 2002). In order to illustrate the range of possible properties of this gas element that lies just beyond the shock front, we use the results mentioned above as a guide and consider two cases for the shock radius $`R_{\mathrm{sh}}`$: the smaller radius $`R_{\mathrm{sh}}=R_{\mathrm{vir}}`$ or the larger radius $`R_{\mathrm{sh}}=R_{90.8}`$. We assume that the pressure gradient is negligible in the pre-shock gas compared to gravity, and therefore the gas density and velocity are equal to those of the dark matter at distances beyond the shock front.

Figure 5 shows the density and total infall velocity of the gas that lies just beyond the shock front. The density and the magnitude of the infall velocity both tend to increase with redshift, and follow only weak trends with halo mass at a given redshift. The differences between the two redshifts shown are fairly significant for low mass halos but rather small for high mass halos. In general, the infalling gas is significantly overdense even for halo masses corresponding to $`z=0`$ dwarf galaxies, especially at the high redshifts which are increasingly becoming the focus of observations. The figure allows us to compare two sources of scatter, namely the varying position of $`R_{\mathrm{sh}}`$ and the variability in the initial density profile. The two effects produce a comparable scatter in the pre-shock density, while the infall velocity, which depends on the enclosed mean density rather than on the density itself, is primarily sensitive to the position of the shock.

## 4 Conclusions

We have developed a model for calculating the initial density profile around overdensities that later collapse to form virialized halos. This is the first such model to account for the possibility that an overdensity on a given scale may be contained inside a large overdensity on an even bigger scale. As a result, we have found that the mean expected profile \[eqs. (18) and (23)\] depends on both the mass of the halo and its formation redshift. Starting from this mean initial profile, we have used spherical collapse to obtain the final density profile surrounding the virialized halo.

In reality, there will be some variation in the initial profiles, with each leading to a different final profile. We have estimated the scatter within our spherical model, and found that a sample of around a few dozen halos is required in order to obtain an accurate mean profile by averaging. Halo samples derived from numerical simulations or from observations can be used to test our predictions for the mean and for the scatter of the density and velocity profiles of infalling matter around virialized halos.

Once infall is included, the standard derivation of the mean density enclosed within the virial radius fails. Redefining the virial radius as the radius enclosing a density higher than the cosmic mean by the standard value of $`18\pi ^2`$, we have found that the initial overdensity needed for a halo to form at some final redshift is lower than it would have to be for an isolated tophat perturbation \[compare eqs. (24) and (25)\]. Therefore, accounting for infall increases the predicted abundance of rare halos for a given initial power spectrum, but the increase is too small to fully explain the discrepancy between the Press-Schechter prediction for the number of high-mass halos and the number seen in numerical simulations.

We have confirmed that rare halos at a given redshift are surrounded by large, overdense regions of infall, but we have found that the more numerous halos that correspond to low-sigma peaks are surrounded by small infall regions that lie within large voids. However, the infalling gas just beyond the cosmological accretion shock is in general much denser than the cosmic mean; a density larger by a factor of at least 10 occurs around the most massive halos at $`z2`$. Possible applications of our results include the Lyman-$`\alpha `$ absorption profiles of high-redshift objects and the proximity effect of quasars at all redshifts.

## Acknowledgments

I am grateful for the hospitality of the Institute for Advanced Study, where this work was begun. I thank Avi Loeb and Evan Scannapieco for related discussions. I acknowledge the support of an Alon Fellowship at Tel Aviv University, Israel Science Foundation grant 28/02, and NSF grant AST-0204514.


format: text

text:
# Viscous cosmological models and accelerated Universes

## Abstract

It is shown that a present acceleration with a past deceleration is a possible solution of the Friedmann equation by considering the Universe as a mixture of a scalar with a matter field and by including a non-equilibrium pressure term in the energy-momentum tensor. The dark energy density decays more slowly with respect to the time than the matter energy density does. The inclusion of the non-equilibrium pressure leads to a less pronounced decay of the matter field with a shorter period of past deceleration.

PACS: 98.80.Cq

According to recent cosmological observations there exists an evidence that the Universe is flat (see Sievers et al ) and expanding with a positive acceleration (see Perlmutter et al  and Riess et al ). The flatness of the Universe is connected with the total density parameter $`\mathrm{\Omega }_{\mathrm{tot}}`$ which is the sum of the density parameters related to vacuum energy $`\mathrm{\Omega }_\mathrm{\Lambda }`$, cold dark matter $`\mathrm{\Omega }_{\mathrm{cdm}}`$ and baryons $`\mathrm{\Omega }_b`$, i.e., $`\mathrm{\Omega }_{\mathrm{tot}}=\mathrm{\Omega }_\mathrm{\Lambda }+\mathrm{\Omega }_{\mathrm{cdm}}+\mathrm{\Omega }_b`$. The sum of the density parameters of cold dark matter and baryons refers to the matter density parameter $`\mathrm{\Omega }_m=\mathrm{\Omega }_{\mathrm{cdm}}+\mathrm{\Omega }_b`$. Most recent values for the density parameters from measurements of the anisotropy of the cosmic microwave background (CMB) are :

$$\mathrm{\Omega }_{\mathrm{tot}}=1.01_{0.06}^{+0.09};\mathrm{\Omega }_m=0.37\pm 0.11\text{and}\mathrm{\Omega }_b=0.060\pm 0.020.$$

Measurements of the redshift of the type Ia supernova SN 1997ff indicate a present acceleration of the Universe and - according to Turner and Riess  \- provide a past decelerating phase of the Universe.

Phenomenological cosmologies play an important role in the understanding of the evolution of the Universe. A remarkable combination of general relativity and thermodynamics allow the description of different regimes in cosmological theories. Models that include the recent experimental evidence of a present accelerated Universe therefore deserve an increasing attention.

On the other hand, there exist other regimes where a positive acceleration is present. Indeed, in inflationary models an exponential expansion with positive acceleration for the early Universe is proposed where a hypothetical particle, the inflaton, is at the core of that mechanism. These ideas play the role of a precursor for the relation between vacuum (or dark) energy and a cosmological constant that depends on time. For recent reviews on dark energy one is referred to the works by Peebles and Ratra  and Turner .

What we know today is that dark energy dominates the composition of our Universe and is the responsible for the positive accelerated expansion. It interacts weakly with ordinary matter (baryonic matter represents at maximum 5% of the whole mass/energy composition of our Universe).

The point we would like to discuss in this work is that a past decelerated with a present accelerated expansion of the Universe can be found as a solution of the Friedmann equation. The Universe is modeled as a mixture of two constituents namely, a scalar field that represents the dark energy and a matter field that describes the baryonic matter and the cold dark matter; in our approach the interaction between the different constituents refers only to the one between the dark energy and the matter via the gravitational field which is represented by the cosmic scale factor $`a(t)`$. We model this interaction by using the ideas of a thermodynamic theory, by means of the inclusion of a non-equilibrium pressure term in the energy-momentum tensor which represents an irreversible process of energy transfer between the matter and the gravitational fields. Within the framework of first order Eckart thermodynamic theory the non-equilibrium pressure is proportional to the Hubble parameter and its proportionality factor is identified with the bulk viscosity, so that this model is referred in the literature as a viscous cosmological model (see, for example, the works ).

Although the transition from a past deceleration to a present acceleration of the Universe is found as a possible solution of the Friedmann equation independently of the presence of the non-equilibrium pressure term, it is very questionable to model the Universe as a perfect gas mixture of scalar and matter fields evolving without dissipative effects (for further discussions on this subject one is referred to the works ). In this work, it is shown that the dark energy density decays more slowly with respect to the time than the matter energy density does. These conclusions are valid for both cases where the non-equilibrium pressure is present or absent. However, the inclusion of a non-equilibrium pressure leads to a less pronounced decay of the matter field with a shorter period of past deceleration, since the non-equilibrium pressure is the responsible for the energy transfer between the matter and the gravitational fields.

We shall consider a spatially flat, homogeneous and isotropic Universe described by the Robertson-Walker metric. With respect to the four-velocity $`U^\mu `$ the energy-momentum tensor $`T^{\mu \nu }`$ is decomposed as111Units have chosen so that $`c=1`$.

$$T^{\mu \nu }=(\rho _X+\rho _m+p_X+p_m+\varpi )U^\mu U^\nu (p_X+p_m+\varpi )g^{\mu \nu }.$$
(1)
where $`p_X`$ is the pressure of the dark energy, $`p_m`$ is the pressure of the matter field and $`\varpi `$ denotes the non-equilibrium pressure which is coupled to the irreversible process of energy transfer between matter and gravitational fields.

In a comoving frame the conservation law of the energy-momentum tensor $`T_{}^{\mu \nu }{}_{;\nu }{}^{}=0`$ leads to the balance equation for the energy density

$$\dot{\rho }_X+\dot{\rho }_m+3H(\rho _X+\rho _m+p_X+p_m+\varpi )=0,$$
(2)
where $`H=\dot{a}/a`$ is the Hubble parameter and the dot refers to a differentiation with respect to the time.

The identification of the dark energy with a scalar field allow us to assume that the balance equation for the dark energy density is given by

$$\dot{\rho }_X+3H(\rho _X+p_X)=0.$$
(3)
Hence, the balance equation for the dark energy density decouples from that for the energy density of the matter and we have from (2) and (3)

$$\dot{\rho }_m+3H(\rho _m+p_m)=3H\varpi .$$
(4)
The equation that connects the evolution of the cosmic scale factor $`a(t)`$ with the energy densities of the scalar and matter fields is given by the Friedmann equation

$$H^2=\frac{8\pi G}{3}(\rho _X+\rho _m),$$
(5)
where $`G`$ is the gravitational constant.

In a previous work  we have calculated the energy-momentum pseudotensor of the gravitational field in a flat Robertson-Walker metric and found

$$T_G^{00}=\frac{3}{8\pi G}\left(\frac{\dot{a}}{a}\right)^2.$$
(6)
If we identify the $`T_G^{00}`$ with the energy density $`\rho _G`$ of the gravitational field we can get thanks to the Friedmann equation (5) the relationship $`\rho _G=(\rho _X+\rho _m)`$.

Now we differentiate $`\rho _G=(\rho _X+\rho _m)`$ with respect to the time and by using the balance equation for the energy density (4) we get

$$\dot{\rho }_G+3H(\rho _Gp_mp_X)=3H\varpi ,$$
(7)
which can be interpreted as a balance equation for the energy density of the gravitational field. Moreover, if we compare equations (4) and (7) we infer that the right-hand side of both equations has a term (with opposite sign) proportional to the non-equilibrium pressure which is responsible to the transfer of energy between the gravitational and matter fields.

It is usual to assume for the pressure of the dark energy the equation of state $`p_X=w_X\rho _X`$, with $`w_X<1/3`$ (see Peebles and Ratra ). Hence, one can obtain by integration of (3)

$$\frac{\rho _X}{\rho _X^0}=\left(\frac{a_0}{a}\right)^{3(w_X+1)},$$
(8)
where $`\rho _X^0`$ and $`a_0`$ are the values of the dark energy density and of the cosmic scale factor at $`t=0`$ (by adjusting clocks), respectively.

In order to determine the time evolution of the dark energy density one has to know the evolution of the cosmic scale factor which is determined from the system of equations (4), (5) and (8). This system of equations is closed by assuming a relationship between the pressure and the energy density of the matter and a constitutive equation for the non-equilibrium pressure $`\varpi `$. We assume a barotropic equation of state for the pressure $`p_m=w_m\rho _m`$, with $`0w_m1`$ and a linear relationship between the non-equilibrium pressure $`\varpi `$ and the Hubble parameter $`H`$ within the framework of Eckart first order thermodynamic theory

$$\varpi =3\eta H,\eta =\alpha (\rho _X+\rho _m).$$
(9)
In the above equation $`\eta `$ is the coefficient of bulk viscosity which is consider to be proportional to the energy density of the mixture and $`\alpha `$ is a constant.

We differentiate the Friedmann equation (5) with respect to the time and get by the use of (8) and (9)

$$\dot{H}=\frac{3}{2}\left[\frac{(w_mw_X)}{1+\rho _m^0/\rho _X^0}\left(\frac{1}{a}\right)^{3(w_X+1)}+(3\alpha Hw_m1)H^2\right].$$
(10)
Above $`\rho _m^0`$ is the energy density of the matter field at $`t=0`$ (by adjusting clocks). Moreover, all terms in (10) are dimensionless quantities defined by

$$H\frac{H}{H_0},ttH_0,a\frac{a}{a_0},\alpha \alpha H_0,\text{with}H_0=\sqrt{\frac{8\pi G}{3}(\rho _X^0+\rho _m^0)}.$$
(11)
Equation (10) is a second-order differential equation for the cosmic scale factor $`a(t)`$ which is a function of four parameters, namely $`w_m`$, $`w_X`$, $`\alpha `$ and $`\rho _m^0/\rho _X^0`$. The solution of the differential equation (10) is found by specifying values for two initial conditions and for the four parameters.

Instead of using the constitutive equation (9) for the non-equilibrium pressure one may consider it as a variable within the framework of extended (causal or second-order) thermodynamic theory. In this case the evolution equation for the non-equilibrium pressure - in a linearized theory - reads (see, for example )

$$\varpi +\tau \dot{\varpi }=3\eta H,$$
(12)
where $`\tau `$ is a characteristic time. Here we follow the works  and assume that the characteristic time is given by $`\tau =\eta /\rho `$.

Hence, we obtain from the Friedmann equation (5) and from the evolution equation for the non-equilibrium pressure (12) the following system of differential equations

$$\dot{H}=\frac{3}{2}\left[\frac{(w_mw_X)}{1+\rho _m^0/\rho _X^0}\left(\frac{1}{a}\right)^{3(w_X+1)}(w_m+1)H^2+\varpi \right],$$
(13)
$$\varpi +\alpha \dot{\varpi }=3\alpha H^3,$$
(14)
thanks to (9)2. Apart from the dimensionless quantities introduced above, $`\varpi 8\pi G\varpi /(3H_0^2)`$ is a dimensionless non-equilibrium pressure.

If we specify values for three initial conditions and for the four parameters $`w_m`$, $`w_X`$, $`\alpha `$ and $`\rho _m^0/\rho _X^0`$ we can obtain from the system of equations (13) and (14) the time evolution of the cosmic scale factor $`a(t)`$ and of the non-equilibrium pressure $`\varpi (t)`$.

Once $`a(t)`$ is determined from (10) or from the system (13) and (14) one can find the energy densities by using the expressions

$$\rho _X=\left(\frac{1}{a}\right)^{3(w_X+1)},\rho _m=\left(1+\frac{\rho _X^0}{\rho _m^0}\right)H^2\frac{\rho _X^0}{\rho _m^0}\left(\frac{1}{a}\right)^{3(w_X+1)},$$
(15)
where $`\rho _X\rho _X/\rho _X^0`$ and $`\rho _m\rho _m/\rho _m^0`$ are dimensionless quantities.

The initial conditions we choose at the instant of time $`t=0`$ (by adjusting clocks) are: $`a(0)=1`$ for the cosmic scale factor and $`H(0)=1`$ for the Hubble parameter and $`\dot{\varpi }(0)=0`$ for the non-equilibrium pressure. Since the differential equations (10), (13) and (14) do depend on the four parameters $`w_m`$, $`w_X`$, $`\alpha `$ and $`\rho _m^0/\rho _X^0`$ there exists much freedom to find their solutions. In order to obtain the graphics in the Fig. 1 we have chosen:

* $`w_X=0.7`$ so that the condition $`w_X<1/3`$ holds;
* $`w_m=0.2`$, so that $`0<w_m<1/3`$ where the value $`w_m=0`$ refers to dust and $`w_m=1/3`$ to radiation;
* $`\rho _m^0/\rho _X^0=2`$, i.e., the amount of the matter energy density is twice the one of the dark energy density;
* $`\alpha =0`$ when the non-equilibrium pressure is absent and $`\alpha =0.05`$ (say) when it is present.

We have plotted in Fig. 1 the dark energy density $`\rho _X`$, the matter energy density $`\rho _m`$ and the acceleration $`\ddot{a}`$ as functions of the time $`t`$ for the three cases: a) by using a second-order thermodynamic theory (straight lines); b) by using a first-order thermodynamic theory (dashed lines) and c) by neglecting the non-equilibrium pressure (perfect fluid solution). We infer from this figure that for the parameters chosen above the behavior of the solutions of the first-order and second-order thermodynamic theories is practically the same. In the three cases there exist periods where the Universe is decelerating and accelerating. Moreover, the dark energy density decays more slowly than the matter energy density does. Even when a small amount of dark energy density - with respect to the matter energy density - is taken into account, these fields evolve in such a manner that for large times the amount of dark energy density is very large with respect to the energy density of the matter field. The inclusion of the non-equilibrium pressure lead to: i) a slower decay of the matter field with respect to the time, which is a consequence of the irreversible process of energy transfer between the matter and gravitational fields; ii) a smaller interval of past deceleration and iii) a more rapid decay of the dark energy density.

Other conclusions can be obtained by changing the values of the parameters. In the case where there is no dark energy density $`\rho _X=0`$ there is no period of acceleration, i.e., only a period of deceleration is possible. By changing the interval $`w_m`$ to $`1/3<w_m<2/3`$ \- which corresponds to the interval between radiation $`w_m=1/3`$ and non-relativistic matter $`w_m=2/3`$ \- the period of past deceleration increases. This last result is also found when $`w_X`$ increases, i.e., for $`w_X>0.7`$. By increasing the value of the dimensionless constant $`\alpha `$ the effect of the non-equilibrium pressure is more pronounced and the period of past deceleration decreases.