A&A 438, 443-460 (2005)
DOI: 10.1051/0004-6361:20052885

Adhesive gravitational clustering

T. Buchert^1,2,3 - A. Domínguez^4,5,6

1 - Arnold Sommerfeld Center for Theoretical Physics, Ludwig-Maximilians-Universität, Theresienstr. 37, 80333 München, Germany
2 - Theory Division, CERN, 1211 Genève 23, Switzerland
3 - Observatoire de la Côte d'Azur, Lab. G.D. Cassini, BP 4229, 06304 Nice Cedex 4, France
4 - Max-Planck-Institut für Metallforschung, Heisenbergstr. 3, 70569 Stuttgart, Germany
5 - Institut für Theor. und Angew. Physik, Univ. Stuttgart, Pfaffenwaldring 57, 70569 Stuttgart, Germany
6 - Física Teórica, Univ. Sevilla, Apdo. 1065, 41080 Sevilla, Spain

Received 16 February 2005 / Accepted 7 March 2005

Abstract
The notion of adhesion has been advanced for the phenomenon of stabilization of large-scale structure emerging from gravitational instability of a cold medium. Recently, the physical origin of adhesion has been identified: a systematic derivation of the equations of motion for the density and the velocity fields leads naturally to the key equation of the "adhesion approximation'' - however, under a set of strongly simplifying assumptions. In this work, we provide an evaluation of the current status of adhesive gravitational clustering and a clear explanation of the assumptions involved. Furthermore, we propose systematic generalizations with the aim to relax some of the simplifying assumptions. We start from the general Newtonian evolution equations for self-gravitating particles on an expanding Friedmann background and recover the popular "dust model'' (pressureless fluid), which breaks down after the formation of density singularities; then we investigate, in a unified framework, two other models which, under the restrictions referred to above, lead to the "adhesion approximation''. We apply the Eulerian and Lagrangian perturbative expansions to these new models and, finally, we discuss some non-perturbative results that may serve as starting points for workable approximations of non-linear structure formation in the multi-stream regime. In particular, we propose a new approximation that includes, in limiting cases, the standard "adhesion model'' and the Eulerian as well as Lagrangian first-order approximations.

Key words: gravitation - methods: analytical - cosmology: theory - cosmology: large-scale structure of Universe

1 Gravitational clustering

The present work aims to push analytical modeling of structure formation into a regime that may be placed between the formation epoch of large-scale structure and the onset of virialization of gravitationally bound objects. Phenomenologically, this regime is characterized by a stabilization of structures that formed out of gravitational instability of a cold medium. The physical origin of this adhesive clustering effect is the balance of gravitational forces and dynamical stresses in collisionless matter. The latter arise from a subsequently establishing multi-stream hierarchy within collapsing high-density regions. Although multi-stream forces tend to disperse structures, the resulting effect together with gravity tends to stabilize them. This regime is summarized by the term non-dissipative gravitational turbulence advanced by Gurevich & Zybin (1995). The models we investigate are a focus of current research, since efforts to simulate Hubble volumes of the Universe and to understand galaxy halo formation are faced within a single approach. However, we take a more conservative point of view to understand the evolution of structure on galaxy cluster scales (for a recent theoretical attempt to address halo structure within kinetic theory see Ma & Bertschinger 2004).

The current status of analytical models concerning large-scale structure formation may be centred on Zel'dovich's approximation together with its foundations (the Lagrangian perturbation theory) and optimizations using filtering techniques for initial perturbation spectra (Zel'dovich 1970, 1973; Shandarin & Zel'dovich 1989; Buchert 1989, 1992; Melott 1994; Bouchet et al. 1992, 1995; Sahni & Coles 1995; Buchert 1996; Ehlers & Buchert 1997). These schemes are capable of modeling the evolution of generic spectra including Cold-Dark-Matter cosmogonies down to galaxy cluster scales (Bouchet et al. 1995; Buchert et al. 1994; Melott et al. 1995; Weiß et al. 1997; Hamana 1998). Above these scales the optimized Lagrangian schemes roughly reproduce the results of N-body simulations. Since the exact solution is not known, we seek agreement between the two modeling techniques used. Below these scales, both techniques tend to fail as a result of both poorly understood physics and poor resolution power, respectively. Agreement between N-body simulations and simple solutions of Lagrangian perturbation schemes may be considered as supporting analytical models, but they could as well be considered as a drawback of N-body simulations given the simplicity of the analytical approximations compared with the complexity of non-linear self-gravity. A deeper understanding of structure formation below galaxy cluster scales may require more than improving spatial resolution. The use of N-body computing in cosmology is possibly overstated as is indicated by the differing results obtained when using different N-body codes on the scales of interest (Melott et al. 1997; Splinter et al. 1998). Furthermore, particular features show that caution is in order: (i) the validity of the mean field approximation for the gravitational field strength that is commonly used is questionable, especially on the scales of galaxy halos; (ii) the possible emergence of soliton states (Götz 1988) that arise as a result of the non-linear interaction between gravity and pressure-like forces; solitons have special stability properties and may dominate large-scale structure at late times; (iii) the behavior of the gravitational field strength near high-density regions: from a generic integral of the field equations (Buchert 1993a, Sect. 6.1) one clearly infers proportionality of the field strength and the density, whereas Lagrangian perturbation schemes (in the regime where they match N-body runs) show smooth and moderately increased field strengths when crossing high-density regions. In the extreme case of infinite resolution, Zel'dovich's approximation produces caustics in the density field, whereas the field strength remains smooth, although it should blow up at caustics; whether N-body codes treat this correctly is questionable in view of the strongly varying fields on the spatial as well as temporal resolution scales.

In this situation, analytical models that are capable of accessing non-linear scales should be further explored. This task is not easy, since - as explained above - such models may not be fully guided by comparisons with N-body results. The pioneering model for adhesive clustering has been built in view of the shortcomings of Zel'dovich's approximation predicting structure decay after their formation: a Laplacian forcing has been added ad hoc to the evolution equation for the peculiar-velocity. Such a force has the right property of holding structures together after their formation and, moreover, the model can be solved exactly in terms of known solutions of the 3D Burgers' equation (Gurbatov et al. 1989; Weinberg & Gunn 1990; Kofman et al. 1990, 1992, 1994; see also: Gurbatov et al. 1983, 1985, 1991). The phenomenon of non-dissipative gravitational turbulence is mimicked by the so-called "burgulence'' (Frisch & Bec 2001) and roughhewned into an "effective viscosity''. However, the performance of the "adhesion approximation'' does not dramatically improve the Lagrangian perturbation schemes, if the latter are subjected to optimization techniques as mentioned above: even for a strongly hierarchical power spectrum (power law index of n = -1) the optimized Lagrangian schemes yield better cross-correlation statistics for the density fields, compared with numerical simulations (Buchert 1999). However, comparisons have been conducted only on the level of large-scale structure (e.g., Sahni et al. 1994 investigate the evolution of voids); the role that the "adhesion approximation'' could play on subcluster scales has not been explored. Other proposals to analytically model the weakly non-linear regime tried to improve the performance for large-scale structure, e.g. "freezing'' the initial streamlines of the fluid keeping the velocity potential constant (so-called "Frozen Flow Approximation'', Matarrese et al. 1992), or treating the gravitational potential as constant ("Frozen Potential Approximation'', Brainerd et al. 1993; Bagla & Padmanabhan 1994). However, these approximations, by their nature, are unable to model a realistic evolution beyond the epoch when the time-dependence of the velocity or the potential is important, cf. the systematic comparison of different models carried out by Sathyaprakash et al. (1995).

Lying at the physical core of the problem of covering multi-streaming, the "adhesion approximation'' is a promising model due to a strong support given to it by Buchert & Domínguez (1998): a Laplacian forcing appears naturally in a physical description of adhesive gravitational clustering and the "adhesion approximation'' can be derived on the basis of a set of - however, strongly simplifying - assumptions. The Laplacian forcing does not represent a true viscosity, it arises by combining the gravitational field equations with a dynamical pressure gradient due to multi-stream stresses. Consequently, momentum is conserved and energy is not lost: the model is time-reversible. Multi-streaming implies a reshuffling of energy components: bulk kinetic energy is transformed into internal kinetic energy, and the gravitational potential energy of a high-density region gradually dominates over the tidal interaction with its environment, thus effectively isolating the system and - in idealized cases - establishing a dynamical equation of state as a relation between gravitational potential energy and internal kinetic energy.

This reasoning can be made mathematically precise by deriving the equations of motion for the coarse-grained (or filtered) density and velocity fields. Given a smoothing length L, the system is divided into the degrees of freedom of the scales above and below L, respectively. The evolution of the smoothed density and velocity fields is in general dynamically coupled to the degrees of freedom below L. The standard "dust model'' is recovered by assuming that the dynamical effect of this coupling is negligible. When the "dust model'' develops singularities, this assumption breaks down, and an improved model has to be used that accounts for multi-streaming and the coupling to the small-scale degrees of freedom. The model introduced by Buchert & Domínguez (1998), which we call the Euler-Jeans-Newton ("EJN'') model, identifies the main source of the corrections to "dust'' as the gravitational multi-stream (GM) effect. This gives rise to a stress tensor of purely kinetic origin - as in an ideal gas - which is added to the large-scale gravity already considered by the "dust model''. This stress may be idealized in the "EJN model'' phenomenologically as isotropic, and modelled with a pressure  $p(\varrho)$ depending only on the density $\varrho$, i.e. a dynamical equation of state.

Numerical tests (Domínguez 2003; Domínguez & Melott 2004) show that a polytropic equation of state, $p \propto \varrho^{2-\eta}$ in a certain range of densities is characteristic for a variety of initial conditions. The "anomaly'' $\eta$ measures the deviations from the naive virial prediction $p \propto \varrho^2$, depends itself on initial conditions ($\eta$ ranges from 0 to 1 as the small-scale power is reduced), and can be assigned to the fact that the fluid elements are not isolated and in a stationary state, as required when deriving $p \propto \varrho^2$ from the virial theorem. An exact integral of the form $p \propto \varrho^{5/3}$ can be deduced theoretically under certain strong assumptions (Buchert & Domínguez 1998).

Recent developments of the "EJN model'' concern the Lagrangian linear regime (Adler & Buchert 1999), where solutions can be found by extrapolating known results of the Eulerian linear regime; the Lagrangian scheme has been developed to second order (Morita & Tatekawa 2001; Tatekawa et al. 2002) and recently to third order (Tatekawa 2005); numerical tests employing N-body and hydrodynamical simulations have been conducted by Tatekawa (2004a,b).

While the foundations of the "EJN model'' can be derived on the basis of the velocity moment hierarchy of the commonly employed Vlasov-Poisson system, its practical implementation needs (phenomenological) closure conditions to truncate the hierarchy. An attempt to systematically go beyond phenomenology is the Small-Size Expansion ("SSE'') introduced by Domínguez (2000). The mode-mode coupling in the equations for the fields smoothed over the scale L is estimated under the assumption that the most important contribution to this coupling comes from inter-mode coupling at scales $\geq$L. The systematic nature of this scheme is due to a formal expansion of this coupling in powers of L. The "dust model'' is recovered as the lowest order term (formally setting L=0). The corrections also yield a Laplacian forcing as well as the "adhesion model'' under certain simplifying assumptions. However, the corrections account not only for the GM effect like the "EJN model'', but also for the gravitational influence of the density inhomogeneities on scales below L: the stress tensor has an additional potential-energy contribution that can be viewed as a correction to the mean field approximation for the gravitational field-strength and also provides a Laplacian forcing, although subdominant relative to the one due to the GM effect. This effect may become dominant if one aims to describe structure evolution on galaxy halo scales. The stress tensor in the "SSE model'' is also a source of vorticity by tidal torques and shear stretching. This has been studied recently with the Eulerian perturbation expansion (Domínguez 2002).

In the present work the "dust'', the "EJN'' and the "SSE'' models will be presented in a unified framework, as different closure approximations to a hierarchy of equations, and a clearcut description of the weakly non-linear regime of adhesive gravitational clustering will be offered. It is expected that the gravitational multi-stream process (GM-effect for short) as well as stresses arising from corrections to the mean field approximation may help to establish the stationary configurations that are usually studied in stellar systems theory, where "virialized objects'' are, however, considered as isolated entities.

Section 2 applies the coarse-graining method to cosmological structure formation in the phase space of N particles, resulting in a continuum description which features the Vlasov dynamics as a subcase. Section 3 describes the resulting fully non-linear evolution equations for the density and peculiar-velocity fields in Eulerian space and summarizes the assumptions that reduce the general problem to the key equation of the "adhesion approximation''. Section 4 reviews the application of the Eulerian perturbative expansion and Sect. 5 that of the Lagrangian perturbative expansion. Section 6 presents new results beyond the perturbative regime. Section 7 summarizes the results and proposes generalizations of the "adhesion approximation'' as well as future prospects.

Notation: Eulerian coordinates are denoted by ${\vec{x}}$, while Eulerian coordinates that are comoving with the Hubble flow (i.e. the Lagrangian coordinates of the background solution) are denoted by ${\vec{q}}$. In both cases, Lagrangian coordinates are denoted by ${\vec{X}}$ (note that Lagrangian coordinates ${\vec{q}}$ for the homogeneous solution are regarded as Eulerian coordinates in an inhomogeneous setting). Greek indices refer to particles, while Latin indices refer to vector components; a repeated index indicates summation.

   
2 Dynamical equations through coarse-graining

We consider a system of N identical particles evolving under gravitational forces. In the cosmological context, these particles are the constituents of the dark matter. The purpose is to model the dynamical clustering of these particles, which is supposed to lead to the observed large-scale structure. In the Newtonian approximation, the equations of motion for the position and velocity of the particles in phase space are

$$% {\dot{\vec{x}}}^{(\alpha)} = {\vec{v}}^{(\alpha)}, \qquad {\dot{\vec{v}}}^{(\alpha)} = {\vec{g}}^{(\alpha)},$$ (1)

where the gravitational field strength ${\vec{g}}^{(\alpha)}$ is determined self-consistently by the Newtonian field equations,

 

$$\nabla^{(\alpha)} \cdot {\vec{g}}^{(\alpha)} = \Lambda - 4 \pi G ... ...eta \neq \alpha}^N \delta\left({\vec{x}}^{(\alpha)}-{\vec{x}}^{(\beta)}\right),$$

$$\nabla^{(\alpha)} \times {\vec{g}}^{(\alpha)} = {\vec{0}}.$$ (2)

m is the mass of a particle, G denotes Newton's gravitational constant and $\Lambda$ the cosmological constant.

There is an alternative way of writing these equations. We define the Klimontovich density in the one-particle phase space as:

$$% f_{\rm K} ({\vec{x}}, {\vec{v}}, t) := \sum_{\alpha=1}^N \d... ...\right) \; \delta\left({\vec{v}}-{\vec{v}}^{(\alpha)}\right).$$ (3)

The equations of motion (1-2) lead to a dynamical equation for this density:

$$\frac{\partial f_{\rm K}}{\partial t} + {\vec{v}}\cdot \fra... ...{g}}\cdot \frac{\partial f_{\rm K}}{\partial {\vec{v}}} = 0,$$ (4a)

$$\nabla \cdot {\vec{g}}:= \Lambda - 4 \pi G m \int {\rm d}{\ve... ...x}},{\vec{v}}, t), \; \nabla \times {\vec{g}}:= {\vec{0}}.$$ (4b)

This description of the system is exhaustive and much more detailed than one is usually interested in. A coarser description may be obtained by smoothing the Klimontovich density: let ${\cal L}$, ${\cal V}$ denote respectively the spatial and velocity smoothing scales; then a smooth phase-space density can be defined as:

 

$$% f({\vec{x}}, {\vec{v}}, t) := \int \frac{{\rm d}{\vec{x}}'}{{\c... ...c{v}}-{\vec{v}}'}{{\cal V}} \right) f_K \left({\vec{x}}', {\vec{v}}', t\right),$$ (5)

where $W (\cdot)$ is a rotationally symmetric coarsening window normalized to unity. We can write the evolution equation of this new density according to Eqs. (4):

$$% \frac{\partial f}{\partial t} + {\vec{v}}\cdot \frac{\part... ...} -\frac{\partial}{\partial {\vec{v}}} \cdot {\vec S}^{(g)},$$ (6a)

$$% \nabla \cdot \bar{\vec{g}}:= \Lambda - 4 \pi G m \int {\r... ...{x}},{\vec{v}}, t), \nabla \times \bar{\vec{g}}:= {\vec{0}},$$ (6b)

with

$$% {\vec S}^{(v)} ({\vec{x}}, {\vec{v}}, t) := \int \frac{{\rm d}{... ...vec{v}}' - {\vec{v}}\right) ~ f_{\rm K} \left({\vec{x}}', {\vec{v}}', t\right),$$ (6c)

$$% {\vec S}^{(g)} ({\vec{x}}, {\vec{v}}, t) := \int \frac{{\rm d}{... ...ec{g}}({\vec{x}}, t)\right] ~ f_{\rm K} \left({\vec{x}}', {\vec{v}}', t\right).$$ (6d)

The field $\bar{\vec{g}}({\vec{x}},t)$ defined by Eqs. (6b) is the gravitational mean field associated with the phase-space density  $f({\vec{x}},{\vec{v}},t)$. The terms  ${\vec S}^{(v)}$ and  ${\vec S}^{(g)}$represent the dynamical coupling to the degrees of freedom removed by the smoothing procedure: ${\vec S}^{(v)}$ accounts for the velocity dispersion, while  ${\vec S}^{(g)}$ represents departures from the mean field gravity on the smoothing scales. Notice that Eq. (6a) is, like the exact Eq. (4a), a conservation equation in the one-particle phase space, since the total number of particles $N=\int {\rm d}{\vec{x}}~{\rm d}{\vec{v}}f({\vec{x}},{\vec{v}},t)= \int {\rm d}{\vec{x}}' ~{\rm d}{\vec{v}}' f_{\rm K} ({\vec{x}}',{\vec{v}}',t)$.

Equations (6a,b) do not form a closed system of equations unless the sources  ${\vec S}^{(v)}$, ${\vec S}^{(g)}$ can be represented or approximated as functionals of  $f({\vec{x}},{\vec{v}},t)$. The form and success of the approximation depends in general on the initial conditions and the choice of the coarsening scales ${\cal L}$, ${\cal V}$. A frequently used approximation in many examples of cosmological and astrophysical interest consists of neglecting small-scale departures from the coarse variables, ${\vec S}^{(v)} \approx {\vec{0}}$, ${\vec S}^{(g)} \approx {\vec{0}}$, so that Eqs. (6a,b) reduce to the Vlasov-Poisson system of equations for a smooth  $f({\vec{x}},{\vec{v}},t)$ (e.g., Peebles 1980; Binney & Tremaine 1987).

However, to describe the formation of cosmological structures, the continuous phase-space density  $f({\vec{x}},{\vec{v}},t)$ is still too detailed, since the observables we are interested in are the large-scale density and velocity fields. A mass density and a mean fluid velocity can be defined by the velocity moments of  $f({\vec{x}},{\vec{v}},t)$:

$$% \rho ({\vec{x}}, t) := m \int {\rm d}{\vec{v}}\; f ({\vec{... ...\left(\frac{{\vec{x}}-{\vec{x}}^{(\alpha)}}{{\cal L}} \right),$$ (7a)

$$% \rho \bar{\vec{v}}({\vec{x}}, t) := m \int {\rm d}{\vec{v}}... ...{\vec{x}}^{(\alpha)}}{{\cal L}} \right) {\vec{v}}^{(\alpha)},$$ (7b)

where the last expressions follow from inserting successively the definitions of f, Eq. (5), and $f_{\rm K}$, Eq. (3). Equation (7b) defines the Eulerian mean fluid velocity  $\bar{\vec{v}}$.

The evolution equations for these two fields, expressing mass and momentum conservation, follow immediately from Eqs. (6) or from Eqs. (1)-(2):

$$% \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \bar{\vec{v}}) = 0,$$ (8a)

$$% \frac{\partial \bar{\vec{v}}}{\partial t} + (\bar{\vec{v}}\... ... \frac{1}{\rho} (\vec{\cal F} - \nabla \cdot \vec{\cal P}),$$ (8b)

$$% \nabla \cdot \bar{\vec{g}}:= \Lambda - 4 \pi G \rho, \qquad \nabla \times \bar{\vec{g}}:= {\vec{0}}.$$ (8c)

$\bar{\vec{g}}$ is identical to the gravitational mean field strength defined in Eq. (6b). Two new fields have emerged:

 

$$% {\cal F}_i ({\vec{x}}, t) m \int {\rm d}{\vec{v}}\; S_i^{(g)} ({\vec{x}}, {\vec{v}}, t)$$ :=

$$\frac{m}{{\cal L}^3} \sum_\alpha W \left(\frac{{\vec{x}}-{\vec{x}... ...}}{{\cal L}} \right) \left[ g_i^{(\alpha)} - {\bar g}_i ({\vec{x}}, t) \right],$$ = (9a)

$$\begin{eqnarray*}{\cal P}_{ij} ({\vec{x}}, t) &:= & m \int {\rm d}{\vec{v}}\Big... ...\vec{x}}, t) \right] S_i^{(v)}({\vec{x}}, {\vec{v}}, t) \Bigr\} \end{eqnarray*}$$

$$% \hspace*{1.3cm} = \frac{m}{{\cal L}^3} \sum_\alpha W \lef... ...v_j^{(\alpha)} - {\bar v}_i {\bar v}_j ({\vec{x}}, t)\right].$$ (9b)

These new fields, in analogy to ${\vec S}^{(v)}$ and  ${\vec S}^{(g)}$, represent the coupling of the dynamics of the relevant fields $\rho$and  $\bar{\vec{v}}$ to the configurational and kinetic degrees of freedom that have been averaged out. The symmetric second-rank tensor field  $\vec{\cal P}$ accounts for the velocity dispersion around the mean velocity. As Eq. (9b) shows, there are two contributions when working at the level of the smoothed phase-space density  $f({\vec{x}},{\vec{v}},t)$, which are subsumed in a single term at the description level of the point particles. The vector field  $\vec{\cal F}$ contains the deviation from mean field gravity - alternatively, $\bar{\vec{g}}$ is the gravitational field created by the monopole moment of the mass distribution in the subsystem defined by the window, while  $\vec{\cal F}$ is the contribution from the higher mass-multipoles generated by the small-scale structure.

The set of hydrodynamic-like Eqs. (8) can be derived through the intermediate step of a kinetic equation, Eq. (6a), or directly from the microscopic particle dynamics, Eqs. (1). In the latter case, ${\cal V}$ disappears from the definitions (7, 9) after integrating over the velocities, indicating that the precise value of ${\cal V}$ is unimportant for Eqs. (8).

An advantage of the formal derivation through coarse-graining as presented here is that it makes explicit the smoothing scales ${\cal L}$ and ${\cal V}$, which are usually implicit in most applications and not clearly presented. Although these two scales are arbitrary, their choice is usually dictated by the physics of the problem at hand. Typically ${\cal V}=0$, and we know of only one application where  ${\cal V} \neq 0$: the numerical codes to solve the Vlasov equation necessarily have a finite resolution in phase-space and coarse-graining in velocity is a method to treat the filamentation problem in velocity space induced by this equation (Klimas 1987). (In other works (e.g. Lynden-Bell 1967; Shu 1978) phase-space coarsening is invoked but the precise value of ${\cal V}$, being irrelevant, is not addressed).

The spatial scale ${\cal L}$ is usually chosen in concordance with the length scale of the phenomena of interest. In the regimes of validity of the paradigmatic examples of kinetic theory (Boltzmann equation and Vlasov equation), ${\cal L}$ is much larger than the mean interparticle distance $\ell$, so that  $f({\vec{x}},{\vec{v}},t)$ does not resemble the spiky Klimontovich density  $f_{\rm K}({\vec{x}},{\vec{v}},t)$ (Spohn 1991, and refs. therein). The Boltzmann equation for dilute gases with a short range $\sigma$ of interaction follows from Eq. (6a) when the mean field force  $\bar{\vec{g}}$ is negligible compared to the effect of the term  ${\vec S}^{(g)}$, which is dominated by binary "close encounters'' of uncorrelated particles. This holds in the scaling limit $\sigma \rightarrow 0$ with $\ell \sim \sigma^{2/3}$ (dilute limit) and ${\cal L} \sim \sigma^{0}$ (continuum limit), and the mean free path, $\lambda \sim \ell^3/\sigma^2$, remains finite. The Vlasov equation describes the dynamical evolution of  $f({\vec{x}},{\vec{v}},t)$ when the interaction is long-ranged and weak. For the gravitational interaction, Eqs. (2), this holds in the scaling limit $m \rightarrow 0$ with $\ell \sim m^{1/3}$ (finite mass density) and ${\cal L} \sim m^{0}$ (continuum limit), so that the particle distribution is statistically homogeneous on scales below ${\cal L}$ and the mean field  $\bar{\vec{g}}$ dominates over the effect of "close encounters'' contained in  ${\vec S}^{(g)}$.

Ma & Bertschinger (2004) study the Klimontovich equation for particles moving in a given stochastic large-scale gravitational field, intended to be a model of particle evolution in a galaxy halo environment; after averaging over realizations of the large-scale gravitational field a new kinetic equation for the ensemble averaged f is obtained that differs from the Vlasov equation (this ensemble averaging is, in the language of statistical physics, similar to an average over "disorder'' rather than an average over "thermal noise''). For comparison, our Eqs. (6) describe the evolution of a single realization of a coarse-grained distribution (note that the mathematical structure of both averaging approaches are similar with reinterpretation of the window function $W (\cdot)$ as a probability density). Of course, in addition to coarse-graining one could implement ensemble averaging. The addition of a noise term as a model of small-scale degrees of freedom involves interesting physics (e.g., Berera & Fang 1994; Barbero et al. 1997; Domínguez et al. 1999; Buchert et al. 1999; Matarrese & Mohayaee 2002; Antonov 2004).

   
3 Dynamical equations for cosmological structure formation

In this section we apply the method of coarse-graining to the case of structure formation in cosmology by gravitational instability. In this case it is convenient to introduce "comoving'' Eulerian coordinates attached to a homogeneous-isotropic solution (Friedmann-Lemaître backgrounds), characterized by the expansion factor a(t) (and Hubble's function  $H=\dot{a} / a$), with tiny inhomogeneities superimposed as the seeds of structure formation. In order to subtract the homogeneous-isotropic motion, one defines comoving positions  ${\vec{q}}^{(\alpha)}$, peculiar-velocities  ${\vec{u}}^{(\alpha)}$ and gravitational peculiar-accelerations  ${\vec{w}}^{(\alpha)}$ as follows:

 

$${\vec{q}}^{(\alpha)} := \frac{1}{a} {\vec{x}}^{(\alpha)},$$

$${\vec{u}}^{(\alpha)} := {\vec{v}}^{(\alpha)} - H {\vec{x}}^{(\alpha)},$$

$${\vec{w}}^{(\alpha)} := {\vec{g}}^{(\alpha)} + \frac{4 \pi G \varrho_{\rm H} - \Lambda}{3} {\vec{x}}^{(\alpha)},$$ (10)

where $\varrho_{\rm H}(t)$ is the homogeneous background density, satisfying $\dot{\varrho}_{\rm H} + 3 H \varrho_{\rm H} =0$. Friedmann's equation $3 \frac{\ddot a}{a} = \Lambda - 4\pi G \varrho_{\rm H}$ has been used in Eq. (10). In terms of these variables, Eqs. (1) and (2) are:

 

$${\dot{\vec{q}}}^{(\alpha)} = \frac{1}{a} {\vec{u}}^{(\alpha)},\qquad {\dot{\vec{u}}}^{(\alpha)} = {\vec{w}}^{(\alpha)} - H {\vec{u}}^{(\alpha)},$$

$$\nabla_{{\vec{q}}}^{(\alpha)} \cdot {\vec{w}}^{(\alpha)} = - 4 \p... ...\delta\left({\vec{q}}^{(\alpha)}-{\vec{q}}^{(\beta)}\right) -\varrho_H \right],$$

$$\nabla_{{\vec{q}}}^{(\alpha)} \times {\vec{w}}^{(\alpha)} = {\vec{0}}.$$ (11)

We can now repeat the derivations of Sect. 2. We define a mass density and a mean peculiar-velocity as follows:

 

$$\varrho ({\vec{q}}, t) := \frac{m}{(a L)^3} \sum_\alpha W \left(\frac{{\vec{q}}-{\vec{q}}^{(\alpha)}}{L} \right),$$ (12a)

$$\varrho \bar{\vec{u}}({\vec{q}}, t) := \frac{m}{(a L)^3} \sum_\alpha W \left(\frac{{\vec{q}}-{\vec{q}}^{(\alpha)}}{L} \right) {\vec{u}}^{(\alpha)},$$ (12b)

where $L={\cal L}/a$ is the (comoving) coarsening length. This length is in principle arbitrary, its use being motivated by the features of cosmological structure formation one is interested in. In general, one should demand that L will be much larger than the (comoving) mean interparticle distance, since dark matter discreteness is cosmologically irrelevant. In the rest of the paper, when we speak about "small/large scales'', we use L to set the comparison scale.

The equations obeyed by these fields can be derived as before. Note that, in the following we drop the overbar used to distinguish the microscopic velocities and accelerations from the mean velocities and mean field strengths; we always refer to the following equations, so that we no longer need this distinction (the time-derivative is taken at constant ${\vec{q}}$ and L in these equations; see Appendix A for the relation to Eqs. (8)):

$$% \frac{\partial \varrho}{\partial t} + 3 H \varrho + \frac{1}{a} \nabla_{\vec{q}}\cdot (\varrho {\vec{u}}) = 0,$$ (13a)

$$% \frac{\partial {\vec{u}}}{\partial t} + \frac{1}{a} \left({... ...{F}}- \frac{1}{a} \nabla_{\vec{q}}\cdot {\vec{\Pi}}\right),$$ (13b)

$$% \nabla_{\vec{q}}\cdot {\vec{w}}:= - 4 \pi G a \left(\varrho - \varrho_{\rm H}\right),$$ (13c)

$$% \nabla_{\vec{q}}\times {\vec{w}}:= {\vec{0}}.$$ (13d)

Here, ${\vec{w}}$ (with overbar omitted) is the mean field strength of the peculiar-gravity, and we have the "peculiar counterparts'' to the fields (9):

$$% F_i ({\vec{q}}, t) := \frac{m}{(a L)^3} \sum_\alpha W \le... ...\right) \left[ w_i^{(\alpha)} - {w}_i ({\vec{q}}, t)\right],$$ (14a)

$$% \Pi_{ij} ({\vec{q}}, t) := \frac{m}{(a L)^3} \sum_\alpha W ... ...\alpha)} u_j^{(\alpha)} - {u}_i {u}_j ({\vec{q}}, t) \right].$$ (14b)

If one computes dynamical equations for the new fields ${\vec{F}}$ and  ${\vec{\Pi}}$, one obtains further new fields, and so on. In order for Eqs. (13) to become useful, one has to close this infinite hierarchy by resorting to physically motivated approximations that yield the fields ${\vec{F}}$ and  ${\vec{\Pi}}$ as functionals of $\varrho$ and ${\vec{u}}$. That a successful closure exists is not guaranteed and this usually depends on the choice of the scale L and the initial and boundary conditions. In the context of fluids dominated by short-range interactions, there is a set of closed hydrodynamic equations when the fields $\varrho$ and ${\vec{u}}$ vary spatially only over length scales much larger than any microscopic length (mean free path, correlation length, etc.). In such cases, one must take L to lie between these microscopic lengths and the scale of variation of the fields. Closure is achieved with the assumption of local quasi-equilibrium: the coarsening cells of size $\sim$L are idealized as systems in almost thermal equilibrium, the corrections to equilibrium being proportional to the relatively small spatial gradients of the fields. In this manner one obtains the Navier-Stokes and Fourier equations. The existence of a hydrodynamic regime where these equations hold can be supported theoretically for dilute gases described by the Boltzmann equation by means of the Chapman-Enskog expansion (Chapman & Cowling 1991); for dense fluids, the hydrodynamic regime has the status of an experimental fact.

In the context of cosmological structure formation one cannot follow this argument, since the notion of thermal equilibrium is not well-defined. We instead discuss three closures:

1. The "dust model'': one assumes that the large-scale dynamics is dominated by the large-scale forces and neglects the coupling to the small scales altogether: one sets $F_i \equiv 0$, $\Pi_{ij,j} \equiv 0$ in Eqs. (13), and obtains the Euler-Newton system for the peculiar-fields:

$$% \frac{\partial \varrho}{\partial t} + 3 H \varrho + \frac{1}{a} \nabla_{\vec{q}}\cdot (\varrho {\vec{u}}) = 0,$$ (15a)

$$% \frac{\partial {\vec{u}}}{\partial t} + \frac{1}{a} ({\vec{u}}\cdot \nabla_{\vec{q}}) {\vec{u}}+ H {\vec{u}}= {\vec{w}},$$ (15b)

$$% \nabla_{\vec{q}}\cdot {\vec{w}}= - 4 \pi G a (\varrho - \varrho_{\rm H}),$$ (15c)

$$% \nabla_{\vec{q}}\times {\vec{w}}= {\vec{0}}.$$ (15d)

This system of equations has a long history and has been thoroughly studied in the cosmological context (see, e.g., Peebles 1980; Sahni & Coles 1995; Ehlers & Buchert 1997; Bernardeau et al. 2002). Coming from the context of hydrodynamics, the matter model "dust'' is also called a "pressure-less fluid''; from our reasoning it is clear that this model is characterized by the assumptions that (i) small-scale inhomogeneities are irrelevant, so that mean-field gravity is dominant, $F_i \approx 0$, and (ii) multi-streaming is absent and small-scale kinetic degrees of freedom are subdominant too, so that "pressure'' drops, $\Pi_{ij} \approx 0$.

The "dust model'' was originally applied to the "top-down scenario'' of structure formation, e.g., "Hot Dark Matter cosmogony'', in which large-scale density inhomogeneities grow while the small scales remain homogeneous. The model, however, turned out to be a relatively good description also for the "bottom-up scenario'', e.g., "Cold Dark Matter cosmogony'' (Pauls & Melott 1995), in which the small scales are always strongly inhomogeneous. Strictly speaking, the "dust model'' is not properly defined if there is too much initial small-scale power, because of the emergence of singularities at arbitrarily short times, so that a small-scale cutoff is required (this becomes manifest in formal perturbative expansions too, Valageas 2002; Bernardeau et al. 2002). Our derivation of the "dust model'' clearly demonstrates that the smoothing length L is indeed a defining ingredient of the model. However, this fact has not been properly emphasized in the literature - the smoothing length is likely irrelevant in the "top-down scenario'', where the initial conditions have a built-in length of smoothness (free-flight scale of the dark matter particles), but this is not so in a "bottom-up scenario'': in this latter case, it was found that a much better agreement with N-body simulations is achieved if the initial conditions are first smoothed, e.g., the "Truncated Zel'dovich Approximation'' (Coles et al. 1993).

The "dust model'' correctly describes many features of the formation of structure by gravitational instability. It has also an important shortcoming, which is the focus of the present work, namely that it generates caustics, i.e., density singularities, as well as multi-streaming regions (shell-crossing), where the velocity field is multi-valued. This points to a breakdown of the approximation, so that the term  $F_i-a^{-1} \Pi_{ij,j}$ in Eq. (13b) is no longer negligible. In the present work we address two approximations beyond the "dust model'' in which this term is modelled as a function of $\varrho$ and ${\vec{u}}$.

2. The Euler-Jeans-Newton ("EJN'') model: motivated by hydrodynamics and the theory of stellar systems, we have the following phenomenological closure: the corrections to mean field gravity are neglected, and the velocity dispersion is approximated by an isotropic tensor field which, in addition, is commonly modelled as a function of the local density:

$$% F_i = 0\;, \qquad \Pi_{ij} = p(\varrho) \delta_{ij},$$ (16)

with a given function $p(\varrho) > 0$, a kinetic pressure. (Appendix C discusses some examples of the functional dependence  $p(\varrho)$.) This velocity dispersion can be a consequence of the multi-streaming or it can be already present in the initial conditions. This model is applicable to situations in which deviations from the mean field can be ignored. For cosmological structure formation in a "top-down evolution'', it would be safe to ignore deviations from mean field, since then the scales below the collapsing ones are relatively homogeneous (e.g., due to the free-flight of the dark matter particles as in a "Hot-Dark-Matter cosmogony''), and the main correction to "dust'' is due to velocity dispersion. On the other hand, if a scenario that entails "bottom-up'' structure formation is investigated realistically, especially down to galaxy halo scales, then the scales below the collapsing ones are very inhomogeneous due to the previously formed clusters, and both velocity dispersion and stresses due to deviations from mean field gravity can be expected to be relevant.

Given the phenomenological character of the approximation (16), one can think of the relationship  $p(\varrho)$ as "templates'' modeling the overall effect of both velocity dispersion and departure from mean field gravity as long as the mathematical analysis of Eqs. (17) is involved. The only exact constraint required by this interpretation is that ${\vec{F}}$ can be written as the divergence of a stress tensor. This is the case in the other closure ansatz to be discussed (see Eq. (20)).

Inserting the ansatz for F_i and $\Pi_{ij}$ into Eq. (13b), one gets ( ${\Delta}_{\vec{q}}$ is the Laplacian operator in the variable ${\vec{q}}$):

$$% \frac{\partial \varrho}{\partial t} + 3 H \varrho + \frac{1}{a} \nabla_{\vec{q}}\cdot (\varrho {\vec{u}}) = 0,$$ (17a)

$$% \frac{\partial {\vec{u}}}{\partial t} + \frac{1}{a} ({\vec{u}}\cdot \nabla_{\vec{q}}) {\vec{u}}+ H {\vec{u}} {\vec{w}}- \frac{p'(\varrho)}{a \varrho} \nabla_{\vec{q}}\varrho$$ =

$${\vec{w}}+ \frac{\lambda^2 (\varrho)}{a^2} {\Delta}_{\vec{q}}{\vec{w}},$$ = (17b)

$$% \nabla_{\vec{q}}\cdot {\vec{w}}= - 4 \pi G a \left(\varrho - \varrho_{\rm H}\right),$$ (17c)

$$% \nabla_{\vec{q}}\times {\vec{w}}= {\vec{0}},$$ (17d)

where we have defined a density-dependent length (see Eq. (34) for the relation to Jeans' length):

$$% \lambda (\varrho) := \left[ \frac{p'(\varrho)}{4 \pi G \varrho} \right]^{1/2}\cdot$$ (18)

Equations (17a,b) are Euler's equations for a compressible fluid, Eq. (17b) is known as Jeans' equation in stellar systems theory. This model was proposed by Buchert & Domínguez (1998) starting from Vlasov's equation for the one-particle distribution function, so that corrections to mean field gravity were neglected from the outset and only velocity dispersion remained. In that work it was also suggested that the phenomenological ansatz Eq. (16) could be replaced by a different closure condition leading to $p \propto \varrho^{5/3}$ (see Appendix C).

The case of isotropic stresses formally covers the hydrodynamics of a perfect fluid or, as common in studies of stellar systems, polytropic models. However, there is no compelling reason for isotropic stresses in collisionless systems before the onset of virialization. Even "virialized'' systems will in general maintain an anisotropic component. At the moment, we consider the isotropy assumption as a good working hypothesis to understand the role of interaction between multi-stream forces and self-gravity: the velocity dispersion ellipsoid is approximated by a sphere. (Interesting considerations of the influence of the anisotropic part have been reported by Maartens et al. 1999; see also Barrow & Maartens 1999.)

3. The Small-Size Expansion ("SSE''): this is a method proposed by Domínguez (2000, 2002) to formalize the notion that the coupling to the small scales is in some sense weak. One argues that the corrections F_i and  $\Pi_{ij,j}$ are determined mainly by the largest scales contributing to them, i.e., those close to L. For instance, in the bottom-up scenario, the large-scale dynamics is sensitive mainly to the motion of the most recently formed clusters as a whole, and not so much to their internal structural details due to the trapped particles, so that one takes $L \sim$ typical size of these clusters (which play the role of effective particles of size $\sim$L). The mathematical implementation of this idea leads to a formal expansion of F_i and $\Pi_{ij}$ in powers of the smoothing length L (see Appendix B for an outline of the derivation):

$$% F_i \rightarrow B L^2 (\partial_k \varrho) (\partial_k w_i) + {\cal O}\left(L^4\right),$$ (19a)

$$% \Pi_{ij} \rightarrow B L^2 \varrho (\partial_k u_i) (\partial_k u_j) + {\cal O}\left(L^4\right),$$ (19b)

where B>0 is a dimensionless constant of order unity depending on the shape of the smoothing window $W (\cdot)$, Eq. (B.5). The correction to mean field, being proportional to  $\partial_k w_i$, is approximated by the tidal forces induced by the mean field gravity. With Eqs. (13c,d), it can be written also as the divergence of a second-rank tensor:

$$% F_i \rightarrow B L^2 \partial_j \left(\varrho ~ \partial_... ... G a \varrho^2 \delta_{ij}\right) + {\cal O}\left(L^4\right).$$ (20)

The velocity dispersion tensor is represented through the kinematical parts of the peculiar-velocity gradient: expansion rate $\theta:=(1/a) \nabla_{\vec{q}}\cdot {\vec{u}}$, vorticity $\vec{\omega}:=(1/a) \nabla_{\vec{q}}\times {\vec{u}}$, and shear  $\sigma_{ij}$:

$$% \frac{1}{a} \partial_i u_j = \frac{1}{3} \theta \delta_{ij} + \sigma_{ij} + \frac{1}{2} \varepsilon_{ijk} \omega_{k},$$ (21)

$$% \Pi_{ij} \rightarrow a^2 B L^2 \varrho ~ \biggl[ \frac{1}{9} \theta^2 \delta_{ij} + \frac{2}{3} \theta \sigma_{ij} + \sigma_{ik} \sigma_{kj}$$

$$- \frac{1}{2}\left(\varepsilon_{ikn} \sigma_{kj} + \varepsilon_{jkn} \sigma_{ki}\right)\omega_n$$

$$+ \frac{1}{4} \left(\omega^2 \delta_{ij} - \omega_i \omega_j\right) \biggr] + {\cal O}\left(L^4\right).$$ (22)

In the context of the ``SSE model'', we recover the ``dust model'' to order ${\cal O}(L^0)$. To order ${\cal O}(L^2)$, we obtain the following equations:

$$% \frac{\partial \varrho}{\partial t} + 3 H \varrho + \frac{1}{a} \nabla_{\vec{q}}\cdot (\varrho {\vec{u}}) = 0,$$ (23a)

$$% \frac{\partial {\vec{u}}}{\partial t} \frac{1}{a} ({\vec{u}}\cdot \nabla_{\vec{q}}) {\vec{u}}+ H {\vec{u}}= {\vec{w}}$$ +

$$+ \frac{B L^2}{\varrho} \left\{ \left(\nabla_{\vec{q}}\varrho \cd... ...ot \left[\varrho (\partial_k {\vec{u}}) (\partial_k {\vec{u}})\right] \right\},$$ (23b)

$$% \nabla_{\vec{q}}\cdot {\vec{w}}= - 4 \pi G a \left(\varrho - \varrho_{\rm H}\right),$$ (23c)

$$% \nabla_{\vec{q}}\times {\vec{w}}= {\vec{0}}.$$ (23d)

The expansions (19) are formally exact. The assumption underlying the "SSE'' is that the first few terms provide an accurate description of the dynamical effect of the mode-mode coupling for the problem at hand. For comparison, this is not true in the hydrodynamic regime of fluids dominated by short-ranged interactions, where the high-order terms in the expansions have to add up to yield the usual L-independent hydrostatic pressure and viscous force.

The closure assumption (19b) for the velocity dispersion has been also used to model the influence of subresolution degrees of freedom of the velocity field in "Large-Eddy Simulations'' of turbulent flow (Pope 2000). It has been termed the "gradient model'' (see, e.g., Vreman et al. 1997) and it was introduced as part of the "Clark model'' (Clark et al. 1979) (a model which includes an extra additive term in the expression for $\Pi_{ij}$).

Both the "EJN'' and the "SSE'' models provide an extension of the "dust model'' with the potential to improve on it by preventing the formation of singularities. The "EJN model'' relies on a phenomenological assumption concerning the kinetic stress $\Pi_{ij}$, but has the advantage that the resulting equations are well established. In favor of the "SSE'' model is the fact that the corrections to "dust'' can be obtained systematically beyond phenomenology, but the mathematical and physical status of the resulting equations is little explored yet.

   
3.1 Where do these equations hide the "adhesion approximation''?

In this subsection we get a first picture of the effect of the corrections to the "dust model''. In the "adhesion model'' (Gurbatov et al. 1989), the evolution of "dust'' is modelled with Zel'dovich's approximation (Zel'dovich 1970, 1973), a powerful and mathematically manageable description of the exact evolution (first-order Lagrangian perturbation solution, see Sect. 5.1). Accordingly, one assumes in Eq. (15b) the proportionality of peculiar-velocity and -acceleration fields (cf. Peebles 1980; Buchert 1989, 1992; Bildhauer & Buchert 1991; Kofman 1991; Buchert 1993b; Vergassola et al. 1994):

$$% {\vec{w}}= F(t) {\vec{u}}, \qquad F(t) = 4\pi G\varrho_{\rm H} \frac{b(t)}{{\dot b}(t)},$$ (24)

where b(t) is identical to the growing-mode solution of the Eulerian linear theory of gravitational instability for "dust'' (i.e., it solves the equation $\ddot{b} + 2 H \dot{b} - 4 \pi G \varrho_{\rm H} b = 0$). Notice that this "slaving'' assumption implies an irrotational flow: $\nabla_{{\vec{q}}} \times {\vec{u}}={\vec{0}}$. Changing the temporal variable from t to b and defining a rescaled peculiar-velocity field $\hat{\vec{u}}= {\vec{u}}/ a \dot{b}$, Eq. (15b) becomes, after inserting Eq. (24),

$$% \frac{{\rm d} \hat{\vec{u}}}{{\rm d} b} = {\vec{0}}, \qquad... ...}{\partial b} + \hat {\vec{u}}\cdot \nabla_{\vec{q}}\right).$$ (25)

The evolution according to this equation leads (kinematically) to the formation of singularities, where $\varrho \rightarrow +\infty$, and the velocity field becomes multi-valued. Following the idea of the "adhesion model'', the corrections to "dust'' are assumed to be negligible almost everywhere, i.e., except at the localized singularities. The corrections are then estimated with the "dust'' solution (i.e., with Zel'dovich's approximation):

-

The "EJN model'': inserting the parallelism condition (24) into Eq. (17b), and performing the same change of variables as before, we obtain, instead of Eq. (25), the following corrected equation:

$$% \frac{{\rm d} \hat{\vec{u}}}{{\rm d} b} = {\mu} {\Delta}_{\... ... \left(\nabla_{\vec{q}}\times \hat{\vec{u}}= {\vec{0}}\right),$$ (26a)

with the GM or gravitational multi-stream coefficient (Buchert & Domínguez 1998; Buchert et al. 1999)

$$% \mu := \frac{F}{a^2 {\dot b}} \lambda^2 (\varrho) > 0,$$ (26b)

formally playing the role of a "viscosity''.

-

The "SSE model'': the parallelism condition allows us to express ${\vec{w}}$ and $\varrho$ (via Eq. (23c)) in terms of ${\vec{u}}$. The correction term is hypothesized to be relevant only near places where Zel'dovich's approximation yields singularities, where $\vert\nabla {\vec{u}}\vert \to \infty$, so that it can be approximated by neglecting subdominant powers of  $\nabla {\vec{u}}$:

$$% \frac{1}{\varrho} \left[ F_i - \frac{1}{a} \Pi_{ij,j} \righ... ...t(\partial_k u_i\right) \left(\partial_k u_j\right) \right].$$ (27)

To simplify this expression further, one recalls that, within Zel'dovich's approximation, density singularities arise generically due to a locally plane-parallel collapse and the velocity gradient looks locally like $\partial_i u_j \approx (\nabla_{\vec{q}}\cdot {\vec{u}}) n_i n_j$, where the unit vector $\vec{n}$ defines the plane of local collapse. In order to compute the correction term, we assume that the spatial dependence of $\vec{n}$ is negligible compared to the (large) gradient of the fields along $\vec{n}$, i.e., as if it were a globally plane-parallel collapse, and we obtain:

$$\frac{1}{\varrho} \left[ F_i - \frac{1}{a} \Pi_{ij,j} \right]... ...t(\nabla_{\vec{q}}\cdot {\vec{u}}\right) \Delta_{\vec{q}}u_i.$$

Therefore, Eq. (23b) simplifies, after the change of variables, to the same Eq. (26a) with a different GM coefficient:

$$% \mu := \frac{3 B L^2}{b \varrho_{\rm H}} \varrho > 0.$$ (28)

The dominant contribution in Eq. (27) is the one by velocity dispersion. It can be easily checked that the correction to the mean field gives, under the same assumptions, a subdominant but also positive contribution to the GM coefficient.

Equation (26) resembles the well-known key equation of the original "adhesion approximation'' (Gurbatov et al. 1989), to which it reduces when $\mu$ = constant  $\rightarrow 0^+$. For constant $\mu$ it can be solved analytically for curl-free flows (Hopf-Cole transformation of the 3D Burgers' equation); this solution shows that the singularities predicted by the "Zel'dovich approximation'', Eq. (25), are indeed regularized by the $\mu$-term (see, e.g., Vergassola et al. 1994); for an interesting application of Burgers' equation and further insight see also the model of Jones (1996) for a two-component system. In the more general cases of a density-dependent GM coefficient, no analytical solution was found. Nevertheless, the application of boundary-layer theory (Domínguez 2000) shows that density singularities are still regularized (provided $p(\varrho \rightarrow +\infty) \sim \varrho^\gamma$ with $\gamma>1$ in the case of the "EJN model'', so that pressure can oppose gravity successfully). The solutions for the velocity and density fields are qualitatively identical to those of the original "adhesion model'', developing into a shock structure of infinite density in the limit $\mu \rightarrow 0^+$.

This derivation of the "adhesion approximation'' offers insight into the physics involved. The $\mu$-term, which was phenomenologically motivated by analogy with viscosity, is not related to a truly dissipative process. From a mathematical point of view, the starting point ("EJN model'' (17) or "SSE model'' (23)), the approximations, and the final Eq. (26) are formally time-reversible, i.e., invariant under the transformation $t \rightarrow -t$, ${\vec{u}}\rightarrow -{\vec{u}}$. The correction to "dust'' can transform mean kinetic energy ($\propto$ $\varrho \vert{\vec{u}}\vert^2$) into internal kinetic energy ($\propto$$\Pi_{ii}$) and conversely (also into internal gravitational potential energy via F_i, but, as we have shown, the "SSE model'' predicts this contribution to be subdominant near density singularities in the present approximation). The tendency to compression in collapsing regions favors the correction to behave as a drain of mean kinetic energy, mimicking viscous dissipation; it is easy to check that the correction supplies energy in expanding regions (see Sect. 6.4), whose expansion is thus accelerated. This mechanism seems to be universal: this is the reason that "adhesiveness'' arises in the two different models we study, "EJN'' and "SSE'', and that it can be modelled in a qualitatively correct manner by the simplified "adhesion model'' $\mu$ = constant. The GM coefficient depends on time and on density, and this leads to differences between models concerning the inner structure of the regularized density singularities. A density independent GM coefficient (as in the original "adhesion approximation'') is obtained with the imposed relationship $p \propto \varrho^2$ corresponding to the naive application of the virial condition, while the "EJN'' and the "SSE'' models coincide when $p \propto \varrho^3$.

3.2 Relaxing key assumptions

Summarizing, to derive Eq. (26) one makes two important assumptions:

1.

The mean peculiar-gravitational field ${\vec{w}}$ for the evolution of "dust'' is modelled by a first-order Lagrangian perturbation solution which includes parallelism of peculiar-velocity and -acceleration (Zel'dovich's approximation).

2.

The so-approximated evolution of "dust'' is used to estimate the correction to "dust'', which for consistency must be assumed to be negligible everywhere except where singularities would develop. (For the "SSE model'', one makes explicit use of the additional constraint of a locally plane-parallel singularity, which is nevertheless a consequence of Zeldovich's approximation and hence not expected to introduce new features.)

Relaxing (1) means using a better approximation to the exact evolution of "dust''. Zel'dovich's approximation, originally designed as an extrapolation of the standard linear theory of gravitational instability into the weakly non-linear regime, features "slaving'' of peculiar-velocity and -acceleration. This simple, but physically well-motived condition, if assumed initially, is maintained up to second order in Lagrangian perturbation theory (see: Buchert 1989, 1992 for the general first-order solution, Buchert & Ehlers 1993 for explicit discussions within the general second-order solution, and Buchert 1994 for the corresponding third-order solution). Although corrections to the "Zel'dovich approximation'' in general break the parallelism of ${\vec{w}}$ and ${\vec{u}}$, the departures are not expected to be large in the perturbative regime, given that the "Zel'dovich approximation'' performs well before shell-crossing. Nevertheless, using better approximations for the evolution of "dust'' will not prevent singularities from appearing, because this is a general property of the "dust model'' (see Sect. 6).

While relaxing (1) implies construction of sophisticated non-perturbative approximations (since only then can we expect to substantially improve on the performance of Lagrangian perturbation schemes, optimized to match N-body simulations for any kind of dark matter), assumption (2) appears more like an unnecessary simplification. It is crucial for the derivation of equations similar to the key equation of the standard "adhesion approximation'', but it must be relaxed if one intends to improve over the "adhesion model'' in order to describe the subsequent dynamical evolution inside the collapsing structures. Assumption (2) may be rephrased by saying that the evolution is dominated by convection, $({\hat {\vec{u}}}\cdot \nabla_{\vec{q}}) {\hat {\vec{u}}}$, over "viscosity'', $\mu \Delta_{\vec{q}}{\hat {\vec{u}}}$ - in the hydrodynamic jargon, this is the limit "Reynolds number  $\rightarrow +\infty$''. Relaxing this assumption, i.e. finite Reynolds number, means taking account of the back-reaction of the correction to "dust'' on the trajectories of fluid elements away from singularities too. In general, this will also imply corrections to the parallelism of ${\vec{w}}$ and ${\vec{u}}$ as well as those due to the exact "dust'' evolution (see, e.g., Eq. (35)).

Menci (2002) has studied the correction to the original "adhesion approximation'' when assumption (1) is relaxed. The starting point is the set of Eqs. (13) with the correction to "dust'' already modelled as $\mu {\Delta}_{\vec{q}}{\vec{u}}$, $\mu=$ constant $\to 0$. The correction to parallelism, ${\vec{w}}- F {\vec{u}}$, is estimated perturbatively in the limit of small velocities or times (the initial condition is taken to satisfy parallelism) - the procedure employed by Menci is not an expansion in the "viscosity'' $\mu$, although it might formally appear so. The results agree somewhat better with N-body simulations concerning the details of the forming structures at small scales.

The solution to the "adhesion model'', Eq. (26a) with constant $\mu$, cannot be computed as a regular expansion in $\mu$, because the naive zeroth-order term, given by Eq. (25), is undefined after shell-crossing. The exact analytical solution to this model demonstrates that the limit $\mu \rightarrow 0$ yields a weak solution of Eq. (25), exhibiting discontinuities in ${\vec{u}}$. In the pedagogical review by Frisch & Bec (2001) it is illustrated how much the Lagrangian solution of a "multi-streamed dust model'' differs from the Lagrangian solution to the "adhesion model'' after shell-crossing. In particular, it is explained how the emergence of multi-streaming cannot be captured by a Taylor expansion in the time variable of the solution to Eq. (25). As a matter of fact, the corresponding perturbative expansion of the solution to the Vlasov-Poisson system is identical to the one of the "dust model'' when the initial condition is exactly single-streamed (Valageas 2001). This problem can be overridden if a small amount of multi-streaming is allowed in the initial conditions, e.g. in the form of velocity dispersion.

Ribeiro & Peixoto de Faria (2005) also attempted to provide a physical motivation for the "adhesion model'': starting from the velocity potential $\varphi$ (with ${\vec{u}}= \nabla_{\vec{q}}\varphi / a$), it is assumed that the complex variable $\psi=\sqrt{\varrho} \exp{(i \varphi/\nu)}$ obeys a Schrödinger equation, where $\nu$ is a constant with the appropriate dimension. From this assumption, it is found that ${\vec{u}}$ satisfies Eq. (13b) with a correction to "dust'' that can be written as a functional of the density. It is not obvious that with this correction term one can recover the "adhesion model'': it cannot be brought into a simple form $\propto \Delta_{\vec{q}}{\vec{u}}$ after using the parallelism approximation, and the reasoning offered by the authors is wrong from the outset because their Eq. (31) is algebraically false.

In the following sections we study the "EJN'' and the "SSE'' models beyond assumptions (1, 2). First, we consider the standard Eulerian and Lagrangian perturbative expansions, both closely related to the expansion in time just mentioned. One of the key points in the derivation of the "adhesion model'' is the parallelism relation between ${\vec{w}}$ and ${\vec{u}}$, Eq. (24). This relation holds at the (Eulerian as well as Lagrangian) linear order of the "dust model''. Application of the perturbative techniques to the "EJN'' and "SSE'' models will show how this relation is modified by the correction to "dust''. Later we consider non-perturbative approaches: this is a mathematically difficult task, and we are only able to collect some results and provide some hints for future work.

   
4 The Eulerian perturbative expansion

The expansion parameter of the Eulerian perturbative expansion is the amplitude of the departures of the initial conditions from homogeneity. Formally, one writes:

$$% \varrho = \sum_{n=0}^\infty \epsilon^n \varrho^{(n)}, \qquad {\vec{u}}= \sum_{n=0}^\infty \epsilon^n {\vec{u}}^{(n)},$$ (29)

where $\epsilon$ is an expansion parameter, which can be taken to be the variance of the initial density fluctuations, $\epsilon = \sigma(t_0)$, with (note that $\sigma$ is finite because the density field has been smoothed over the scale L, Eq. (12a))

$$% \sigma^2 (t) := \frac{\left\langle \left[\varrho({\vec{0}},... ...{\rm H}(t)\right]^2 \right\rangle}{\varrho^2_{\rm H} (t)}\cdot$$ (30)

Here, $\langle \cdots \rangle$ denotes ensemble average over the initial conditions. Inserting these expansions into Eqs. (13) and equating powers of $\epsilon$, one obtains a set of linear equations for $\varrho^{(n)}({\vec{q}},t)$, ${\vec{u}}^{(n)}({\vec{q}},t)$. The zeroth order is the homogeneous background, $\varrho^{(0)} ({\vec{q}}, t) = \varrho_{\rm H} (t)$, ${\vec{u}}^{(0)} ({\vec{q}}, t) = {\vec{0}}$. Due to the gravitational instability, the fields $\varrho^{(n)}$, ${\vec{u}}^{(n)}$ ($n \geq 1$) eventually grow in time and the expansions (29) (likely to be understood as asymptotic expansions) cease to be a useful representation of the solution after some time T, which can be roughly estimated by the condition  $\sigma(T)=1$.

   
4.1 The linear regime, ${{\cal O}(\epsilon )}$

To this order, one obtains the following equations:

$$% \frac{\partial \varrho^{(1)}}{\partial t} + 3 H \varrho^{(1... ... \varrho_{\rm H} \nabla_{\vec{q}}\cdot {{\vec{u}}}^{(1)} = 0,$$ (31a)

$$% \frac{\partial {{\vec{u}}}^{(1)}}{\partial t} + H {\vec{u}}... ... \frac{1}{a} \nabla_{\vec{q}}\cdot {\vec{\Pi}}^{(1)} \right),$$ (31b)

$$% \nabla_{\vec{q}}\cdot {{\vec{w}}}^{(1)} = - 4 \pi G a \left(\varrho^{(1)} - \varrho_{\rm H}\right),$$ (31c)

$$% \nabla_{\vec{q}}\times {{\vec{w}}}^{(1)} = {\vec{0}}.$$ (31d)

For the "dust'' and the "SSE'' models, ${\vec{F}}^{(1)} = {\vec{0}}$ and ${\vec{\Pi}}^{(1)}=0$, and they do not differ at this order. One can show that the vorticity, $\vec{\omega} = \frac{1}{a}\nabla \times {\vec{u}}$, decays with time, and that the peculiar-velocity becomes parallel to the peculiar-acceleration asymptotically in time,

$$% {\vec{w}}^{(1)} = F(t) {\vec{u}}^{(1)}.$$ (32)

For the "EJN model'', F_i⁽¹⁾ = 0, but:

$$% \Pi_{ij}^{(1)} = p'(\varrho_{\rm H}) ~ \varrho^{(1)} ~ \delta_{ij}.$$ (33)

Then, Eqs. (31) can be simplified to a single equation for the density contrast $\delta := \varrho/\varrho_{\rm H} - 1$:

$$% \frac{\partial^2 \delta^{(1)}}{\partial t^2} + 2 H \frac{\p... ...lambda^2 (\varrho_{\rm H})}{a^2} \Delta_{\vec{q}}\delta^{(1)},$$ (34)

where we have used the length defined in Eq. (18). Notice that in the absence of cosmological expansion, the time-dependence of $\lambda$ induced by  $\varrho_{\rm H}$ drops and it reduces to the usual definition of Jeans' length (up to a factor $2 \pi$, e.g. Peebles 1980). Equation (34) can be solved analytically for particular choices of the cosmological background and the dependence  $p(\varrho)$ (Haubold et al. 1991). Qualitatively, one finds that the pressure term is irrelevant for inhomogeneities on scales $\gg$$\lambda$, which grow according to the "dust model'', while those on scales $\ll$$\lambda$ are damped by pressure and may exhibit sound-like oscillations.

The "EJN model'' provides a correction to the parallelism relation already at the linear order. Following Buchert et al. (1999), the solution to Eq. (34) and the departure from parallelism can be computed as a series in powers of $\lambda$ (unlike in the fully non-linear problem, the limit $\lambda \rightarrow 0$ in Eq. (34) is regular). Asymptotically in time, we may write ( ${\cal C}(t)$ is a dimensionless time-dependent coefficient):

 

$$% {\vec{w}}^{(1)} F (t) \left[ {\vec{u}}^{(1)} - \frac{{\cal C}(t)}{a^2} \lambda^2\... ...right) \Delta_{\vec{q}}{\vec{u}}^{(1)} + {\cal O}\left(\lambda^4\right) \right]$$ =

$$F(t) {\vec{u}}^{(1)} + {\cal C}(t) \frac{p'\left(\varrho_{\rm H}\... ...varrho_{\rm H}} \nabla_{\vec{q}}\varrho^{(1)} + {\cal O}\left(\lambda^4\right),$$ = (35)

where the last equality follows by using Eqs. (31c,18). The correction is proportional to the extra acceleration due to the pressure gradient. Since  ${\cal C}(t) \neq 1$in general, the peculiar-velocity is not parallel to the total peculiar-acceleration (gravity+pressure).

   
The non-linear orders, ${{\cal O}(\epsilon ^n)}$, n>1

The Eulerian perturbative expansion beyond the linear order has been intensively studied for the "dust model'' in recent years (see, e.g., Bernardeau et al. 2002 and references therein). The predictions of the "SSE model'' differ from those of the "dust model'' at non-linear orders. An important difference concerns the vorticity. With full generality (i.e., to all orders of the perturbative expansion), one can easily show from Eq. (13b) by application of the vector identity $({\vec{u}}\cdot \nabla_{\vec{q}}) {\vec{u}}= (1/2) \nabla_{\vec{q}}\vert{\vec{u}}\vert^2 - {\vec{u}}\times (\nabla_{\vec{q}}\times {\vec{u}})$ that (here, a comma denotes partial derivative with respect to comoving Eulerian coordinates)

 

$$% \frac{\partial \omega_i}{\partial t} + 2 H \omega_i + \frac{1}{... ...l_j \left[ \frac{1}{\varrho} \left(F_k - \frac{1}{a} \Pi_{km,m}\right) \right].$$ (36)

For the "dust model'' (vanishing r.h.s.), this equation can be integrated exactly along the flow lines, Eq. (39) (e.g., Serrin 1959; Buchert 1992 Eq. (26)ff):

$$% \omega_i ({\vec{X}},t) = \left(\frac{a_0}{a}\right)^2 \; \f... ...ial P_i ({\vec{X}},t)}{\partial X_j} \right) \right\vert}\cdot$$ (37)

We see that only the corrections to "dust'' can be a source of vorticity: if the initial vorticity vanishes (as is usually assumed in cosmological models of structure formation by virtue of the decay of the linear vorticity), then the evolution according to the "dust model'' predicts that the vorticity remains zero (Kelvin's circulation theorem). This remains true for the "EJN model'', because the pressure is represented as a function of density only and the r.h.s. of Eq. (36) vanishes. The "SSE model'', however, does incorporate the non-linear generation of vorticity: when Eq. (36) is evaluated to order ${\cal O}(\epsilon^2)$ with the expansions (29), one obtains:

$$% \frac{\partial \omega_i^{(2)}}{\partial t} + 2H \omega_i^{(... ...{1}{a} \left(u_{m,l}^{(1)} u_{k,l}^{(1)}\right)_{,jm} \right].$$ (38)

The source terms of vorticity are the tidal and shear stretching and can be evaluated as a functional solely of the linear density field $\varrho^{(1)}$ by resorting to Eqs. (31c), (32). A detailed study of the growth of the ensemble average $\langle \omega^2 \rangle$ and its dependence on the initial conditions can be found in (Domínguez 2002). As soon as the velocity field has a non-vanishing vortical component, parallelism with the gravitational peculiar-acceleration is rigorously ruled out. (But this does not mean that deviations from parallelism must be necessarily large or important.) Vorticity acts against collapse (see Sect. 6) and is likely to play a role in shaping the internal structure of collapsed regions, which in the "adhesion model'' are structureless.

   
5 The Lagrangian perturbative expansion

To perform the Lagrangian perturbative expansion, we have to transform the equations with respect to Lagrangian coordinates and follow the field quantities along comoving trajectories of the fluid elements:

$$% {\vec{q}}= {\vec{X}}+ {\vec{P}}({\vec{X}},t) , \qquad {\vec... ...{\dot{\vec{P}}} , \qquad {\vec{P}}({\vec{X}},t_{0})={\vec{0}},$$ (39)

with the Lagrangian time-derivative operator $\dot{ } := \frac{\partial}{\partial t}\vert_{\vec{X}}=$  $\frac{\partial}{\partial t}\vert_{\vec{q}}+$ $\frac{u_i}{a} \frac{\partial}{\partial q_{i}}$. (Note that ${\vec{X}}+ {\vec{P}}$ are integral curves of the mean peculiar-velocity, scaled by a(t)). ${\vec{X}}$ are Lagrangian coordinates and  ${\vec{P}}({\vec{X}},t)$ is the displacement of the fluid volume elements away from their initial positions; this is the expansion parameter in the Lagrangian perturbative approach. This method has been applied successfully to study the "dust model'', the "Zel'dovich approximation'' being a subcase of the first-order solution, as well as higher-order solutions. It has been recently applied to the "EJN model'' too (see the references below). Concerning the "SSE model'', there is no difference with the "dust model'' to linear order.

In view of later considerations, for the application of the Lagrangian methods we take an unconventional route and first rewrite Eqs. (13) as follows. We obtain an evolution equation for the peculiar-gravitational field ${\vec{w}}$ by eliminating the density in Eq. (13a) with the help of Eq. (13c). By formally integrating the divergence one gets (Buchert 1989):

 

$$\dot{{\vec{w}}} + 2H{\vec{w}}- 4\pi G \varrho_{\rm H} {\vec{u}}= \vec{\cal R},$$ (40a)

$$\vec{\cal R}: = \frac{1}{a} \left[\left({\vec{u}}\cdot {\nabla}_{... ...a}_{\vec{q}}\cdot {\vec{w}}\right) + {\nabla}_{\vec{q}}\times {\vec{T}}\right].$$

The vector ${\vec{T}}$ is essentially fixed by the other equation satisfied by ${\vec{w}}$, Eq. (13d): acting with $\nabla_{\vec{q}}\times$ on Eq. (40a) and subjecting ${\vec{T}}$ to the Coulomb gauge condition, ${\nabla}_{\vec{q}}\cdot {\vec{T}}= 0$, one finds

$$% \Delta_{\vec{q}}{\vec{T}} 4\pi Ga {\nabla}_{\vec{q}}\times (\varrho {\vec{u}}) = 4\pi Ga \l... ...a}_{\vec{q}}\times {\vec{u}}+ {\nabla}_{\vec{q}}\varrho \times {\vec{u}}\right]$$ =

$$\left(4\pi Ga \varrho_{\rm H} - {\nabla}_{\vec{q}}\cdot {\vec{w}}... ...{\nabla}_{\vec{q}}\times {\vec{u}}+ {\vec{u}}\times{\Delta}_{\vec{q}}{\vec{w}}.$$ = (40b)

That is, ${\vec{T}}$ is the vector potential (up to an unimportant factor) of the peculiar-current density $\vec{j}: =\varrho {\vec{u}}$ (i.e., $4 \pi G a \vec{j}=-\nabla_{\vec{q}}\times {\vec{T}}-\nabla_{\vec{q}}~\frac{1}{a}\partial_{\rm t} (a\phi)$, with $\phi$denoting the scalar peculiar-gravitational potential). (If ${\vec{w}}= :-\frac{1}{a}\nabla_{\vec{q}}\phi$ is initially irrotational, Eqs. (40a,b) propagate this condition in time.) Finally, Eq. (13b) gives the evolution of ${\vec{u}}$:

$$% \dot{{\vec{u}}} + H {\vec{u}}= {\vec{w}}+ \frac{1}{\varrho}... ...vec{F}}- \frac{1}{a} \nabla_{\vec{q}}\cdot {\vec{\Pi}}\right).$$ (40c)

Thus, we have replaced Eqs. (13) for $\varrho$, ${\vec{u}}$, and ${\vec{w}}$ by the new set of Eqs. (40) for ${\vec{u}}$ and ${\vec{w}}$. They hold independently of what the correction to "dust'' looks like. These equations still retain the nonlocal nature imposed by gravity through the vector potential ${\vec{T}}$ given by Eq. (40b); we obtain "local'' evolution equations only if the source of this equation vanishes. Physically this means that the peculiar-current density has to be irrotational, and, for irrotational flows, that the mean flow follows the gradient of the density field. In the following, we will not be further concerned with the residual vector field  $\vec{\cal R}$. We just acknowledge that it is a functional of the dynamical fields ${\vec{u}}$ and ${\vec{w}}$ so that the system of Eqs. (40a-c) is closed for given ${\vec{F}}$ and  ${\vec{\Pi}}$.

  
5.1 The single-streamed case of "dust matter''

For the "dust model'' the system of Eqs. (40a-c) has been studied thoroughly in (Buchert 1989). There, the currently used notion of the weakly non-linear regime has been defined as "Lagrangian linearization''. The comoving trajectory field ${\vec{P}}$ of Zel'dovich's approximation (Zel'dovich 1970, 1973) can be understood as a subclass of the linear solutions in Lagrangian space (Buchert 1989, 1992). For an Einstein-de Sitter cosmology with initial peculiar-velocity ${\vec{u}}({\vec{X}}, t_{0})$, the subclass corresponding to Zel'dovich's model is ${\vec{P}}({\vec{X}}, t) = \frac{3}{2} [b(t)-b(t_{0})] {\vec{u}}({\vec{X}}, t_{0}) ~ t_{0}$, where b(t)=a(t) here. For all background cosmologies including "curvature'' and a cosmological constant Zel'dovich's approximation is given in (Bildhauer et al. 1992; for the flat models with cosmological constant, see also the supplement by Chernin et al. 2003).

This well-known result can be easily obtained from the general system of Eqs. (40a-c) by the following reasoning: in the Lagrangian picture the Eulerian dynamical variables  ${\vec{u}}({\vec{q}},t)$ and  ${\vec{w}}({\vec{q}},t)$ are replaced by the single dynamical variable  ${\vec{P}}({\vec{X}},t)$. Lagrangian linearization means that the equations have to be linearized with respect to ${\vec{P}}$. Equations (39) and (40c) (without corrections to "dust'') express the fields ${\vec{u}}$ and ${\vec{w}}$ as linear functions of ${\vec{P}}$. It follows that Eqs. (40a,b) have to be linearized in the fields ${\vec{u}}$ and ${\vec{w}}$. For simplicity we restrict the initial conditions to irrotational flow, so that the velocity remains irrotational through the "dust'' evolution, Eq. (37). Then, Eq. (40b) implies, upon linearization, ${\vec{T}}={\vec{0}}$. From Eq. (40a) we can therefore immediately drop the manifestly non-linear terms in the residual vector field  $\vec{\cal R}$. Introducing the longitudinal and transversal parts of ${\vec{u}}$ and ${\vec{w}}$ with respect to Lagrangian coordinates, e.g., ${\nabla}_{\bf0}\cdot{\vec{w}}^T = 0$, ${\nabla}_{\bf0}\times{\vec{w}}^L = {\bf0}$, where  ${\nabla}_{\bf0}$ denotes the nabla operator with respect to Lagrangian coordinates, the remaining equations for "dust'' read:

 

$${\dot{\vec{w}}^L} + 2H {{\vec{w}}^L} - 4\pi G \varrho_{\rm H} {\vec{u}}^L = {\bf0},\;\; {\dot{\vec{u}}^L} + H {{\vec{u}}^L} - {\vec{w}}^L = {\bf0},$$

$${\vec{u}}^T = {\vec{w}}^T = {\vec{0}}.$$ (41)

These equations are satisfied if the inhomogeneous displacement field ${\vec{P}}$, split into its longitudinal and transversal parts, obeys:

 

$${\ddot{\vec{P}}^L} + 2 H {\dot{\vec{P}}^L} - 4\pi G \varrho_{\rm H} {\vec{P}}^L = {\bf0},$$ (42)

$${\vec{P}}^T = {\bf0},$$

which can be easily demonstrated by reinserting the definitions of the fields ${\vec{u}}$ and ${\vec{w}}$ in favor of ${\vec{P}}= {\vec{P}}^L + {\vec{P}}^T$. If we restrict the general solution of Eqs. (41) to the initial data ${\vec{w}}({\vec{X}},t_{0}) = ({\rm const.}/t_{0}) ~ {\vec{u}}({\vec{X}},t_{0}) \;\;({\rm const}. = 1$ for an Einstein-de Sitter cosmology), we obtain Zel'dovich's approximation. The rigorous calculation, starting from the fully transformed system of equations, may be found in (Buchert 1989, 1992).

5.2 The multi-streamed case modelled by dynamical stresses

In the regime when the "dust model'' predicts multi-streaming, the corrections to "dust'' become relevant. We look for the Lagrangian linear form by also dropping the residual vector field  $\vec{\cal R}$ (since initial irrotationality is preserved to linear order by the evolution in both the "EJN'' and the "SSE model''). The "SSE'' correction to "dust'' is also non-linear, Eq. (23b), so that this model does not differ from the "dust model'' to linear order. For the "EJN model'', Eq. (17b), however, we have to linearize additionally the (Eulerian) Laplacian together with the GM coefficient along the comoving trajectory field. Using transformation tools developed by Adler & Buchert (1999, Appendix A) we find (for isotropic  ${\vec{\Pi}}$):

$$% \frac{\lambda^2 (\varrho)}{a^2} {\Delta}_{\vec{q}}{\vec{w}}... ...\Delta}_{\vec{0}}{\vec{w}}+ {\cal O}\left({\vec{P}}^2\right).$$ (43)

The complete set of Lagrangian linear equations may then be cast into a form that reproduces the "dust'' result (41) on their left-hand-sides:

 

$${\dot{\vec{w}}^L} + 2H {{\vec{w}}^L} - 4\pi G \varrho_{\rm H} {\vec{u}}^L = {\bf0},$$

$${\dot{\vec{u}}^L} + H {{\vec{u}}^L} - {\vec{w}}^L = \frac{1}{a^2}\lambda^2 \left(\varrho_{\rm H}\right) {\Delta}_{\bf0} {\vec{w}}^L,$$ (44)

$${\vec{u}}^T = {\vec{w}}^T = {\vec{0}}.$$

Again, the rigorous way of deriving the equation obeyed by Lagrangian linear displacements is the following: we first transform the basic equations into Lagrangian form, then linearize the Lagrangian system with respect to ${\vec{P}}$. This has been done in (Adler & Buchert 1999) with the result that ${\vec{P}}$ obeys:

 

$${\ddot{\vec{P}}^L} + 2 H {\dot{\vec{P}}^L} - 4\pi G \varrho_{\rm ... ... \frac{\lambda^2 \left(\varrho_{\rm H}\right)}{a^2} {\Delta}_{\bf0}{\vec{P}}^L,$$

$${\ddot{\vec{P}}^T} + 2 H {\dot{\vec{P}}^T} ={\bf0};\;\;\;{\mbox{here:}}\;\;\;{\vec{P}}^T ={\bf0}.$$ (45)

After some (in this case more involved) manipulations one can find the equations for ${\vec{u}}$ and ${\vec{w}}$ above by inserting their definitions in terms of ${\vec{P}}$ into Eqs. (45). The differential operator acting on  ${\vec{P}}^L$ is formally identical to Eqs. (34), and Eqs. (45) can thus be solved by extrapolation of the known Eulerian linear solutions (Adler & Buchert 1999).

Evidently, the adhesive term consists of a Lagrangian Laplacian and, therefore, it involves, in addition to the convective non-linearities hidden in the overdot, non-linearities in Eulerian space when mapping it back using the inverse solution of the mapping ${\vec{q}}= {\vec{X}}+ {\vec{P}}({\vec{X}},t)$.

In general, parallelism of peculiar-velocity and -acceleration is violated. By computing the time derivative of the difference  ${\vec{w}}- F(t) {\vec{u}}$ with Eqs. (45), one can show that in general it does not decay to zero asymptotically for large times. The correction is proportional to  ${\Delta}_{{\vec{0}}} {\vec{u}}$ to lowest order in $\lambda^2$. This correction is similar to the one found in the Eulerian linear approximation, Eq. (35), but in terms of the Lagrangian gradients.

5.3 Beyond first order

The Lagrangian perturbation scheme including pressure (Adler & Buchert 1999) has been pushed to second order (Morita & Tatekawa 2001; Tatekawa e al. 2002) and recently to third order (Tatekawa 2005). Results on cross-correlation statistics with N-body simulations (Tatekawa 2004a) and a comparison of the corresponding density fluctuations with results of hydrodynamical simulations (Tatekawa 2004b) has demonstrated that multi-stream forces tend to be underestimated by Lagrangian perturbative models. In the next section we shall investigate why this happens.

   
6 The non-perturbative regime

   
6.1 Exact evolution equations I: General case

Following a systematic perturbative approach we could provide models for adhesive gravitational clustering in the weakly non-linear regime. This has allowed us also to relax the parallelism assumption (24) that was a necessary ingredient of the standard "adhesion model'' but, as we have shown, is not expected to hold in the multi-streamed regime.

Finding an extension into the non-perturbative (both Eulerian and Lagrangian) regimes is an involved mathematical task. Notwithstanding, it should be attempted in view of the fact that even the Lagrangian perturbation approach falls short of capturing the action of multi-stream forces, since a genuine property of adhesive models in the simplest cases involves an Eulerian Laplacian. Thus, we expect the best approximations to be of hybrid Lagrangian/Eulerian type like the standard "adhesion model'', meaning a non-linear equation both in Eulerian space (due to the convective non-linearities) as well as in Lagrangian space (due to the Eulerian gradients in the forcing).

An exact equation for the comoving deviations of the trajectory field, ${\vec{P}}({\vec{X}},t)$ defined by Eq. (39), can be obtained from Eqs. (40): formal integration of Eq. (40a) yields

$$% {\vec{w}}({\vec{X}}, t) = 4 \pi G a \varrho_{\rm H} {\vec{P... ... \!\!\! {\rm d}\tau \; a^2(\tau) \vec{\cal R}({\vec{X}},\tau),$$ (46)

with the integration constant chosen such that the Lagrangian coordinates are defined by the initial datum for ${\vec{w}}$, $4 \pi G a_{0}\varrho_{\rm H} (t_{0}) {\vec{P}}({\vec{X}}, t_{0}) = {\vec{w}}({\vec{X}}, t_{0})$. Then, from Eq. (40c) one obtains the general result:

 

$$% \ddot{{\vec{P}}} + 2 H \dot{{\vec{P}}} - 4 \pi G \varrho_{\rm H... ...varrho} \left[ {\vec{F}}- \frac{1}{a} \nabla_{\vec{q}}\cdot {\vec{\Pi}}\right].$$ (47)

The left-hand-side of this equation reproduces the "dust model''. If the correction to "dust'' is expressed in terms of the fields $\varrho$, ${\vec{u}}$, ${\vec{w}}$, then the equation can be written entirely in terms of ${\vec{P}}$ and  $\vec{\cal R}$ by virtue of Eqs. (39), (46), (51). Thus, the difficulty in deriving a closed, local differential equation for the deviation field ${\vec{P}}$ lies in the term  $\vec{\cal R}$, which is a nonlocal (in space and time) and non-linear functional of ${\vec{P}}$ given by Eq. (40b).

As a further illustration of the general case we write down the corresponding equation for the density. For this purpose we introduce the contrast density field $\Delta = (\varrho_{\rm H} - \varrho)/\varrho$, $-1 < \Delta < \infty$, which is more adapted to the non-linear situation than the conventional density contrast $\delta : = (\varrho - \varrho_{\rm H} )/\varrho_{\rm H}$ (for the "dust'' case see: Buchert 1989 Eq. (7)ff, 1992 Eq. (31)ff, and 1996 Eq. (25)ff (with a different sign convention for $\Delta$)). The exact evolution equation for this variable can be obtained from Eqs. (13a,b):

 

$$\ddot \Delta + 2 H \dot \Delta - 4\pi G \varrho_{\rm H} \Delta = ... ...vec{F}}- \frac{1}{\varrho a}{\nabla}_{{\vec{q}}}\cdot{\vec{\Pi}}\right)\right],$$ (48)

with the second scalar invariant of the peculiar-velocity gradient,

$$II^{\rm pec}:= \frac{1}{2 a^2}\left[(u_{i,i})^2 - u_{i,j}u_{j... ...ec{u}}- {\vec{u}}\cdot {\nabla}_{{\vec{q}}} {\vec{u}}\right).$$

For the purpose of evaluating the Eulerian spatial derivatives in the right-hand-side of the above general equations, we employ the following transformation:

$$% \frac{\partial}{\partial q_i} = \left(J^{-1}\right)_{ji} ~ \frac{\partial}{\partial X_j},$$ (49)

where J^-1 is the inverse to the following matrix:

$$% J_{ij}({\vec{X}},t) := \delta_{ij} + \frac{\partial P_i({\vec{X}},t)}{\partial X_j}\cdot$$ (50)

In the rest of this section we assume always that the map ${\vec{X}}\to {\vec{q}}$is invertible, i.e., no shell-crossing, so that $\det J_{ij} >0$. Mass conservation, Eq. (13a), relates the density field to ${\vec{P}}$ via the matrix J_ij:

 

$$% \varrho({\vec{X}}, t) \left[\frac{a_{0}}{a}\right]^3 \varrho ({\vec{X}}, t_{0}) ~ \frac{\det J_{ij} ({\vec{X}},t_{0})}{\det J_{ij} ({\vec{X}},t)}$$ =

$$\approx \varrho_{\rm H}(t) / \det J_{ij} ({\vec{X}},t),$$ (51)

where we take initial quasi-homogeneity for simplicity ( $\varrho ({\vec{X}},t_{0}) \approx \varrho_{\rm H} (t_{0})$, $J_{ij}({\vec{X}},t_{0}) \approx \delta_{ij}$).

   
6.2 Exact evolution equations II: Plane symmetry

The problem is simplified if one considers particular geometric settings with high symmetry. Although this leads to not fully realistic models, it can provide hints of the mathematical and physical properties of the models in the non-perturbative regime. Plane-symmetric models provide an excellent test case for numerical simulations of multi-stream systems concerning a variety of properties including spatial scaling (e.g. Doroshkevich et al. 1980; Gouda & Nakamura 1988, 1989; Yano & Gouda 1998; Fanelli & Aurell 2002; Aurell et al. 2003; Yano et al. 2003; the assumption of spherical symmetry plays a complementary role, e.g. Padmanabhan 1996.) In such configurations, the fields vary spatially along a single direction, and so render the Lagrangian approach particularly suitable for an extension into the non-perturbative regime. Plane-symmetric models simplify the problem enormously, since the residual term  $\vec{\cal R}$ vanishes, and the "dust'' evolution remains in the Lagrangian linear regime at all times (implying in particular that the parallelism relation (24) - stating a spatially constant factor of proportionality - holds asymptotically at large times).

In a plane-symmetric configuration, the fields vary only along a single direction, say i=1. Equations (46), (47), (49), (51) then are simplified to (we hereafter drop the index i=1):

 

$$\varrho(X, t) = \varrho_{\rm H}(t) ~ \left[ 1 + \frac{\partial P(X,t)}{\partial X} \right]^{-1},$$

$$\frac{\partial}{\partial q} = \left[ 1 + \frac{\partial P(X,t)}{\partial X} \right]^{-1} \frac{\partial}{\partial X},$$

$$w (X, t) = 4 \pi G a \varrho_{\rm H} P(X,t),$$ (52)

$$\ddot{P} + 2 H \dot{P} - 4 \pi G \varrho_{\rm H} P = \frac{1}{a \varrho} F - \frac{1}{a^2 \varrho_{\rm H}} \frac{\partial \Pi}{\partial X}\cdot$$

Therefore, depending on the specific model, we get different equations:

• "Dust model'':
  

  $$% \ddot{P} + 2 H \dot{P} - 4 \pi G \varrho_{\rm H} P = 0.$$ (53a)

  
  

• "EJN model'': inserting Eqs. (16) particularized to the plane-symmetric configuration, one obtains:
  

  $$% \ddot{P} + 2 H \dot{P} - 4 \pi G \varrho_{\rm H} P = - \f... ... \frac{\varrho_{\rm H}}{1+\partial P/\partial X} \right)\cdot$$ (53b)

  
  

• "SSE model'': using Eqs. (19),
  

  $$% \ddot{P} + 2 H \dot{P} - 4 \pi G \varrho_{\rm H} P = - B L^2 \d... ...t(\partial \dot{P}/\partial X\right)^2}{\left(1+\partial P/\partial X\right)^3}$$

  $$+ 4 \pi G \varrho_{\rm H} \left[ \displaystyle\frac{1}{2 \left(1+... ...right)^2} - \displaystyle\frac{1}{1+\partial P/\partial X} \right] \Bigg\}\cdot$$ (53c)

  
  

Equations (53) are closed partial differential equations for the function P(X,t). Notice that the equation for the "dust model'' is linear, so that Lagrangian linearization already yields the exact solution. The equation for P(X,t) in the "EJN model'' becomes also a linear equation for the special dependence $p \propto \varrho^{-1}$ (Chaplygin gas, App. C; cf. Adler & Buchert 1999). This case is formally of interest because it provides a subcase in which the Lagrangian linear solution exactly solves the plane-symmetric problem.

Equations (53) should allow one to explore the role of the corrections to "dust'' beyond the assumptions employed in Sect. 3.1, in particular, the back-reaction on the "dust'' trajectories by the correction terms. It should be noted that Eqs. (53) are of second order in time, in contrast to the generalized "adhesion model'', Eq. (26a) - the relevance of this fact concerning the formation of singularities is discussed in Sect. 6.4. In this connection Eqs. (53b,c) have features of a quasilinear hyperbolic ("wave'') equation, while Eq. (26a) is a quasilinear parabolic ("diffusion'') equation. Götz (1988) considered the "EJN model'' with the "isothermal'' relationship $p \propto \varrho$ in the absence of background cosmological expansion ( $a(t) \equiv 1$), and he was able to show that it can be mapped exactly to the Sine-Gordon equation, which is known to admit soliton solutions.

Fanelli & Aurell (2002) numerically studied the one-dimensional Vlasov-Poisson system (see Sect. 2), in particular the time-dependence of a collapsing isolated perturbation. Fanelli & Aurell tried to interpret the numerical findings in the framework of the simplified "EJN model'' (26) with a polytropic relationship for the pressure, $p(\varrho) \propto \varrho^\gamma$. A simple scaling argument led them to the value $\gamma=7/9$ (the original "adhesion model'' corresponding to $\gamma=2$). This conclusion is, however, questionable because boundary-layer theory cannot be applied to Eqs. (26) when $\gamma<1$. A more general analysis (Domínguez, unpublished) shows that a shock can form in this case only when the initial velocity is smaller than the maximum velocity, and this maximum velocity vanishes in the limit $\mu \rightarrow 0$. That is, when $\gamma<1$, the adhesive term is too weak to prevent singularities for arbitrary initial velocities. Thus, it might be the case that Eqs. (26) cannot describe even qualitatively the numerical experiment by Fanelli & Aurell.

6.3 A new non-perturbative approximation

As remarked, one of the difficulties with the exact Eq. (47) is the nonlocal nature introduced by  $\vec{\cal R}$. A possible approximation consists therefore in setting $\vec{\cal R}={\vec{0}}$, so that Eq. (47) simplifies to

$$% \ddot{{\vec{P}}} + 2 H \dot{{\vec{P}}} - 4 \pi G \varrho_{\... ...vec{F}}- \frac{1}{a} \nabla_{\vec{q}}\cdot {\vec{\Pi}}\right].$$ (54)

As long as ${\vec{F}}$ and ${\vec{\Pi}}$ are modelled as local functions of $\varrho$, ${\vec{u}}$ and ${\vec{w}}$ (as in the models we consider in this work), this is a partial differential equation for ${\vec{P}}$ by virtue of Eqs. (39), (46), (51), albeit a difficult one.

One can view the approximation $\vec{\cal R}={\vec{0}}$ in the same spirit as the parallelism assumption (24): one gets rid of the problem of nonlocality by approximating the exact peculiar-gravitational field by a relationship which holds in the (both Eulerian and Lagrangian) linear regime of the "dust model'' as well as exactly in the plane-symmetric configuration:

$$% {\vec{w}}= 4 \pi G a \varrho_{\rm H} {\vec{P}},$$ (55)

to be compared with ${\vec{w}}= 4 \pi G a \varrho_{\rm H} (b/\dot{b}) \dot{{\vec{P}}}$corresponding to Eq. (24). The relation (55) replaces the "slaving'' assumption of parallelism of ${\vec{w}}$ to the peculiar-velocity. It is more general and simpler, because it does not involve specification of the function b(t). Unlike the condition of parallelism, this approximation also retains features of the exact field equations for ${\vec{w}}$: e.g., at any given instant of time, the approximated field ${\vec{w}}$ depends only on the displacement (i.e., the distribution of mass) at this time, and not on its time-derivative (the peculiar-velocity field). This renders the differential equation obeyed by ${\vec{P}}$second-order in time (see Sect. 6.4 in this respect). A proper understanding of the extrapolation of the linear relationship between ${\vec{w}}$ and ${\vec{P}}$ into the non-perturbative regime requires that the residual term  $\vec{\cal R}$ be analyzed, a task beyond the goals of this work.

To get a glimpse of the mathematical difficulties posed by Eq. (54), we particularize it to the "EJN model'' of ${\vec{F}}$and  ${\vec{\Pi}}$, Eq. (16). The right-hand-side can be written, in analogy with Eq. (17b), as:

$$% \frac{1}{a \varrho} \left[ {\vec{F}}- \frac{1}{a} \nabla_{... ...rm H} ~ p'(\varrho)}{a^2 \varrho} {\Delta}_{\vec{q}}{\vec{P}}.$$ (56)

(In the expressions for ${\vec{F}}$ and ${\vec{\Pi}}$, the fields $\varrho$ and ${\vec{w}}$ are interchangeable via Eqs. (13c,d), but the resulting different equations are no longer exactly equivalent when the approximation $\vec{\cal R}={\vec{0}}$ is used.) This term can be simplified with the "virial equation of state'' $p = \kappa \varrho^2$, so that the prefactor of the Laplacian is density-independent. This expression for  $p(\varrho)$ is certainly a good phenomenological model consistent with the results of numerical simulations (Domínguez 2003; Domínguez & Melott 2004). The results presented in Sect. 3.1 also suggest that the precise value of the polytropic exponent does not alter the qualitative behavior in high-density regions. With this simplification, Eq. (54) becomes:

$$% \ddot{{\vec{P}}} + 2 H \dot{{\vec{P}}} - 4 \pi G \varrho_{\... ...ac{2 \kappa \varrho_{\rm H}}{a^2} {\Delta}_{\vec{q}}{\vec{P}}.$$ (57)

This equation, being non-linear in Eulerian and in Lagrangian space, belongs to the class of hybrid Lagrangian/Eulerian models mentioned in Sect. 6.1. As such, the approximation may be viewed - in contrast to the Lagrangian perturbation scheme - in closer correspondence to the former "adhesion approximation''. By construction, this approximation reduces to the full Lagrangian linear part of the "dust model'', Eq. (42), when the correction to "dust'' is dropped. Upon linearizing the Eulerian Laplacian, we recover the Lagrangian linear approximation, Eq. (45). (In other words, our non-perturbative approximation extrapolates the Lagrangian linear model by the replacement $\Delta_{{\vec{X}}} \rightarrow \Delta_{{\vec{X}}+ {\vec{P}}}$.) Linearization of the convective non-linearities yields the Eulerian linear approximation, Eq. (34) (for this, note that linearization of Eq. (51) gives $\delta = - \nabla_{\vec{0}}\cdot {\vec{P}}$).

In conclusion, we propose the approximation (55) as a first extrapolation into the non-perturbative regime of the Eulerian and Lagrangian linear approximations without resort to the constraining assumption of parallelism. It should be viewed as a simplification of the exact Eqs. (13) in the sense that the approximation yields a local equation, which should facilitate the theoretical analysis. An attempt to solve the approximate equation could be successful, but lies beyond the scope of the present work.

   
6.4 Regularization of singularities

A point one would like to be able to prove rigorously beyond the argument offered in Sect. 3.1 is whether the "EJN and SSE models'' do indeed prevent the formation of singularities given initial conditions of cosmological relevance. The absence of singularities when initially smooth data are propagated by the Vlasov-Poisson system has been proved mathematically (Lions & Perthame 1991; Schaeffer 1991; Pfaffelmoser 1992; Rein & Rendall 1994). The original "adhesion model'' reduces to the 3D Burgers' equation (Eq. (26) with constant $\mu$), which has the unusual property of being solvable analytically. From there, one can demonstrate that no singularity arises for any sufficiently smooth initial conditions (see, e.g., Frisch & Bec 2001). For other models this is a very difficult question (e.g., it is still an open problem for the 3-dimensional incompressible Navier-Stokes equations (Frisch 1995, Sect. 9.3)), and here we can only provide some general remarks.

Starting from Eqs. (13), (21), one can derive the following evolution equations for the density $\varrho$ and the peculiar-expansion rate $\theta$:

$$% \dot{\varrho} = - 3 H \varrho - \frac{1}{a} \varrho ~ \theta,$$ (58a)

$$% \dot{\theta} = - 2 H \theta - \frac{1}{3} \theta^2 - \sigma_{ij... ...arrho} {\vec{F}}- \frac{1}{a \varrho} \nabla_{\vec{q}}\cdot {\vec{\Pi}}\right).$$ (58b)

The continuity Eq. (58a) can be integrated in Lagrangian coordinates:

$$% \varrho({\vec{X}}, t) = \varrho({\vec{X}}, t_{0}) ~ \left[... ...}\tau \; \frac{\theta ({\vec{X}}, \tau)}{a(\tau)} \right]\cdot$$ (59)

Starting from bounded initial conditions, a singularity in the density field can occur only when $\theta \rightarrow -\infty$ at some finite time. From Eq. (58b) we see that only three terms can actually prevent $\theta$ from becoming too negative: the background expansion, the vorticity, and perhaps the model-dependent correction to "dust''. In particular, for the "dust model'' with vanishing initial vorticity, only the term linear in $\theta$in Eq. (58b) can oppose collapse, so that a density singularity is unavoidable provided there is an initial overdensity.

For the "EJN model'', the correction to "dust'' reads $-\nabla_{\vec{q}}\cdot (\varrho^{-1} \nabla_{\vec{q}}p)$, where p is modelled as a density-dependent pressure, so that Eqs. (17) have the structure of the equations of inviscid fluid dynamics. In the absence of self-gravity, it is known that the solution becomes non-differentiable for a generic class of smooth initial conditions, so that the Lagrangian-to-Eulerian map  ${\vec{q}}({\vec{X}},t)$ ceases to be defined. The map can remain uni-valued if we allow for shocks (moving discontinuities in the fields): this is called a weak solution, and they have to be understood as solutions of the differential Eqs. (17) with an additional viscous term  $\mu {\Delta}_{\vec{q}}{\vec{u}}$ in the limit  $\mu \to 0^+$, in the same way as we discussed for the adhesion model. The presence of self-gravity is not expected to alter this conclusion; the question is rather whether the initial conditions of cosmological interest belong to the class of initial conditions inducing shock formation. The reason for the emergence of shocks can be ultimately traced back to the fact that density and peculiar-velocity are independent variables: when convection by the peculiar-velocity field occurs at such a high rate that the density gradients are not enough to build up a counteracting pressure, a singularity arises. However, in the cosmological context the peculiar-velocity field becomes "slaved'' to the density field in the linear regime (see Sect. 4.1). The assumption of parallelism (24) extrapolates this "slaving'' relationship into the non-linear regime: this allows one to express the density-dependent correction to "dust'' in terms of the velocity gradient, and it seems sufficient to avoid singularities (see Sect. 3.1). But we have seen that a departure from the condition of parallelism is generated during the evolution: an open question is therefore whether the assumption of parallelism is qualitatively good or instead the dynamically generated departure from it can lead to the formation of singularities.

The mathematical properties of the equations underlying the "SSE model'' have been barely studied in the literature. The correction to "dust'' is a complicated function of the fields. Unlike the "EJN model'', the corrections depend explicitly on the gradients of the peculiar-velocity field too, and the analogy with shock formation in inviscid fluids is not valid. Another important difference to the "dust'' and the "EJN'' models is that vorticity is generated by the "SSE'' corrections to "dust'' (see Sect. 4.2), which, as pointed out, tends to oppose collapse. Vreman et al. (1996, 1997) have studied numerically the hydrodynamical equations of an ideal gas with an additional stress tensor given by Eq. (19b) when the initial condition is a perturbation to a smooth profile. An instability was discovered for which the "SSE''-term is made responsible. Vreman et al. (1996) also undertook a theoretical analysis of the following simplified model (one-dimensional forced Burgers equation with a simplified, "SSE''-like correction; $\eta>0$):

$$% \frac{\partial u}{\partial t} + u \frac{\partial u}{\partia... ...al x} \left( \frac{\partial u}{\partial x} \right)^2 +f(x).$$ (60)

It is found that a sinusoidal profile (which fixes the form of the forcing f(x) by demanding stationarity) is linearly unstable to perturbations if $2 \eta > \mu-1$. Indeed, as the discussion of Eq. (62) below exemplifies, it is not true that the "SSE''-corrections will behave like viscous terms under all conditions. As with the "EJN model'', the question is whether the "SSE model'' does indeed prevent the formation of singularities when the initial conditions are restricted to be of cosmological relevance.

An issue closely related to this problem is the behavior of the correction to "dust'' concerning the local "energy budget''. The density of mean peculiar-kinetic energy is defined as $K:=(1/2) \varrho \vert{\vec{u}}\vert^2$. From Eqs. (13) one finds:

 

$$% \frac{\partial K}{\partial t} + \frac{1}{a} \nabla_{{\vec{q}}} ... ...{\vec{u}}) - \frac{1}{a} \nabla_{{\vec{q}}} \cdot ({\vec{\Pi}}\cdot {\vec{u}}),$$ (61)

where : indicates contraction of two indices. The first term on the right-hand-side represents the work done by the total (mean field and corrections thereof) gravitational force, the second term is the work done by the velocity dispersion, and the third term is the contribution due to the flux of velocity dispersion. In the "EJN model'', Eqs. (16), only velocity dispersion does work: ${\vec{\Pi}}: (\nabla_{{\vec{q}}} {\vec{u}})=p(\varrho) \theta$. Given that p>0, this term tries to reduce the kinetic energy in (relative to the Hubble expansion) locally collapsing regions ( $\theta < 0$) and to increase it in locally expanding ones ( $\theta > 0$).

In the "SSE model'', the sign of the work done by the deviations from mean field gravity, ${\vec{u}}\cdot {\vec{F}}= B L^2 ({\vec{u}}\nabla_{{\vec{q}}} \varrho) : (\nabla_{{\vec{q}}} {\vec{w}})$, depends on the orientation of the tidal tensor and, like the work done by mean field gravity, is in general unrelated to whether the fluid is locally compressing or expanding. On the other hand, the work done by velocity dispersion is most easily expressed in the basis that diagonalizes the symmetric part of $(1/a) \partial_i u_j$ (this is $\sigma_{ij} + (1/3) \theta ~ \delta_{ij}$, according to the decomposition (21)):

 

$$% \frac{1}{a} {\vec{\Pi}}: (\nabla_{{\vec{q}}} {\vec{u}}) \frac{1}{a} B L^2 \varrho (\partial_k u_i) (\partial_k u_j) (\partial_i u_j)$$ =

$$2 B (a L)^2 \varrho \sum_i s_{(i)} \left[ s_{(i)}^2 + \frac{1}{4} \left(\omega^2-\omega_i^2\right) \right],$$ = (62)

where s_(i) are the eigenvalues of the tensor $\sigma_{ij} + (1/3) \theta ~ \delta_{ij}$, satisfying $\sum_i s_{(i)}=\theta$. The work is the sum of contributions along the three eigendirections: the contribution of the direction of local compression ( s_(i) < 0) tends to reduce the kinetic energy, that of the direction of local expansion ( s_(i) > 0) tends to increase it. The net sign is the outcome of this competition. In this sense, the "SSE'' expression for the velocity dispersion can be viewed as a "directional'' generalization of the case of an isotropic pressure like in the "EJN model''. It must be noticed, however, that this competing effect between local expansion and compression is effective only in the fully three-dimensional case: it can be easily checked that in the restricted cases of one or two-dimensional motion, expression (62) has the same sign as $\theta$.

   
7 Summary and prospects

Starting from the equations of motion for point particles, we have implemented a smoothing procedure in phase space leading to equations for the smoothed density and peculiar-velocity fields in Eulerian space. These equations are the first ones in an infinite hierarchy of equations, and approximations to close the hierarchy are required. A closure approximation corresponds physically to an ansatz on how the evolution of the filtered fields is coupled to the degrees of freedom below the filtering scales. In this unified context, we have addressed three different closures, leading to three different models proposed earlier in the literature: the "dust model'', the "EJN model'', and the "SSE model''. The "dust model'' describes the evolution of fluid volume elements under the sole influence of the self-consistently generated gravitational field. The corrections to "dust'' take into account the dynamical influence of the small-scale degrees of freedom. The "EJN model'' approximates this influence as a density-dependent pressure. The "SSE model'' expresses this influence as a functional of density and velocity in a systematic expansion.

We have provided a clear explanation of how the corrections to the "dust model'' by both the "EJN closure'' and the "SSE closure'' lead to models that are qualitatively equivalent to the "adhesion model''. The key hypothesis is the parallelism relation (24) between peculiar-velocity and peculiar-acceleration. This already holds in Zel'dovich's approximation to the exact "dust'' evolution - the assumption underlying the "adhesion approximation'' consists of extrapolating it to also evaluate the corrections to "dust''.

We have shown that within both the Eulerian and Lagrangian perturbation frameworks, the corrections to "dust'' violate this assumption in general. Going beyond the phenomenology of the "adhesion model'' therefore requires relaxing of the parallelism assumption. When this is done, differences between the "EJN model'' and the "SSE model'' appear: deviations from the parallelism relationship arise at linear (both Eulerian and Lagrangian) perturbative order in the "EJN model'', but not in the "SSE model''; vorticity is generated by the "SSE''-corrections to "dust'', but not by the "EJN''-corrections.

The non-perturbative theoretical investigation raises considerable difficulties. We suggested the plane-symmetric case as the "exact body'' of the non-perturbative equations that may serve as a starting point for a deeper understanding of the presented models. This case can be reduced to the study of a single partial differential equation for the displacement field, offering the possibility to study the "EJN'' and "SSE'' models beyond any approximation. This simplification misses some features that are expected to be of relevance for the detailed inner structure of collapsing high-density regions (e.g., generation of vorticity or the fact that velocity and gravity have different directions). We suggested a new non-perturbative approximation, which yields a hybrid Eulerian/Lagrangian partial differential equation for the displacement field that contains as limits (i) the Eulerian and Lagrangian linear approximations, and (ii) the standard "adhesion approximation'' if one imposes the assumption of parallelism of peculiar-velocity and -acceleration. When particularized to the "EJN model'', the "virial equation of state'' $p \propto \varrho^2$ provides substantial simplifications and, at the same time, is a reasonable choice. We argued how this approximation can be viewed as the next step beyond the parallelism assumption required by the "adhesion model''. Finally, we provided a qualitative discussion of the exact "EJN'' and "SSE'' models concerning the formation of singularities.

An issue we have not addressed in this work is the choice of the smoothing scale L. This must be guided by the expected applicability of the corresponding closure approximations, in turn likely related to the initial conditions. The investigation of this point seems to require first the ability to extract predictions from the new models in order to validate them with respect to observations and viable N-body simulations.

Further exploration and numerical realization of the models presented will hopefully lead us not only to approaches capable of following large-scale structure into the non-linear regime, but will also offer an interesting way of understanding the role of previrialization (Peebles & Groth 1976; Davis & Peebles 1977; Peebles 1990; okas et al. 1996), and the onset of virialization in connection with galaxy and galaxy cluster formation. Furthermore such models can be employed to improve generic inhomogeneous collapse models (Buchert et al. 2000; Kerscher et al. 2001), and to enhance the power of reconstructions of initial data from present-day density and peculiar-velocity fields (e.g., Croft & Gaztañaga 1997; Brenier et al. 2003; Mohayaee et al. 2003 for recent efforts). In the latter context, the proposed non-perturbative approximation, based on the relation ${\vec{w}}\propto {\vec{P}}$, Eq. (55), may provide an alternative to the assumption of parallelism of peculiar-acceleration and peculiar-velocity, ${\vec{w}}\propto {\vec{u}}$, Eq. (24).

Acknowledgements

This work was motivated by an exciting workshop at Observatoire de la Côte d'Azur (OCA) on the present subject with financial support by OCA, the French Ministry of Education, the Programme National de Cosmologie and the Laboratoire G. D. Cassini. Particular thanks to Uriel Frisch for organizing this workshop. Fruitful discussions were held with him, Erik Aurell, Stéphane Colombi, José Gaite, Sergei Gurbatov, Roman Juszkiewicz, Roya Mohayaee, Sergei Shandarin, Andrei Sobolevski, Alexei Starobinsky, Takayuki Tatekawa, Roland Triay, Patrick Valageas and Barbara Villone. T.B. acknowledges support by the "Sonderforschungsbereich SFB 375 für Astro-Teilchenphysik der Deutschen Forschungsgemeinschaft'', and by CERN, Geneva, where a preliminary manuscript was written during a visit of TB in 1998. A.D. acknowledges financial support by the project "Formación de estructuras en astrofísica y cosmología'' (BFM2002-01014) of the Spanish government.

References

• Adler, S., & Buchert, T. 1999, A&A, 343, 317 [NASA ADS] (In the text)
• Antonov, N. V. 2004, Phys. Rev. Lett., 92, 161101 [NASA ADS] [CrossRef] (In the text)
• Aurell, E., Fanelli, D., Gurbatov, S. N., & Moshkov, A. Yu. 2003, Physica D, 186, 171 [NASA ADS] (In the text)
• Bagla, J. S., & Padmanabhan, T. 1994, MNRAS, 266, 227 [NASA ADS] (In the text)
• Barbero, F., Domínguez, A., Goldman, T., & Pérez-Mercader, J. 1997, Europhys. Lett., 38, 637 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Barrow, J. D., & Maartens, R. 1999, Phys. Rev. D, 59, 043502 [NASA ADS] [CrossRef] (In the text)
• Berera, A., & Fang, L.-Z. 1994, Phys. Rev. Lett., 72, 458 [NASA ADS] [CrossRef] (In the text)
• Bernardeau, F., Colombi, S., Gaztañaga, E., & Scoccimarro, R. 2002, Phys. Rep., 367, 1 [NASA ADS] [CrossRef] (In the text)
• Bharadwaj, S. 1996, ApJ, 472, 1 [NASA ADS] [CrossRef]
• Bildhauer, S., & Buchert, T. 1991, Prog. Theor. Phys., 86, 653 [NASA ADS] (In the text)
• Bildhauer, S., Buchert, T., & Kasai, M. 1992, A&A, 263, 23 [NASA ADS] (In the text)
• Binney, J., & Tremaine, S. 1987, Galactic Dynamics (Princeton Univ. Press) (In the text)
• Bouchet, F. R., Juszkiewicz, R., Colombi, S., & Pellat, R. 1992, ApJ, 394, L5 [NASA ADS] [CrossRef] (In the text)
• Bouchet, F. R., Colombi, S., Hivon, E., & Juszkiewicz, R. 1995, A&A, 296, 575 [NASA ADS] (In the text)
• Brainerd, T. G., Scherrer, R. J., & Villumsen, J. V. 1993, ApJ, 418, 570 [NASA ADS] [CrossRef] (In the text)
• Brenier, Y., Frisch, U., Hénon, M., et al. 2003, MNRAS, 346, 501 [NASA ADS] [CrossRef] (In the text)
• Buchert, T. 1989, A&A, 223, 9 [NASA ADS] (In the text)
• Buchert, T. 1992, MNRAS, 254, 729 [NASA ADS] (In the text)
• Buchert, T. 1993a, Habilitation Thesis (Ludwig-Maximilians-Universität München; link on ADS) (In the text)
• Buchert, T. 1993b, A&A, 267, L51 [NASA ADS] (In the text)
• Buchert, T. 1994, MNRAS, 267, 811 [NASA ADS] (In the text)
• Buchert, T. 1996 in Proc. IOP Enrico Fermi, Course CXXXII (Dark Matter in the Universe), Varenna 1995, ed. S. Bonometto, J. Primack, & A. Provenzale (IOS Press Amsterdam), 543 (In the text)
• Buchert, T. 1999 in Proc. MPG-CAS Workshop, Shanghai 1998 [arXiv:astro-ph/9901002] (In the text)
• Buchert, T., & Bartelmann, M. 1991, A&A, 251, 389 [NASA ADS]
• Buchert, T., & Domínguez, A. 1998, A&A, 335, 395 [NASA ADS] (In the text)
• Buchert, T., Domínguez, A., & Pérez-Mercader, J. 1999, A&A, 349, 343 [NASA ADS] (In the text)
• Buchert, T., & Ehlers, J. 1993, MNRAS, 264, 375 [NASA ADS] (In the text)
• Buchert, T., & Ehlers, J. 1997, A&A, 320, 1 [NASA ADS]
• Buchert, T., Karakatsanis, G., Klaffl, R., & Schiller, P. 1997, A&A, 318, 1 [NASA ADS]
• Buchert, T., Kerscher, M., & Sicka, C. 2000, Phys. Rev. D, 62, 043525-1 [NASA ADS] [CrossRef]-21 (In the text)
• Buchert, T., Melott, A. L., & Weiß, A. G. 1994, A&A, 288, 349 [NASA ADS] (In the text)
• Chandrasekhar, S. 1967, An introduction to the study of stellar structure (Dover Pub.) (In the text)
• Chandrasekhar, S., & Lee, E. S. 1968, MNRAS, 139, 135 [NASA ADS] (In the text)
• Chapman, S., & Cowling, T. G. 1991 The Mathematical Theory of Non-uniform Gases (Cambridge Univ. Press) (In the text)
• Chernin, A. D., Nagirner, D. I., & Starikova, S. V. 2003, A&A, 399, 19 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Clark, R. A., Ferziger, J. H., & Reynolds, W. C. 1979, J. Fluid Mech., 91, 1 [NASA ADS] (In the text)
• Coles, P., Melott, A. L., & Shandarin, S. F. 1993, MNRAS, 260, 765 [NASA ADS] (In the text)
• Croft, R. A. C., & Gaztañaga, E. 1997, MNRAS, 285, 793 [NASA ADS] (In the text)
• Davis, M., & Peebles, P. J. E. 1977, ApJS, 34, 425 [NASA ADS] [CrossRef] (In the text)
• Domínguez, A. 2000, Phys. Rev. D, 62, 103501 [NASA ADS] [CrossRef] (In the text)
• Domínguez, A. 2002, MNRAS, 334, 435 [NASA ADS] [CrossRef] (In the text)
• Domínguez, A. 2003, Astron. Nachr., 324, 560 [NASA ADS] [CrossRef] (In the text)
• Domínguez, A., & Gaite, J. 2001, Europhys. Lett., 55, 458 [EDP Sciences] [CrossRef] (In the text)
• Domínguez, A., Hochberg, D., Martín-García, J. M., Pérez-Mercader, J., & Schulman, L. 1999, A&A, 344, 27 [NASA ADS]; Erratum: A&A, 363, 373 (In the text)
• Domínguez, A., & Melott, A. L. 2004, A&A, 419, 425 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Doroshkevich, A. G., Kotok, E. V., Novikov, I. D., et al. 1980, MNRAS, 192, 321 [NASA ADS] (In the text)
• Ehlers, J., & Buchert, T. 1997, Gen. Rel. Grav., 29, 733 [NASA ADS] [CrossRef] (In the text)
• Ehlers, J., & Rienstra, W. 1969, ApJ, 155, 105 [NASA ADS] [CrossRef] (In the text)
• Fanelli, D., & Aurell, E. 2002, A&A, 395, 399 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Frisch, U. 1995, Turbulence (Cambridge Univ. Press) (In the text)
• Frisch, U., & Bec, J. 2001 in New Trends in Turbulence (Les Houches 2000), ed. M. Lesieur, A. Yaglom, & F. David (Springer EDP Sciences), 341 [arXiv:nlin.CD/0012033] (In the text)
• Götz, G. 1988, Class. Quant. Grav., 5, 743 [NASA ADS] (In the text)
• Gorini, V., Kamenshchik, A., Moschella, U., & Pasquier, V. 2004, [arXiv:gr-qc/0403062] (In the text)
• Gouda, N., & Nakamura, T. 1988, Prog. Theor. Phys., 79, 765 [NASA ADS] (In the text)
• Gouda, N., & Nakamura, T. 1989, Prog. Theor. Phys., 81, 633 [NASA ADS] (In the text)
• Gurbatov, S. N., Malakhov, A. N., & Saichev, A. I. 1991, Nonlinear Random Waves and Turbulence in Nondispersive Media (Manchester Univ. Press) (In the text)
• Gurbatov, S. N., Saichev, A. I., & Yakushkin, I. G. 1983, Sov. Phys. Usp., 26, 857 (In the text)
• Gurbatov, S. N., Saichev, A. I., & Shandarin, S. F. 1985, Sov. Phys. Dokl., 30, 921 (In the text)
• Gurbatov, S. N., Saichev, A. I., & Shandarin, S. F. 1989, MNRAS, 236, 385 [NASA ADS] (In the text)
• Gurevich, A. V., & Zybin, K. P. 1995, Sov. Phys. Usp., 38(7), 687 (In the text)
• Hamana, T. 1998, ApJ, 507, L1 [NASA ADS] [CrossRef] (In the text)
• Haubold, H. J., Mathai, A. M., & Mücket, J. P. 1991, Astron. Nachr., 312, 1 [NASA ADS] (In the text)
• Jones, B. J. T. 1996, in Proc. IOP Enrico Fermi, Course CXXXII (Dark Matter in the Universe), Varenna 1995, ed. S. Bonometto, J. Primack, & A. Provenzale (IOS Press Amsterdam), 661 (In the text)
• Kerscher, M., Buchert, T., & Futamase, T. 2001, ApJ, 558, L79 [NASA ADS] [CrossRef] (In the text)
• Klimas, A. J. 1987, J. Comput. Phys., 68, 202 [NASA ADS] [CrossRef] (In the text)
• Kofman, L. A. 1991 in IUAP Proceedings on Nucleosynthesis in the Universe, ed. K. Sato (Dordrecht: Kluwer), 495 (In the text)
• Kofman, L. A., Bertschinger, E., Dekel, A., Gelb, J. M., & Nusser, A. 1994, ApJ, 420, 44 [NASA ADS] [CrossRef] (In the text)
• Kofman, L. A., Pogosyan, D. Yu., & Shandarin, S. F. 1990, MNRAS, 242, 200 [NASA ADS] (In the text)
• Kofman, L. A., Pogosyan, D. Yu., Shandarin, S. F., & Melott, A. L. 1992, ApJ, 393, 437 [NASA ADS] [CrossRef] (In the text)
• Lions, P. L., & Perthame, B. 1991, Invent. Math., 105, 415 [NASA ADS] [CrossRef] (In the text)
• Lokas, E. L., Juskiewicz, R., Bouchet, F. R., & Hivon, E. 1996, ApJ, 467, 1 [NASA ADS] [CrossRef] (In the text)
• Lynden-Bell, D. 1967, MNRAS, 136, 101 [NASA ADS] (In the text)
• Ma, C.-P., & Bertschinger, E. 2004, ApJ, 612, 28 [NASA ADS] (In the text)
• Maartens, R., Triginer, J., & Matravers, D. R. 1999, Phys. Rev. D, 60, 103503 [NASA ADS] [CrossRef] (In the text)
• Matarrese, S., Lucchin, F., Moscardini, L., & Sáez, D. 1992, MNRAS, 259, 437 [NASA ADS] (In the text)
• Matarrese, S., & Mohayaee, R. 2002, MNRAS, 329, 37 [NASA ADS] (In the text)
• Melott, A. L. 1994, ApJ, 426, L19 [NASA ADS] [CrossRef] (In the text)
• Melott, A. L, Buchert, T., & Weiß, A. G. 1995, A&A, 294, 345 [NASA ADS] (In the text)
• Melott, A. L., Shandarin, S. F., Splinter, R. J., & Suto, Y. 1997, ApJ, 479, L79 [NASA ADS] [CrossRef] (In the text)
• Menci, N. 2002, MNRAS, 330, 907 [NASA ADS] [CrossRef] (In the text)
• Morita, M., & Tatekawa, T. 2001, MNRAS, 328, 815 [NASA ADS] [CrossRef] (In the text)
• Mohayaee, R., Frisch, U., Matarrese, S., & Sobolevskii, A. 2003, A&A, 406, 393 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Padmanabhan, T. 1996, MNRAS, 278, L29 [NASA ADS] (In the text)
• Pauls, J. L., & Melott, A. L. 1995, MNRAS, 274, 99 [NASA ADS] (In the text)
• Peebles, P. J. E. 1980, The Large-scale Structure of the Universe (Princeton Univ. Press) (In the text)
• Peebles, P. J. E. 1990, ApJ, 365, 27 [NASA ADS] [CrossRef] (In the text)
• Peebles, P. J. E., & Groth, E. J. 1976, A&A, 53, 131 [NASA ADS] (In the text)
• Pfaffelmoser, K. 1992, J. Diff. Eqs., 95, 281 [CrossRef] (In the text)
• Pope, S. B. 2000, Turbulent Flows (Cambridge Univ. Press) (In the text)
• Rein, G., & Rendall, A. D. 1994, Arch. Rational Mech. Anal., 126, 183 [CrossRef] (In the text)
• Ribeiro, A. L. B., & Peixoto de Faria, J. G. 2005, Phys. Rev. D, 71, 067302 [NASA ADS] [CrossRef] (In the text)
• Sahni, V., & Coles, P. 1995, Phys. Rep., 262, 1 [NASA ADS] [CrossRef] (In the text)
• Sahni, V., & Sathyaprakash, B. S., & Shandarin, S. F. 1994, ApJ, 431, 20 [NASA ADS] [CrossRef] (In the text)
• Sathyaprakash, B. S., Sahni, V., Munshi, D., Pogosyan, D., & Melott, A. L. 1995, MNRAS, 275, 463 [NASA ADS] (In the text)
• Schaeffer, J. 1991, Part. Diff. Eqns., 16, 1313 (In the text)
• Serrin, I. 1959, in Encyclopedia of Physics Vol. III.1, ed. S. Flügge (Berlin: Springer) (In the text)
• Shandarin, S. F., & Zel'dovich, Ya. B. 1989, Rev. Mod. Phys., 61, 185 [NASA ADS] [CrossRef] (In the text)
• Shu, F. H. 1978, ApJ, 225, 83 [NASA ADS] [CrossRef] (In the text)
• Splinter, R. J., Melott, A. L., Shandarin, S. F., & Suto, Y. 1998, ApJ, 497, 38 [NASA ADS] [CrossRef] (In the text)
• Spohn, H. 1991, Large Scale Dynamics of Interacting Particles (Springer-Verlag) (In the text)
• Tatekawa, T., Suda, M., Maeda, K., Morita, M., & Anzai, H. 2002, Phys. Rev. D, 66, 064014 [NASA ADS] [CrossRef] (In the text)
• Tatekawa, T. 2004a, Phys. Rev. D, 69, 084020 [NASA ADS] [CrossRef] (In the text)
• Tatekawa, T. 2004b, Phys. Rev. D, 70, 064010 [NASA ADS] [CrossRef] (In the text)
• Tatekawa, T. 2005, Phys. Rev. D, 71, 044024 [NASA ADS] [CrossRef] (In the text)
• Trümper, M. 1967, Z. Astrophys., 66, 215 [NASA ADS] (In the text)
• Valageas, P. 2001, A&A, 379, 8 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Valageas, P. 2002, A&A, 382, 477 [EDP Sciences] [NASA ADS] [CrossRef] (In the text)
• Vergassola, M., Dubrulle, B., Frisch, U., & Noullez, A. 1994, A&A, 289, 325 [NASA ADS] (In the text)
• Vreman, B., Geurts, B., & Kuerten, H. 1996, Theor. Comp. Fluid Dyn., 8, 309 (In the text)
• Vreman, B., Geurts, B., & Kuerten, H. 1997, J. Fluid Mech., 339, 357 [NASA ADS] [CrossRef] (In the text)
• Weinberg, D. H., & Gunn, J. 1990, MNRAS, 247, 260 [NASA ADS] (In the text)
• Weiß, A. G., Gottlöber, S., & Buchert, T. 1998, MNRAS, 278, 953 [NASA ADS]
• Yano, T., & Gouda, N. 1998, ApJS, 118, 267 [NASA ADS] [CrossRef] (In the text)
• Yano, T., Koyama, H., Buchert, T., & Gouda, N. 2003, ApJS, 151, 185 [NASA ADS] (In the text)
• Zel'dovich, Ya. B. 1970, A&A, 5, 84 [NASA ADS] (In the text)
• Zel'dovich, Ya. B. 1973, Astrophys., 6, 164 [CrossRef] (In the text)
• Zel'dovich, Ya. B., & Shandarin, S. F. 1982, Sov. Astron. Lett., 8, 139 [NASA ADS]

 

  
Online Material

   
Appendix A: On the transformation to comoving coordinates

To derive Eqs. (13) we first introduced comoving coordinates, Eqs. (10), and then the coarse-graining procedure was applied. It is legitimate to ask whether the order can be exchanged, i.e., what are the transformations that lead from Eqs. (8) to Eqs. (13)? They can be obtained easily by inserting the transformations (10) into the definitions (12), (14) and using the definitions (7), (9). One gets (here we write explicitly the dependence on the smoothing length):

$$% {\vec{q}}= \frac{1}{a} {\vec{x}}, \qquad L = \frac{1}{a} {\cal L},$$ (A.1a)

$$% \varrho ({\vec{q}}, t; L) = \rho ({\vec{x}}, t; {\cal L}),$$ (A.1b)

$$% {\vec{u}}({\vec{q}}, t; L) = {\vec{v}}({\vec{x}}, t; {\cal L}) - H \vec{R} ({\vec{x}}, t; {\cal L}),$$ (A.1c)

$$% {\vec{w}}({\vec{q}}, t; L) = {\vec{g}}({\vec{x}}, t; {\cal L}) + \frac{4 \pi G \varrho_{\rm H} - \Lambda}{3} {\vec{x}},$$ (A.1d)

$$% {\vec{F}}({\vec{q}}, t; L) = \vec{\cal F} ({\vec{x}}, t; {\... ...ambda}{3} (\vec{R}-{\vec{x}}) \rho ({\vec{x}}, t; {\cal L}),$$ (A.1e)

$$% \Pi_{ij} ({\vec{q}}, t; L) = {\cal P}_{ij} ({\vec{x}}, t; {\cal... ...} ({\vec{x}}, t; {\cal L})\right] + H^2 {\cal I}_{ij} ({\vec{x}}, t; {\cal L}).$$ (A.1f)

The following new quantities had to be introduced:

$$% \rho \vec{R} ({\vec{x}}, t; {\cal L}) := \frac{m}{{\cal L}^... ...\vec{x}}^{(\alpha)}}{{\cal L}} \right) {\vec{x}}^{(\alpha)} ,$$ (A.2a)

$$% {\cal D}_{ij} ({\vec{x}}, t; {\cal L}) := \frac{m}{{\cal L}^3} ... ...\right) v_i^{(\alpha)} x_j^{(\alpha)} - \rho v_i R_j ({\vec{x}}, t; {\cal L}) ,$$ (A.2b)

$$% {\cal I}_{ij} ({\vec{x}}, t; {\cal L}) := \frac{m}{{\cal L}^3} ... ... \right) x_i^{(\alpha)} x_j^{(\alpha)} - \rho R_i R_j ({\vec{x}}, t; {\cal L}).$$ (A.2c)

$\vec{R}$ is the center-of-mass position of the subsystem defined by the smoothing window, and  $\vec{\cal I}$ is the inertia tensor of this same subsystem with respect to $\vec{R}$. Concerning  $\vec{\cal D}$, we note that the antisymmetrized combination  $\epsilon_{ijk} {\cal D}_{kj}$ is the angular momentum of this subsystem with respect to $\vec{R}$; the transformation (A.1f) involves only the symmetrized tensor. Therefore, to obtain a self-contained set of transformations, new dynamical equations for $\vec{R}$, $\vec{\cal I}$, $\vec{\cal D}$ would be required, including also further approximations to close the corresponding hierarchy. The simplest one is provided by the "dust model'': the small-scale degrees of freedom are neglected altogether as irrelevant for the dynamics, so that $\vec{R}={\vec{x}}$, $\vec{\cal I}={\vec{0}}$, $\vec{\cal D}={\vec{0}}$, and the transformation (A.1) reduces to the standard one (e.g., Peebles 1980).

Compared to this "standard transformation'', transformation (A.1) differs in two points: (i) the additional transformation rules (A.1e,f) involving the fields ${\vec{F}}$ and  ${\vec{\Pi}}$, and (ii) the transformation between ${\vec{u}}$ and ${\vec{v}}$, because in general $\vec{R} ({\vec{x}}, t) \neq {\vec{x}}$. When the particle distribution looks homogeneous at all cosmological scales, the deviation  $\vec{R}-{\vec{x}}$ is irrelevant. In the highly non-linear regime of large inhomogeneities, one expects $\vert\vec{R} -{\vec{x}}\vert \sim a L$, but then $\vert{\vec{v}}\vert, \vert{\vec{u}}\vert \gg H a L$ nevertheless, and the cosmological expansion is also irrelevant (Domínguez & Gaite 2001). The difference may be important in the most interesting case when the structure is beginning to enter the non-linear regime, $\vert{\vec{v}}\vert, \vert{\vec{u}}\vert \sim H a L$.

Although the change to comoving coordinates does not affect the physical content of the resulting equivalent Eqs. (8) or (13), this equivalence can be broken after approximations are made. It seems more straightforward to "subtract'' the known background cosmological expansion before performing any smoothing/averaging, as we did to derive Eqs. (13), so that the approximations are concerned only with the "unknown'' part of the evolution of the inhomogeneities. However, the above remarks show that this procedure "hides'' additional physics and requires a careful definition of "background'' at the coarse-graining scale, which is only implicit by coarse-graining in the pre-defined comoving space.

   
Appendix B: Derivation of the "Small-Size-Expansion (SSE)''

In this Appendix we recall the derivation of the "SSE'' closure (19) (Domínguez 2000). In analogy with definitions (12), one can formally define microscopic fields by integrating the Klimontovich density (3) over the velocity variable:

$$% \varrho_{{\rm mic}} ({\vec{q}}, t) m \int {\rm d}{\vec{u}}\; f_{\rm K} (a {\vec{q}}, {\vec{u}}-\dot{a} {\vec{q}}, t)$$ :=

$$\frac{m}{a^3} \sum_\alpha \delta \left({\vec{q}}- {\vec{q}}^{(\alpha)}\right),$$ = (B.1a)

$$% \varrho_{{\rm mic}} \bu_{{\rm mic}} ({\vec{q}}, t) m \int {\rm d}{\vec{u}}\; {\vec{u}}~ f_{\rm K} (a {\vec{q}}, {\vec{u}}-\dot{a} {\vec{q}}, t)$$ :=

$$\frac{m}{a^3} \sum_\alpha {\vec{u}}^{(\alpha)} ~ \delta \left({\vec{q}}- {\vec{q}}^{(\alpha)}\right),$$ = (B.1b)

 

$$% \varrho_{{\rm mic}} {\vec{w}}_{{\rm mic}} ({\vec{q}}, t) m \int {\rm d}{\vec{u}}\; {\vec{w}}({\vec{q}}, t) ~ f_{\rm K} (a {\vec{q}}, {\vec{u}}-\dot{a} {\vec{q}}, t)$$ :=

$$\frac{m}{a^3} \sum_\alpha {\vec{w}}^{(\alpha)} ~ \delta \left({\vec{q}}- {\vec{q}}^{(\alpha)}\right),$$ = (B.1c)

already expressed in comoving coordinates. In Eq. (B.1c), ${\vec{w}}({\vec{q}}, t) := {\vec{g}}(a{\vec{q}}, t) + (4 \pi G \varrho_{\rm H}(t) - \Lambda) a {\vec{q}}/3$, with ${\vec{g}}$ given by Eqs. (4b). Note that ${\vec{w}}_{{\rm mic}}({\vec{q}}, t) \equiv {\vec{w}}({\vec{q}},t)$, so that by virtue of Eqs. (4b) (when expressed in comoving coordinates), we can write:

$$% \nabla_{\vec{q}}\cdot {\vec{w}}_{{\rm mic}} = - 4 \pi G a \... ...ad \nabla_{\vec{q}}\times {\vec{w}}_{{\rm mic}} = {\vec{0}}.$$ (B.2)

In terms of the microscopic fields, definitions (12), (14) can be rewritten as follows:

$$% \varrho ({\vec{q}}, t) = \int \frac{{\rm d}{\vec{q}}'}{L^3}... ...{L} \right) \varrho_{{\rm mic}} \left({\vec{q}}', t\right),$$ (B.3a)

$$% \varrho \bar{\vec{u}}({\vec{q}}, t) = \int \frac{{\rm d}{\v... ...arrho_{{\rm mic}} \bu_{{\rm mic}} \left({\vec{q}}', t\right),$$ (B.3b)

$$% {\vec{F}}({\vec{q}}, t) = \int \frac{{\rm d}{\vec{q}}'}{L^3... ...t({\vec{q}}', t\right) - \varrho \bar{\vec{w}}({\vec{q}}, t),$$ (B.3c)

$$% {\vec{\Pi}}({\vec{q}}, t) = \int \frac{{\rm d}{\vec{q}}'}{L^3} ... ... \left({\vec{q}}', t\right) - \varrho \bar{\vec{u}}\bar{\vec{u}}({\vec{q}}, t),$$ (B.3d)

where $\bar{\vec{w}}$ is the mean-field strength of the peculiar-gravity, defined by Eqs. (13c,d). The introduction of the mic-fields is a formal step to express the definitions in terms of volume integrals over fields. Notice in particular that ${\vec{u}}$ is defined by a volume-average of the peculiar-momentum density (i.e., by a mass-average of the velocity) and thus it has the physical meaning of center-of-mass velocity of the coarsening cell. These fields obey Eqs. (B.3a,b), which can be formally inverted to yield the following expansion in L (this is most easily derived in Fourier space; this formal inversion is defined provided the Fourier transform of the window $W (\cdot)$ has no zeros):

$$% \varrho_{{\rm mic}} = \left[ 1 - \frac{1}{2} B ~ \left(L ~ ... ...al O}\left(L ~ \nabla_{{\vec{q}}}\right)^4 \right] ~ \varrho,$$ (B.4a)

$$% \bu_{{\rm mic}} = \left[ 1 - \frac{1}{2} B ~ \left(L ~ \nab... ...ft(L ~ \nabla_{{\vec{q}}}\right)^4 \right] ~ \bar{\vec{u}},$$ (B.4b)

with

$$% \int {\rm d} \vec{z} \; \vec{z} ~ W(\vec{z}) = {\vec{0}}, \qquad B := \frac{1}{3} \int {\rm d} \vec{z} \; z^2 ~ W(\vec{z}).$$ (B.5)

These expansions show how the microscopic fields can be recovered by taking into account finer and finer details in the spatial distribution of the coarse-grained fields. To derive Eq. (19b), one inserts Eqs. (B.4) into the definition (B.3d) encompassing the coupling to the small-scale kinetic degrees of freedom:

$$\begin{eqnarray*}{\vec{\Pi}}({\vec{q}}, t) & = & \int \frac{{\rm d}{\vec{q}}'}{... ...\bar{\vec{u}}\right) ({\vec{q}}, t) + {\cal O}\left(L^4\right), \end{eqnarray*}$$

where the fields in the integrand have been expanded about ${\vec{q}}'={\vec{q}}$consistently to a given order in L. The expansion (19a) for the correction to mean field is derived in the same manner (for this, notice that  ${\vec{w}}_{{\rm mic}}$ is a linear functional of  $\varrho_{{\rm mic}}$ given by Eqs. (B.2)).

   
Appendix C: Dynamical equations of state and corresponding model features

An ingredient of the "EJN'' model is an "equation of state'' $p(\varrho)$ for the trace of the velocity dispersion tensor as function of the mass density, Eq. (16). In this appendix we summarize the evidence in favor of this assumption and the physical meaning of some particular functional dependences.

Numerical simulations

The density-dependence of the trace of the velocity dispersion (14b) was addressed systematically with N-body simulations of CDM scenarios by Domínguez (2003) and Domínguez & Melott (2004) and can be summarized in a (sometimes piecewise) polytropic dependence $p \propto \varrho^\gamma$, which should be understood as an average, because the numerical data scatter around this dependence. The index $\gamma$ is always found to lie in the range $1 \leq \gamma \leq 2$. The more power there is in the small scales of the initial mass distribution relative to the large scales, the smaller is $\gamma$.

Special values of $\gamma$

There are values of the index $\gamma$ which have a particular significance (see also Buchert et al. 1999):

$\bullet$

$\gamma=2$: "Virialized state''

This model is frequently encountered in our paper. From a mathematical point of view, the index $\gamma=2$ entails a simplification because the length  $\lambda(\varrho)$, Eq. (18), becomes $\varrho$-independent. Thus, the standard "adhesion model'' is recovered in this case because the GM coefficient $\mu$, Eq. (26), is also $\varrho$-independent; see also the derivation of Eq. (57).

Physically, one can argue on the basis of the virial theorem (e.g. Chandrasekhar & Lee 1968): a well-known consequence of it is that, if a system dominated by gravity is isolated and in a stationary state, then the total kinetic energy ${\cal U}$ (in the frame of the system's center of mass) and the total (gravitational) potential energy ${\cal W}$ are related by the equation $2 {\cal U} + {\cal W} = 0$. If our system is a coarsening cell of size $\sim$L defined by the smoothing window $W (\cdot)$, we have ${\cal U} \propto p$ and ${\cal W} \sim G L^5 \varrho^2$ assuming that ${\cal W}$ can be estimated reliably within a mean field approximation (i.e., not too much small-scale structure). To the extent that a coarsening cell can be considered approximately isolated and stationary, the virial theorem implies $p \propto \varrho^2$. The use of this relationship in the evolution Eq. (17b) always supposes an approximation, because the term $-\nabla p$ represents the flux of momentum by the exchange of particles between neighboring coarsening cells, in contradiction with the assumption of isolation required to derive the dependence $p \propto \varrho^2$.

$\bullet$

$\gamma=1$: "Isothermal state''

In this case, the plane-symmetric problem of the "EJN model'', in the absence of cosmological expansion, can be mapped exactly to the "Sine-Gordon equation'' (Götz 1988).

Physically, this dependence implies that the kinetic energy per particle, which is proportional to $p/\varrho$, is density-independent. In the context of a fluid in thermal equilibrium, the reason is that, concerning any single particle, the rest of the system behaves like a thermostat.

$\bullet$

$\gamma=5/3$: "Adiabatic state''

Starting from the definition (14b), one can compute the following evolution equation for the tensor $\Pi_{ij}$:

 

$$% \frac{\partial \Pi_{ij}}{\partial t} + 5 H \Pi_{ij} - \frac{1}{a} \Pi_{ik} ~ \left(\partial_k u_j\right) - \frac{1}{a} \Pi_{jk} ~ (\partial_k u_i)$$ =

$$- \frac{1}{a} \partial_k \left(u_k ~ \Pi_{ij}\right) - \frac{1}{a} \partial_k {\cal L}_{ijk} + {\cal P}_{ij}.$$ (C.1)

In this equation, ${\cal L}_{ijk}$ is the third reduced velocity moment (representing the heat flux in dilute gases) and  ${\cal P}_{ij}$ involves the departures from mean field. (The precise expressions for these two tensors can be found in Domínguez 2000.) One may take the closure assumption that these two tensors are negligible, meaning physically that corrections to mean field are unimportant, and that velocity dispersion is relatively small (and with it the "heat flux''  $\Rightarrow$ "adiabatic evolution''). If one further makes the isotropy ansatz $\Pi_{ij} = p \delta_{ij}$, this equation in combination with the continuity Eq. (13a) then yields $p = \kappa \varrho^{5/3}$ and the strong kinematical constraint of vanishing shear, $\partial_i u_j + \partial_j u_i \propto \delta_{ij}$ (Buchert & Domínguez 1998). Note that this implies that the only consistent exact solution is spatially homogeneous (Trümper 1967; Ehlers & Rienstra 1969; the homogeneous model also satisfies "adiabaticity''.)

The coefficient of proportionality $\kappa$ is a local function of the Lagrangian coordinates given by the initial data; if we assume that, initially, this function is a spatial constant, then we are entitled to talk about a "state'', since the relation between dynamical "pressure'' and mass density is then global.

$\bullet$

$\gamma=3$: Plane symmetric "adiabatic state''

In a plane-symmetric configuration, Eq. (C.1) can be integrated under the assumption that  ${\cal L}_{ijk}$ and  ${\cal P}_{ij}$ are negligible without the need for the additional isotropy ansatz for $\Pi_{ij}$. Thus, the equation for the single component $\Pi_{xx}$ becomes

$$\frac{\partial \Pi_{xx}}{\partial t} + \frac{1}{a} u \frac{\p... ...} = - \frac{3}{a} \Pi_{xx} ~ \frac{\partial u}{\partial x},$$

from which one obtains $\Pi_{xx} \propto \varrho^3$.

Without regard to this derivation, in the context of the generalized "adhesion models'', the $\varrho$-dependence of the GM coefficient derived with the "SSE model'', Eq. (28), can be obtained with the "EJN model'' and a dependence $p \propto \varrho^3$, Eq. (26).

$\bullet$

$\gamma=4/3$: "Cosmogenetic state''

If one looks for a solution of Eqs. (17) in "comoving hydrostatic equilibrium'', i.e., ${\vec{u}}={\vec{0}}$ but $\nabla \varrho \neq {\vec{0}}$, with an "equation of state'' $p = \kappa \varrho^\gamma$, then one can easily check that $\gamma=4/3$ is the only possible index for which a solution exists when $\kappa$ is time-independent (Chandrasekhar 1967).

$\bullet$

$\gamma=-1$: "Chaplygin state''

In our context, this formal case is of interest because the particularization of the "EJN model'' to a plane-symmetric configuration gives a linear equation for the Lagrangian displacement field ${\vec{P}}$ exactly, see Sect. 6.2.

More generally, the index $\gamma=-1$ introduces mathematical simplifications in the integration of Euler's equation. Recently, the Chaplygin gas has received attention in cosmological research as one possible model of dark energy (see e.g. the recent preprint by Gorini et al. 2004, and refs. therein).

$\bullet$

$\gamma<1$:

Application of boundary-layer theory to Eqs. (26) (Domínguez 2000) shows that density singularities are regularized when $\gamma>1$. On the other hand, a similar analysis (Domínguez, unpublished) demonstrates that "pressure'' is too weak to avoid the occurrence of singularities when $\gamma<1$.