text:
# The growth of structure in the intergalactic medium

## 1 Introduction

The analysis of the spectra of distant QSOs blueward of the Ly$`\alpha `$ emission, the so-called Ly$`\alpha `$ forest, has revealed the presence of large coherent structures in the cosmic distribution of the neutral hydrogen, left over by the reionization process (e.g. Rauch 1998 and references therein). The low-column density ($`N_{\mathrm{HI}}10^{14.5}`$ cm<sup>-2</sup>) Ly$`\alpha `$ forest is thought as being generated by unshocked gas in voids and mildly overdense fluctuations of the photoionized intergalactic medium (IGM). The IGM is, in turn, expected to trace the underlying dark matter (DM) distribution on large scales, while the thermal pressure in the gas smears small-scale DM nonlinearities. According to this picture, the large coherent spatial structures in the IGM, which are probed by quasar absorption spectra, reflect the underlying network of filaments and sheets in the large-scale DM distribution (e.g. Efstathiou, Schaye & Theuns 2000), the so-called ‘cosmic web’ (Bond, Kofman & Pogosyan 1996). Thanks to this simple picture, a number of semi-analytical, or pseudo-hydrodynamical particle-mesh codes have been developed, which all aim at simulating the IGM distribution, with enough resolution to predict the main features of the Ly$`\alpha `$ forest (e.g. Bi & Davidsen 1997; Croft et al. 1998; Gnedin & Hui 1996, 1998; Hui, Gnedin & Zhang 1997; Petitjean, Mücket & Kates 1995; Gaztañaga & Croft 1999; Meiksin & White 2000; Viel et al. 2001). This physical interpretation of the Ly$`\alpha `$ forest is, moreover, largely corroborated by numerical simulations based on fully hydrodynamical codes (e.g. Cen et al. 1994; Zhang, Anninos & Norman 1995, 1997; Miralda-Escudé et al. 1996; Hernquist et al. 1996; Theuns et al. 1998).

Although the above scenario captures the main physical processes which determine the coarse-grained dynamics of the baryons, it neglects a number of fine-grained details which distinguish the baryon distribution from the underlying DM one. There has been, recently, growing observational and theoretical interest on the low-column density Ly$`\alpha `$ forest, as an important indicator of the state of the Universe in the redshift range $`24`$. We then believe it has become necessary to develop new and more sophisticated techniques, enabling to understand the IGM dynamics and accurately simulate its clustering, without resorting to the full machinery of hydro-simulations. The model presented here aims at providing an accurate scheme for the growth and evolution of mildly nonlinear structures in the baryon gas, accounting for the detailed thermal history of the IGM.

The stochastic adhesion model introduced here is based on applying the forced Burgers equation of nonlinear diffusion (Forster, Nelson & Stephen 1977; Kardar, Parisi & Zhang 1986; Barabási & Stanley 1995; E et al. 1997; Frisch & Bec 2000 and references therein) to the cosmological framework, where it can be given the general form

$$\frac{D𝐮(𝐱,\tau )}{D\tau }=\nu ^2𝐮(𝐱,\tau )\eta (𝐱,\tau ).$$
(1)
Here $`𝐱`$ are comoving coordinates, the time variable $`\tau `$ is chosen to coincide with the cosmic scale-factor $`a`$ (although, for the pure DM evolution, a better choice would be the linear growth factor of density fluctuations), $`𝐮\mathrm{𝐝𝐱}/𝐝\tau `$ is the (suitably rescaled) peculiar velocity field, which is here assumed to be irrotational, $`𝐮=\mathrm{\Phi }`$, and $`D/D\tau /\tau +𝐮`$ denotes the convective time-derivative. The coefficient of kinematical viscosity, $`\nu `$, is always assumed to be small. In standard applications of the forced Burgers equation, $`\eta `$ is a random external potential. For a general cosmological fluid, such as the dark matter or the baryons, $`\eta =(3/2a)(\mathrm{\Phi }+\phi +w)`$, where $`\phi `$ is the local gravitational potential (up to an appropriate rescaling, whose precise form will be given in Section 2), which must be consistently determined via the Local Poisson equation, and $`w`$ is the specific enthalpy of the fluid (up to the same rescaling), which vanishes in the DM case. Therefore, our potential $`\eta `$ in the RHS of eq. (1) is far from being ‘external’, as it is, for instance, dynamically related to the velocity field which appears on the LHS of the same equation.

The main idea of this paper is that, if $`\eta `$ is given an approximate expression, e.g. by using the results of perturbation theory, it can be legitimately treated as a truly external random potential. The choice of variables in eq. (1) is such that, if this approximation technique is applied to the DM evolution, to first order it reduces to the well-known adhesion model (Gurbatov, Saichev & Shandarin, 1985, 1989; Kofman & Shandarin 1988), which was introduced in cosmology to extend the validity of the Zel’dovich approximation (Zel’dovich 1970) beyond the epoch of first caustic formation. In the DM case, therefore, only a second-order calculation would produce a non-vanishing external force, whose presence would then affect rather small scales. Quite different is the case of the collisional baryon component, where, already at first order in perturbation theory (both Eulerian and Lagrangian) a non-zero contribution to $`\eta `$ generally comes out: it is originated from the unbalanced composition of Hubble drag, local gravity and gas pressure, thus representing a genuine baryonic feature.

The presence of the kinematical viscosity term requires some explanation. Its role in the present context is twofold. First, it prevents the formation of multi-streams, which are well-known to affect the DM dynamics, but may also appear as a spurious effect in the collisional case, when pressure gradients are given an approximate form in terms of linear theory. Second, it allows to transform the problem into a linear one, through the so-called Hopf-Cole substitution (e.g. Burgers 1974) $`𝐮=2\nu _𝐱\mathrm{ln}𝒰`$, where the ‘expotential’ $`𝒰`$ obeys the random heat, linear diffusion, equation

$$\frac{𝒰(𝐱,\tau )}{\tau }=\nu ^2𝒰(𝐱,\tau )+\frac{\eta (𝐱,\tau )}{2\nu }𝒰(𝐱,\tau ),$$
(2)
whose solution is expressible in terms of path-integrals (e.g. Feynman & Hibbs 1965). As we will see in Sections 5 and 6, a suitable approximation technique, valid in the small $`\nu `$ case, allows to give the velocity field a simple and finite form, more convenient for practical applications.

The forced Burgers equation has been used to describe a variety of different physical problems, ranging from interfacial growth in condensed matter physics, where it is known as the KPZ model (Kardar, Parisi & Zhang 1986), to fully developed turbulence (e.g. Bouchaud, Mézard & Parisi 1995). Later in the paper we will come back to this interesting connection and discuss both the analogies and the peculiarities of the cosmological application of this equation.

The stochastic adhesion approximation provides an analytical description of the generation, and subsequent merging of shocks, which give rise to the thin network of filaments and sheets in the IGM spatial distribution. Their existence is clearly observed in the spectra of high-redshift QSOs, e.g. through the presence of common absorption lines in the Ly$`\alpha `$ forest of multiple QSOs, with lines of sight separated by several comoving Mpc at redshifts $`z24`$ (e.g. Rauch 1998 and references therein). An important property of our model is that it allows to draw the skeleton of the IGM distribution through a straightforward extension of the geometrical technique applied in the free adhesion model (e.g. Sahni & Coles 1995 and references therein). The present paper will be mostly devoted to provide the physical and mathematical bases for our stochastic adhesion model. Simulations of the IGM large-scale structure will be obtained only from the simplified inviscid ($`\nu =0`$) model. In a subsequent paper we will implement our algorithm to produce numerical simulations of the IGM distribution and to study the statistical properties of the IGM density and velocity fields. Let us stress that approximation schemes like the present one can be particularly useful, as they allow to better account for the cosmic variance of large-scale modes, which is poorly probed (especially at low redshifts) by the existing hydro-simulations, which are forced to adopt small computational boxes to increase the small-scale resolution \[e.g. the discussion in (Viel et al. 2001)\]. Even more interesting is the possibility to combine our scheme with a hydro-code, using the former to provide the large-scale skeleton of the IGM and the latter to achieve the required resolution on small scales.

The approach most closely related to ours is that recently proposed by Jones (1996, 1999), which aims at modelling the nonlinear clustering of the baryonic material. In Jones’ model, however, a different set of variables is adopted, which does not allow a direct comparison neither with the present scheme, nor with the Zel’dovich and the adhesion approximations, in the collisionless limit. The most important difference is that, in Jones’ model, the external random term is identified with the local gravitational potential $`\phi `$, which is treated as an ‘external’ one, as it is essentially generated by the dominant DM component; moreover, no explicit account for the gas pressure is given. In our model instead, the external potential $`\eta `$ is obtained by linearly approximating the composition of Hubble drag, local gravity and thermal pressure which act on the IGM fluid elements; in our case $`\eta `$ is an external potential because it is determined by a convolution of the initial gravitational potential with the IGM linear filter.

It should be clear from this introduction that the forced Burgers equation might have wider applications in the cosmological structure formation problem. Its relevance (in terms of the closely related random heat equation) in the cosmological framework has been first advocated by Zel’dovich and collaborators in the mid eighties (Zel’dovich et al. 1985, 1987), as a means to describe the possible origin of intermittency in the matter distribution. The intermittency phenomenon consists in the appearance of rare high peaks, where most of the matter is concentrated, separated by vast regions of reduced intensity. From the statistical point of view, intermittency in a stochastic process is signalled by an anomalous scaling of e.g. structure functions (moments of velocity increments): moments of order $`p>2`$, made dimensionless by suitable powers of the second-order moment, grow without bound on small scales (e.g. Gärtner J., Molchanov 1990; Frisch 1995). This may be viewed as increased non-Gaussianity on small scales, which, in Fourier space, appears as a slow decrease in the amplitude of Fourier modes with increasing wavenumber and by a definite phase relation between them (Zel’dovich et al. 1985, 1987). Actually, the occurrence of this form of intermittency, which seems too extreme and far from our present understanding of the large-scale structure of the Universe, is usually obtained under special properties of the noise and only appears at asymptotically late times. Much more interesting for the structure formation problem is a second phenomenon, called intermediate intermittency, also described by the forced Burgers equation, which consists in the formation of a cellular, or network structure, with “thin channels of raised intensity (the rich phase), separating isolated islands of the poor phase” (Zel’dovich et al. 1985). This second phenomenon is expected to arise as an intermediate asymptotic situation.

Can one take advantage of the description of the structure formation process in terms of intermediate intermittency to predict the specific non-Gaussian statistics which characterizes the nonlinear density field? This is, we believe, a challenging issue, which would deserve further analysis.

The prototype distribution that describes intermittency is the Lognormal one, which naturally arises in multiplicative processes, through the action of the Central Limit Theorem (e.g. Shimizu & Crow 1988). In the DM case, Coles & Jones (1991) proposed that a local Lognormal mapping of the linear density field can describe the nonlinear evolution of structures in the Universe. Detailed comparison with N-body simulations showed that this model fits very well the bulk of the probability density function (PDF) of the mass density field in N-body simulations, for moderate values of the rms overdensity $`\sigma `$ (e.g. Bernardeau & Kofman 1995). Where the Lognormal fails is in reproducing the correct PDF for the strong clustering regime ($`\sigma 1`$), as well as the high- and low-density tails of the PDF, even on mildly nonlinear scales. Moreover, the predicted spatial pattern is too clumpy and poorly populated of extended structures to reproduce that of simulations (Coles, Melott & Shandarin 1993). A ‘skewed’ Lognormal PDF is proposed by Colombi (1994), to follow the transition from the weakly to the highly nonlinear regime.

Quite noticeably, the Lognormal model has also been used in connection with the IGM dynamics. Indeed, the semi-analytical model proposed by Bi and collaborators (Bi 1993; Bi et al. 1995; Bi & Davidsen 1997; see also Feng & Fang 2000; Roy Choudhury, Padmanabhan & Srianand 2000; Roy Choudhury, Srianand & Padmanabhan 2000; Viel et al. 2001), to simulate the low-column density Ly$`\alpha `$ forest, is based on a local Lognormal model, similar to the one of Coles & Jones (1991). Comparison with the results of more refined techniques has shown that it provides a good fit of the column density distribution in a wide range of values, but it tends to underestimate the abundance of lower column density systems (Hui, Gnedin & Zhang 1997). Moreover, the Lognormal model for the IGM tends to produce an excess of saturated absorption lines in the simulated transmitted flux, compared with real QSO spectra.

The plan of the paper is as follows. In Section 2 we review the equations which govern the Newtonian dynamics of a two-component fluid of dark matter (DM) and baryons, in the expanding Universe; these are solved in Section 3 at the Eulerian linear level and under fairly general assumptions on the baryon equation of state. This allows to obtain the IGM filter connecting the linear baryon density fluctuations to the DM ones. A simplified version of our model is presented in Section 4: it enables one to follow the combined dark matter and baryon dynamics on weakly nonlinear scales, within the laminar regime. In Section 4 we also give the first-order Lagrangian solution for the baryon dynamics, which is then used to perform numerical simulations of the IGM distribution. In Section 5 we discuss the modifications introduced in the baryon dynamics by our improved final model, where we add a kinematic viscosity term, to avoid the occurrence of shell-crossing singularities (arising also in the baryon fluid, because of the approximate treatment of pressure gradients). This leads to our stochastic adhesion model for the dynamics of the intergalactic medium. According to this model, the approximate IGM dynamics is governed by the forced Burgers equation for the baryon peculiar velocity field, whose solution can be expressed as a path-integral. In the physically relevant limit of small viscosity, a solution can be found through the standard saddle-point technique. In Section 6 we discuss how to implement our solution in terms of the first-order Lagrangian particle trajectories previously obtained. A geometrical algorithm is also outlined, which allows to draw the skeleton of the IGM distribution, given a realization of the gravitational potential and the linear IGM filter. A preliminary analysis of the statistical properties of the density field obtained through the stochastic adhesion model is given in Section 7. The concluding Section 8 contains a brief discussion on the possible applications of our model.

## 2 Dynamics of dark matter and baryons in the expanding Universe

The Newtonian dynamics of a self-gravitating two-component fluid, made of collisionless dark matter and collisional baryonic gas is governed by the continuity, Euler and Poisson equations. The continuity equation for the dark matter component (indicated by a subscript $`\mathrm{DM}`$) reads

$$\frac{\rho _{\mathrm{DM}}}{t}+3H\rho _{\mathrm{DM}}+\left(\rho _{\mathrm{DM}}𝐯_{\mathrm{DM}}\right)=0,$$
(3)
where $`\rho `$ is the mass density, $`𝐯ad𝐱/dt`$ the peculiar velocity and $`H`$ the Hubble parameter at time $`t`$; the Euler equation reads

$$\frac{(a𝐯_{\mathrm{DM}})}{t}+\left(𝐯_{\mathrm{DM}}\right)𝐯_{\mathrm{DM}}=\varphi ,$$
(4)
where $`\varphi `$ is the peculiar gravitational potential. For the baryon fluid (indicated by a subscript $`b`$), we have

$$\frac{\rho _b}{t}+3H\rho _\mathrm{b}+\left(\rho _\mathrm{b}𝐯_\mathrm{b}\right)=0,$$
(5)
and

$$\frac{(a𝐯_\mathrm{b})}{t}+\left(𝐯_\mathrm{b}\right)𝐯_\mathrm{b}=\varphi \frac{p_\mathrm{b}}{\rho _\mathrm{b}},$$
(6)
where $`p`$ is the pressure. The peculiar gravitational potential obeys the Poisson equation $`^2\varphi =4\pi Ga^2\delta \rho `$, where the mass-density fluctuation $`\delta \rho `$ takes contribution both from dark matter and baryons. Let then $`f_{\mathrm{DM}}`$ and $`f_\mathrm{b}=1f_{\mathrm{DM}}`$ be the mean mass fraction of these two types of matter. We have

$$^2\varphi =\frac{3}{2}H_0^2\mathrm{\Omega }_{0\mathrm{m}}\left(f_{\mathrm{DM}}\delta _{\mathrm{DM}}+f_b\delta _\mathrm{b}\right),$$
(7)
where $`H_0`$ is the Hubble constant, $`\mathrm{\Omega }_{0\mathrm{m}}`$ is the closure density of non-relativistic matter (both dark matter and baryons) today, $`\delta _{\mathrm{DM}}`$ and $`\delta _\mathrm{b}`$ are the fractional DM and baryon overdensities. In the analytical calculations which follow we will neglect, for simplicity, the self-gravity of the baryons, i.e. we will put $`f_{\mathrm{DM}}=1`$.

To close the system we need the IGM equation of state, which can be taken of the polytropic form

$$p_b=\frac{\rho _\mathrm{b}k_BT}{\mu m_p}=\frac{\overline{\rho }_\mathrm{b}k_BT_0(1+\delta _\mathrm{b})^\gamma }{\mu m_p},$$
(8)
where $`k_B`$ is Boltzman’s constant, $`\gamma `$ the adiabatic index, $`\mu `$ the mean molecular weight (for fully ionized gas with primordial abundances it is about $`0.6`$), $`m_p`$ the proton mass and $`\overline{\rho }_\mathrm{b}`$ the baryon mean density. In writing the pressure term we have assumed the power-law temperature-density relation $`T=T_0(1+\delta _b)^{\gamma 1}`$ where $`T_0=T_0(z)`$ is the IGM temperature at mean density, at redshift $`z`$ (e.g. Hui & Gnedin 1997; Schaye et al. 1999; McDonald et al. 2000). This is adequate for low to moderate baryon overdensity ($`\delta _\mathrm{b}<10`$), where the temperature is locally determined by the interplay between photoheating by the UV background and adiabatic cooling due to the Universe expansion. The underlying assumption for this ‘equation of state’ is a tight local relation between the temperature and the baryon density, which is only true for unshocked gas (e.g. Efstathiou et al. 2000).

In what follows it will prove convenient to change time variable from the cosmic time $`t`$ to the scale-factor $`a`$ (e.g. Shandarin & Zel’dovich 1989; Matarrese et al. 1992; Sahni & Coles 1995). This defines new peculiar velocity fields $`𝐮d𝐱/da=𝐯/(a^2H)`$. Let us also introduce the dimensionless comoving densities $`\rho /\overline{\rho }=1+\delta `$ and a scaled gravitational potential $`\phi =2\varphi /(3a^3H^2)`$. For the redshift range of interest here ($`z24`$), one can write $`a^3H^2H_0^2\mathrm{\Omega }_{0\mathrm{m}}`$. A more exact treatment of the DM component would require the use of the growing mode of linear density perturbations, $`D_+(t)`$, as time variable (e.g. Gurbatov et al. 1989; Catelan et al. 1995); once again, for the range of redshifts of interest here one can safely assume $`D_+(t)a(t)`$. In the Einstein-de Sitter case the present treatment becomes exact.

We have the following set of equations for the DM component

$$\frac{D𝐮_{\mathrm{DM}}}{Da}=\frac{3}{2a}(𝐮_{\mathrm{DM}}+\phi );\frac{D\delta _{\mathrm{DM}}}{Da}=(1+\delta _{\mathrm{DM}})𝐮_{\mathrm{DM}},$$
(9)
where $`D/Da/a+𝐮`$ denotes the convective, or Lagrangian, derivative w.r.t. our new time variable $`a(t)`$. For the baryons, we have

$`{\displaystyle \frac{D𝐮_b}{Da}}`$ $`=`$ $`{\displaystyle \frac{3}{2a}}\left[𝐮_b+\phi +\left({\displaystyle \frac{2\gamma k_BT_0}{3H_0^2\mathrm{\Omega }_{0\mathrm{m}}\mu m_p}}\right)(1+\delta _\mathrm{b})^{\gamma 2}\delta _\mathrm{b}\right]`$ (10)
$`{\displaystyle \frac{D\delta _\mathrm{b}}{Da}}`$ $`=`$ $`(1+\delta _\mathrm{b})𝐮_\mathrm{b}.`$ (11)
The Poisson equation also gets simplified:

$$^2\phi =\frac{\delta _{\mathrm{DM}}}{a}.$$
(12)
## 3 Linear theory

### 3.1 Dark matter in the linear regime

Let us start by writing the above set of equations for the DM component in the (Eulerian) linear approximation. The continuity and Euler equations, respectively, simplify to

$$\dot{\delta }_{\mathrm{DM}}=𝐮_{\mathrm{DM}};\dot{𝐮}_{\mathrm{DM}}=\frac{3}{2a}(𝐮_{\mathrm{DM}}+\phi ),$$
(13)
where, from now on, dots will denote partial differentiation w.r.t. the scale-factor $`a`$. The solutions are well known, and we will simply report them here. Keeping only the growing mode terms we have

$$\delta _{\mathrm{DM}}(𝐱,a)=a\delta _0(𝐱);\phi (𝐱,a)=\phi _0(𝐱);𝐮_{\mathrm{DM}}(𝐱,a)=\phi _0(𝐱),$$
(14)
where $`\delta _0=^2\phi _0`$.

### 3.2 Baryons in the linear regime

In order to make a similar analysis for the baryons we need to specify the time (or redshift) dependence of the baryon mean temperature $`T_0`$, which generally depends upon the thermal history of the IGM, as well as on the spectral shape of the UV background.<sup>1</sup><sup>1</sup>1Similar reasonings would actually also apply to the adiabatic index $`\gamma `$, which we will, however, approximate as being constant in what follows. We will here consider a simple, but fairly general law of the type $`T_0(z)(1+z)^\alpha `$ which will allow to get exact solutions of the linearized baryon equations.

At high redshifts, before decoupling, the mean baryon temperature drops like $`T_0(1+z)`$; when Compton scattering becomes inefficient adiabatic cooling of the baryons implies $`T_0(1+z)^2`$, which makes it practically vanish before reionization \[this will justify our initial conditions for the evolution of the baryon overdensity in eqs. (18)\]. As the Universe reionizes, the IGM temperature rises and a different redshift dependence takes place. Various types of dependence have been considered in the literature. The linear law $`T_0(1+z)`$ is often assumed for simplicity. According to Miralda-Escudé & Rees (1994) and Hui & Gnedin (1997), long after hydrogen reionization has occurred the diffuse IGM settles into an asymptotic state where adiabatic cooling is balanced by photoheating, leading to the power-law $`T_0(1+z)^\alpha `$, with $`\alpha =1/1.76`$. Much steeper an exponent, $`\alpha =1.7`$, has been recently adopted by Bryan & Machacek (2000) and Machacek et al. (2000). More complex is the picture which emerges from hydrodynamical simulations of the IGM, where the mean gas temperature appears to retain some memory of when and how it was reionized (e.g. Schaye et al. 2000). Observational constraints on $`T_0(z)`$ should also be taken into account (e.g. Shaye et al. 1999, 2000; Bryan & Machacek 2000; Ricotti et al. 2000; McDonald et al. 2000).

In view of this variety of assumptions and results it seems reasonable to look for solutions of our equations assuming a general exponent $`\alpha `$ (although, in practical applications we will focus on cases with $`\alpha 1`$). Moreover, as we will see later, our model can be straightforwardly extended to any redshift dependence of the mean IGM temperature.

Solutions of the hydrodynamical equations for general values of $`\alpha `$ have been obtained by Bi, Börner and Chu (1992). Many authors (Peebles 1984; Soloveva & Starobinskii 1985; Shapiro, Giroux & Babul 1994; Nusser 2000) have given solutions for the case $`\alpha =1`$. In particular, Nusser (2000) has obtained the linear IGM overdensity in the case $`\alpha =1`$, extending his calculation to the case of non-negligible baryon fraction (i.e. for $`0f_\mathrm{b}1`$), i.e. accounting for the baryon self-gravity. We will give here an extensive presentation of this problem for general values of $`\alpha `$, both because we are going to use the linear baryon overdensity in our nonlinear model for the IGM dynamics, and because of the specific form taken by the solutions for our set of initial conditions, which had not been previously obtained in the literature. The particular case $`\alpha =1`$ needs to be studied separately; this will be done in the next subsection.

Let us start by writing the linearized continuity and Euler equations for the baryonic component,

$$\dot{\delta }_\mathrm{b}=𝐮_\mathrm{b};\dot{𝐮}_\mathrm{b}=\frac{3}{2a}(𝐮_\mathrm{b}+\phi +\frac{2\gamma k_BT_0}{3H_0^2\mathrm{\Omega }_{0\mathrm{m}}\mu m_p}\delta _\mathrm{b}).$$
(15)
Combining these equations together, using the linear solutions for the DM and our temperature-redshift relation we obtain, in Fourier space,

$$\ddot{\delta }_\mathrm{b}(𝐤,a)+\frac{3}{2a}\dot{\delta }_\mathrm{b}(𝐤,a)+\frac{3}{2a^{\alpha +1}}\frac{k^2}{\stackrel{~}{k}_J^2}\delta _\mathrm{b}(𝐤,a)=\frac{3}{2a}\delta _{\mathrm{DM}}(𝐤).$$
(16)
The redshift-independent wavenumber $`\stackrel{~}{k}_J`$ is related to the comoving Jeans wavenumber $`k_J`$ through

$$k_J^2=a^{\alpha 1}\stackrel{~}{k}_J^2=\frac{3H_0^2\mathrm{\Omega }_{0\mathrm{m}}\mu m_p(1+z)}{\gamma k_BT_0}.$$
(17)
Only for $`\alpha =1`$ the two wavenumbers coincide and the Jeans length becomes a constant.

Equation (16) will be solved with the initial, or, more precisely, matching conditions at $`a=a_i`$,

$$\delta _\mathrm{b}(𝐤,a_i)=\delta _{\mathrm{DM}}(𝐤,a_i);\dot{\delta }_\mathrm{b}(𝐤,a_i)=\dot{\delta }_{\mathrm{DM}}(𝐤,a_i)=\frac{\delta _{DM}(𝐤,a_i)}{a_i}=\frac{\delta _{DM}(𝐤,a)}{a},$$
(18)
which are appropriate if the IGM undergoes sudden reionization at $`z_i`$ (Nusser 2000).

#### 3.2.1 Case $`T_0(z)(1+z)`$

In the case $`\alpha =1`$, we try a solution of the type $`\delta _\mathrm{b}=Aa^n+Ba`$, where $`n`$ and $`B`$ can be found by substitution back into equation (16). This gives

$$\delta _\mathrm{b}(𝐤,a)=a^{1/4}\left[\stackrel{~}{A}_1a^{Q/4}+\stackrel{~}{A}_2a^{Q/4}\right]+\frac{\delta _{\mathrm{DM}}(𝐤,a)}{(1+k^2/k_{J}^{}{}_{}{}^{2})},$$
(19)
where $`Q=\sqrt{124k^2/k_J^2}`$ and $`\stackrel{~}{A}_1`$ and $`\stackrel{~}{A}_2`$ are integration constants. From the latter expression we see that on large scales the homogeneous part only contains decaying modes; on smaller scales, instead, the homogeneous part is characterized by an oscillatory behavior with decaying amplitude. Asymptotically in time one recovers the well known solution (e.g. Peebles 1993) $`\delta _\mathrm{b}(𝐤,a)\delta _{\mathrm{DM}}(𝐤,a)/(1+k^2/k_J^2)`$.

In spite of the absence of growing modes, the homogeneous part of our solution does not necessarily become negligible at the time scales relevant to our problem. It is therefore important, to exactly evaluate the homogeneous part, by relating it to our initial conditions above. For the two constants of integration we find $`\stackrel{~}{A}_1=\stackrel{~}{B}_1\delta _{\mathrm{DM}}(𝐤,a_i)`$ and $`\stackrel{~}{A}_2=\stackrel{~}{B}_2\delta _{\mathrm{DM}}(𝐤,a_i)`$, where

$$\stackrel{~}{B}_1=\frac{5+Q}{2Q}\left(\frac{1}{1+k^2/k_J^2}\right)\left(\frac{k}{k_J}\right)^2;\stackrel{~}{B}_2=\frac{5Q}{2Q}\left(\frac{1}{1+k^2/k_J^2}\right)\left(\frac{k}{k_J}\right)^2.$$
(20)
The baryon peculiar velocity immediately follows from the linear continuity equation. We obtain

$$𝐮_\mathrm{b}=\left[1+\left(\frac{k}{k_J}\right)^2\left(\frac{a}{a_i}\right)^{5/4}\left(\frac{5+Q}{2Q}\left(\frac{a}{a_i}\right)^{Q/4}\frac{5Q}{2Q}\left(\frac{a}{a_i}\right)^{Q/4}\right)\right]\frac{𝐮_{\mathrm{DM}}}{1+k^2/k_J^2}.$$
(21)
#### 3.2.2 General case: $`T_0(z)(1+z)^\alpha `$

We will now look for the solution of eq. (16), for arbitrary values of $`\alpha `$, excluding the case $`\alpha =1`$. In order to find the full solution to the inhomogeneous equation (16), we first notice that its homogeneous counterpart is solved in terms of Bessel functions, namely, $`\delta _{\mathrm{hom}}(𝐤,a)=y^{\stackrel{~}{\alpha }}[A_1𝒥_{\stackrel{~}{\alpha }}(y)+A_2𝒥_{\stackrel{~}{\alpha }}(y)]`$, where $`\stackrel{~}{\alpha }=1/[2(1\alpha )]`$ and we introduced the independent variable $`y\sqrt{24}\stackrel{~}{\alpha }(k/\stackrel{~}{k}_J)a^{1/4\stackrel{~}{\alpha }}=\sqrt{24}\stackrel{~}{\alpha }k/k_J`$. The full solution of eq. (16) for any $`\alpha `$, obtained by the standard Green’s method, reads

$$\delta _\mathrm{b}(𝐤,a)=y^{\stackrel{~}{\alpha }}\left[A_1𝒥_{\stackrel{~}{\alpha }}(y)+A_2𝒥_{\stackrel{~}{\alpha }}(y)\right]+\left[24\stackrel{~}{\alpha }^2y^{5\stackrel{~}{\alpha }}S_{5\stackrel{~}{\alpha }1,\stackrel{~}{\alpha }}(y)\right]\delta _{\mathrm{DM}}(𝐤,a),$$
(22)
where $`S_{\mu ,\stackrel{~}{\alpha }}(y)`$ is the Lommel function.

By imposing the initial conditions (18) on our solution, we obtain $`A_1=B_1\delta _{\mathrm{DM}}(𝐤,a_i)`$ and $`A_2=B_2\delta _{\mathrm{DM}}(𝐤,a_i)`$, where

$`B_1`$ $`=`$ $`{\displaystyle \frac{\pi y_i^{1+\stackrel{~}{\alpha }}}{2\mathrm{s}\mathrm{i}\mathrm{n}(\pi \stackrel{~}{\alpha })}}\left[\left(124\stackrel{~}{\alpha }^2{\displaystyle \frac{S_{5\stackrel{~}{\alpha }1,\stackrel{~}{\alpha }}}{y_i^{5\stackrel{~}{\alpha }}}}\right)𝒥_{(\stackrel{~}{\alpha }+1)}{\displaystyle \frac{4\stackrel{~}{\alpha }}{y_i}}\left(112\stackrel{~}{\alpha }(2\stackrel{~}{\alpha }1)y_i^{15\stackrel{~}{\alpha }}S_{5\stackrel{~}{\alpha }2,\stackrel{~}{\alpha }+1}\right)𝒥_{\stackrel{~}{\alpha }}\right],`$
$`B_2`$ $`=`$ $`{\displaystyle \frac{\pi y_i^{1+\stackrel{~}{\alpha }}}{2\mathrm{s}\mathrm{i}\mathrm{n}(\pi \stackrel{~}{\alpha })}}\left[\left(124\stackrel{~}{\alpha }^2{\displaystyle \frac{S_{5\stackrel{~}{\alpha }1,\stackrel{~}{\alpha }}}{y_i^{5\stackrel{~}{\alpha }}}}\right)𝒥_{\stackrel{~}{\alpha }+1}+{\displaystyle \frac{4\stackrel{~}{\alpha }}{y_i}}\left(112\stackrel{~}{\alpha }(2\stackrel{~}{\alpha }1)y_i^{15\stackrel{~}{\alpha }}S_{5\stackrel{~}{\alpha }2,\stackrel{~}{\alpha }+1}\right)𝒥_{\stackrel{~}{\alpha }}\right],`$ (23)
and the Bessel and Lommel functions are evaluated at $`y=y_i=\sqrt{24}\stackrel{~}{\alpha }(k/\stackrel{~}{k}_J)a_i^{1/4\stackrel{~}{\alpha }}`$.

Once again, by replacing this expression in the continuity equation, we obtain the linear peculiar velocity of the baryons,

$$𝐮_\mathrm{b}(𝐤,a)=\left[\frac{y_i^{4\stackrel{~}{\alpha }}}{4\stackrel{~}{\alpha }}\left(B_1𝒥_{\stackrel{~}{\alpha }+1}B_2𝒥_{(\stackrel{~}{\alpha }+1)}\right)+12\stackrel{~}{\alpha }(2\stackrel{~}{\alpha }1)S_{5\stackrel{~}{\alpha }2,\stackrel{~}{\alpha }+1}\right]y^{15\stackrel{~}{\alpha }}𝐮_{\mathrm{DM}}(𝐤).$$
(24)
### 3.3 IGM linear filter

The results of the previous section indicate that one can always express the linear baryon density field as a linear convolution of the linear DM one. In Fourier space, one can therefore define an IGM linear filter, $`W_\mathrm{b}(k,a)`$, such that

$$\delta _\mathrm{b}(𝐤,a)=W_\mathrm{b}(k,a)\delta _{\mathrm{DM}}(𝐤,a).$$
(25)
For the case $`\alpha =1`$ we obtain

$$W_\mathrm{b}(k,a)=\left\{1+\left(\frac{k}{k_J}\right)^2\left(\frac{a}{a_i}\right)^{5/4}\left[\frac{5+Q}{2Q}\left(\frac{a}{a_i}\right)^{Q/4}\frac{5Q}{2Q}\left(\frac{a}{a_i}\right)^{Q/4}\right]\right\}\frac{1}{1+k^2/k_J^2},$$
(26)
which corresponds to eq. (A7) of (Gnedin & Hui 1998). Let us look at the main features of the IGM filter, at a finite time after reionization. It is immediate to check that on large scales $`W_\mathrm{b}`$ tends to unity: on scales much larger than the baryon Jeans length, the baryon distribution always traces the DM, due to the negligible role of pressure gradients. On scales much smaller than the Jeans length, instead, the baryon overdensity undergoes rapid oscillations. As noticed by previous authors (e.g. Gnedin & Hui 1998; Nusser 2000), and contrary to naive expectations, $`W_\mathrm{b}(k,a)`$ has no power-law asymptote on small scales.

In the general case $`\alpha 1`$, instead we have

$`W_\mathrm{b}(k,a)`$ $`=`$ $`24\stackrel{~}{\alpha }^2y^{5\stackrel{~}{\alpha }}S_{5\stackrel{~}{\alpha }1,\stackrel{~}{\alpha }}(y)`$ (27)
$``$ $`\left({\displaystyle \frac{y}{y_i}}\right)^{5\stackrel{~}{\alpha }}{\displaystyle \frac{\pi y_i}{2\mathrm{s}\mathrm{i}\mathrm{n}(\pi \stackrel{~}{\alpha })}}\left[\left(𝒥_{\stackrel{~}{\alpha }}(y)𝒥_{(\stackrel{~}{\alpha }+1)}(y_i)+𝒥_{\stackrel{~}{\alpha }}(y)𝒥_{\stackrel{~}{\alpha }+1}(y_i)\right)\left(124\stackrel{~}{\alpha }^2y_i^{5\stackrel{~}{\alpha }}S_{5\stackrel{~}{\alpha }1,\stackrel{~}{\alpha }}(y_i)\right)\right]`$
$``$ $`\left({\displaystyle \frac{y}{y_i}}\right)^{5\stackrel{~}{\alpha }}{\displaystyle \frac{2\pi \stackrel{~}{\alpha }}{\mathrm{sin}(\pi \stackrel{~}{\alpha })}}\left[\left(𝒥_{\stackrel{~}{\alpha }}(y)𝒥_{\stackrel{~}{\alpha }}(y_i)𝒥_{\stackrel{~}{\alpha }}(y)𝒥_{\stackrel{~}{\alpha }}(y_i)\right)\left(112\stackrel{~}{\alpha }(2\stackrel{~}{\alpha }1)y_i^{15\stackrel{~}{\alpha }}S_{5\stackrel{~}{\alpha }2,\stackrel{~}{\alpha }+1}(y_i)\right)\right].`$
It is worth mentioning that for half-integer values of $`\stackrel{~}{\alpha }`$ the Bessel and the Lommel functions in the above equation can be written in terms of trigonometric functions.

Analogous relations hold for the baryon peculiar velocity in terms of the DM one \[see eqs. (21) and (24) above\]. Also in this case, of course, $`W_\mathrm{b}1`$ in the large-scale limit, as it would be easy to see. The oscillating behavior of $`\delta _\mathrm{b}`$ on small scales is in general more complicated and depends upon the value of $`\alpha `$.

Plots of the IGM linear filter for some values of $`\alpha `$ are given in Figure 1, at various redshifts. Hydrogen reionization is assumed to have occurred at a redshift $`z_i=7`$. Note that, on scales smaller than the baryon Jeans length, acoustic oscillations, with decaying amplitude, persist even long after reionization. Note also that different equations of state lead to a different small-scale behaviour even if the IGM density is plotted vs. the scaled wavenumber $`\kappa k/\stackrel{~}{k}_J`$. This is because eq. (16), which governs the linear evolution of baryon density fluctuations, retains a dependence on the slope $`\alpha `$ in the last term on the LHS, which describes the effect of pressure gradients.

## 4 Modelling the weakly nonlinear IGM dynamics

With our choice of time variable, the Euler equation, both for the DM and the baryons, takes the simple form

$$\frac{D𝐮}{Da}=\eta ,$$
(28)
where the potential $`\eta `$ takes a different expression for the two fluids. As the force acting on both fluids is conservative, in the absence of initial vorticity the flow remains irrotational (actually, this is only true prior to shell-crossing, for the DM component), and we can express the peculiar velocity in terms of a velocity potential, $`𝐮\mathrm{\Phi }`$. The potential $`\eta `$ in the RHS of the Euler equation reads

$$\eta _{\mathrm{DM}}=\frac{3}{2a}(\mathrm{\Phi }_{\mathrm{DM}}+\phi );\eta _\mathrm{b}=\frac{3}{2a}[\mathrm{\Phi }_\mathrm{b}+\phi +\frac{1}{(\gamma 1)ak_J^2}(1+\delta _\mathrm{b})^{\gamma 1}],$$
(29)
for the DM and IGM components respectively.

The dynamical model we introduce here is largely inspired by the Zel’dovich approximation (Zel’dovich 1970), and is based on replacing the potential $`\eta `$, which should be consistently calculated using the full set of equations for the two fluids, by a mock external potential, obtained by evaluating its expression within linear (Eulerian) perturbation theory. Indeed, this model might be seen as only the first step of a more general, iterative approximation scheme.

In order to evaluate to linear order the quantity $`\eta `$ one can either compute to first order any single contribution or compute its Laplacian by taking the divergence of the LHS of the linearized Euler equation. This yields $`^2\eta \dot{𝐮}`$, and, using the linearized continuity equation, $`^2\eta \ddot{\delta }.`$

We thus conclude that, to first order in perturbation theory the potential $`^2\eta `$ is given by the second ‘time’-derivative of the linear overdensity. This fact immediately implies that, to first order, $`\eta _{\mathrm{DM}}=0`$ (neglecting the contribution from decaying modes), while

$$\eta _\mathrm{b}(k,a)=\left[a\ddot{W}_\mathrm{b}(k,a)+2\dot{W}_\mathrm{b}(k,a)\right]\phi _0(k)(k,a)\phi _0(k).$$
(30)
In particular, for $`\alpha =1`$ only the homogeneous part of the $`\delta _\mathrm{b}`$ solution contributes, and we obtain

$$(k,a)=\left(\frac{a}{a_i}\right)^{5/4}\left[\frac{(Q1)(Q5)}{16a}\left(\frac{a}{a_i}\right)^{Q/4}\stackrel{~}{B}_1+\frac{(Q+1)(Q+5)}{16a}\left(\frac{a}{a_i}\right)^{Q/4}\stackrel{~}{B}_2\right],$$
(31)
with $`Q=\sqrt{124k^2/k_J^2}`$ and the integration constants $`\stackrel{~}{B}_1`$ and $`\stackrel{~}{B}_2`$ as in eq. (20). For $`\alpha 1`$ we get

$`(k,a)`$ $`=`$ $`y_i^{4\stackrel{~}{\alpha }}\left(\sqrt{24}\stackrel{~}{\alpha }{\displaystyle \frac{k}{\stackrel{~}{k}_J}}\right)^{4\stackrel{~}{\alpha }}{\displaystyle \frac{y^{19\stackrel{~}{\alpha }}}{16\stackrel{~}{\alpha }^2}}\left(B_1\left[(4\stackrel{~}{\alpha }2)𝒥_{\stackrel{~}{\alpha }+1}(y)+y𝒥_{\stackrel{~}{\alpha }+2}(y)\right]+B_2\left[(24\stackrel{~}{\alpha })𝒥_{(\stackrel{~}{\alpha }+1)}(y)+y𝒥_{(\stackrel{~}{\alpha }+2)}(y)\right]\right)`$ (32)
$`+`$ $`\left(\sqrt{24}\stackrel{~}{\alpha }{\displaystyle \frac{k}{\stackrel{~}{k}_J}}\right)^{4\stackrel{~}{\alpha }}y^{19\stackrel{~}{\alpha }}(12\stackrel{~}{\alpha }6)\left[2(\stackrel{~}{\alpha }1)yS_{5\stackrel{~}{\alpha }3,\stackrel{~}{\alpha }+2}+(12\stackrel{~}{\alpha })S_{5\stackrel{~}{\alpha }2,\stackrel{~}{\alpha }+1}\right],`$
where $`B_1`$ and $`B_2`$ are given by eq. (23). Plots of the function $`(k,a)`$ are given in Figure 2, at different redshifts after reionization, for various values of $`\alpha `$.

Therefore, our model is described by the two Euler equations:

$$\frac{D𝐮_{\mathrm{DM}}}{Da}=0,$$
(33)
for the DM component, and

$$\frac{D𝐮_\mathrm{b}}{Da}=\eta _\mathrm{b},$$
(34)
for the baryons.

The solution of the above DM equation of motion is well known: it corresponds to the Zel’dovich approximation, according to which mass elements move along straight lines, with constant ‘velocity’ impressed by local linear fluctuations of the gravitational force at their initial Lagrangian location $`𝐪`$: $`𝐮_{\mathrm{DM}}(𝐱(𝐪,a),a)=_𝐪\phi _0(𝐪)`$. We therefore have

$$𝐱_{\mathrm{DM}}(𝐪,a)=𝐪a_𝐪\phi _0(𝐪).$$
(35)
The trajectories of the IGM fluid elements are instead more complex, as our equation implies a non-zero force acting on them,

$$𝐱_\mathrm{b}(𝐪,a)=𝐪a_𝐪\phi _0(𝐪)_0^a𝑑\tau _0^\tau 𝑑\tau ^{}_𝐱\eta (𝐱(𝐪,\tau ^{}),\tau ^{}),$$
(36)
where we have formally extended the time integration from $`0`$ to $`a`$, as, before reionization, when the Jeans length is negligible, $`\eta _\mathrm{b}(𝐱,a<a_i)=0`$.

It will also prove convenient in the following to transform the baryon Euler equation (34) into the following Bernoulli, or Hamilton-Jacobi, equation for the velocity potential:

$$\frac{\mathrm{\Phi }_\mathrm{b}}{a}+\frac{1}{2}(\mathrm{\Phi }_\mathrm{b})^2=\eta _\mathrm{b}.$$
(37)
The extreme simplicity of our scheme is shown by the fact that the only information needed to evolve the baryon distribution in the weakly nonlinear regime is the IGM filter $`W_\mathrm{b}`$. It is immediate to realize that, because of our derivation of $`\eta _\mathrm{b}`$, its expression, $`^2\eta \ddot{\delta }`$, in terms of the linear baryon overdensity is fully general, i.e. it would apply to general reionization histories \[i.e. to general $`T_0(z)`$ relations\], general baryon equations of state and to the case in which the baryon self-gravity is properly taken into account (i.e. the general case $`0<f_\mathrm{b}1`$). Note that the IGM linear filter is simply the ratio of the baryon to DM transfer function at the given time $`W_\mathrm{b}(k,a)=T_\mathrm{b}(k,a)/T_{\mathrm{DM}}(k,a)`$ (provided the reionization process has been taken into account in evaluating the baryon transfer function $`T_\mathrm{b}`$).

One might think that there is some degree of arbitrariness in the particular choice we made of which terms to linearize \[those on the RHS of eq. (28)\] and which ones to treat exactly (those on its LHS). The strongest motivation for such a choice is its analogy with the procedure leading to the Zel’dovich approximation for the DM component. In this sense, the choice is unique; as we will see in the next subsection, in fact, the external force $`\eta _\mathrm{b}`$ obtained above and its effect on the IGM trajectories are closely connected to the results of first-order Lagrangian perturbation theory: no other choice would have provided such a connection. It should also be stressed that other successful approximation schemes for the evolution of the DM component, such as the ‘frozen flow’ (Matarrese et al. 1992) and ‘frozen potential’ (Brainerd, Scherrer & Villumsen 1993; Bagla & Padmanabhan 1994) ones, are indeed based on the same choice. The physical reason which makes a linear theory evaluation of $`\eta `$ reasonably accurate is that this quantity, for both fluids, contains terms that receive their dominant contribution from small wavenumbers. This was indeed the original motivation which led Zel’dovich to obtain his celebrated algorithm, although in modern language the Zel’dovich approximation is more commonly explained within the first-order Lagrangian approximation (e.g. Sahni & Coles 1995 and references therein). The link between these two alternative derivations is discussed in the next subsection.

### 4.1 Lagrangian approximation to the baryon trajectories

The trajectory of any mass element can be written in the general form

$$𝐱(𝐪,a)=𝐪+𝐒(𝐪,a),$$
(38)
where $`𝐒`$ is the ‘displacement vector’. As far as the evolution of the fluid is far from the strongly nonlinear regime, the displacement vector can be considered to be small, i.e. the Eulerian and Lagrangian positions of each particle are never too far apart. This consideration is at the basis of Lagrangian approximation methods. The basics of the first-order Lagrangian scheme applied to our two-component fluid are reported in Appendix A. Applying these ideas to our external force, we obtain

$$_𝐱\eta _\mathrm{b}(𝐱(𝐪,a),a)=_𝐪\eta _\mathrm{b}(𝐪,a)+𝒪\left(𝐒^2\right).$$
(39)
To first order in the displacement vector we can therefore replace the Eulerian force $`_𝐱\eta _\mathrm{b}(𝐱,a)`$ by its Lagrangian counterpart $`_𝐪\eta _\mathrm{b}(𝐪,a)`$ in the baryon trajectories, which leads to the much simpler form

$$𝐱_\mathrm{b}(𝐪,a)=𝐪a_𝐪\psi _\mathrm{b}(𝐪,a),$$
(40)
where the ‘baryon potential’ $`\psi _\mathrm{b}`$ is defined by

$$_𝐪^2\psi _\mathrm{b}(𝐪,a)=\frac{\delta _\mathrm{b}(𝐪,a)}{a}$$
(41)
and is related to the peculiar gravitational potential by the Fourier-space expression $`\psi _\mathrm{b}(𝐤,a)=W_\mathrm{b}(k,a)\phi _0(𝐤)`$.

Moreover, as shown in Appendix A, the trajectories described by eq. (40) represent the result of first-order Lagrangian perturbation theory applied to our baryon gas. Similar results have been recently obtained by Adler & Buchert (1999) for the case of a single self-gravitating collisional fluid (i.e. for $`f_\mathrm{b}=1`$). There is an important difference between these trajectories and the DM ones. According to the Zel’dovich approximation, DM particles move along straight lines with constant velocity, whereas the baryons are generally accelerated along curved paths; this is due to the non-zero force, resulting from the composition of three terms: the local Hubble drag, the local gravitational force caused by the dominant DM component and the gradient of the gas pressure.

The baryon peculiar velocity which follows from this first-order Lagrangian approximation is

$$𝐮_\mathrm{b}(𝐱(𝐪,a),a)=_𝐪\left(a\dot{\psi }_\mathrm{b}(𝐪,a)+\psi _\mathrm{b}(𝐪,a)\right),$$
(42)
which, unlike the DM one, deviates from the initial velocity $`𝐮_{\mathrm{b0}}(𝐪)𝐮_\mathrm{b}(𝐱(𝐪,0),0)=_𝐪\phi _0(𝐪)`$.

Let us mention an important property of the approximate expression for the velocity field that we obtained using first-order Lagrangian perturbation theory: unlike what happens in the DM case, the vector $`𝐮_\mathrm{b}`$ of eq. (42) is irrotational in Eulerian space only if we restrict its validity to first-order in the displacement vector; a non-zero velocity curl component, $`_𝐱\times 𝐮_\mathrm{b}0`$, in fact arises beyond this limit. This feature is not expected to imply serious problems as long as the system has not evolved too deeply into the nonlinear regime, i.e. until the baryon overdensity has not reached values $`1`$. This is an unphysical feature deriving from the extrapolation of the first-order results to a regime where higher-order terms should be taken into account. A similar feature appears in the Lagrangian perturbation approach to the DM dynamics, when both the growing and decaying solutions are included (e.g. Buchert 1992). No problems of this type, however, occur when the Eulerian scheme above is adopted, i.e. when eq. (34) and its solution (36) are assumed. For this reason we will consider our ‘Eulerian model’ of eqs. (34) and (36) as the correct one and the Lagrangian scheme leading to eq. (40) as being essentially a shortcut to get approximated baryon trajectories.

The slight discrepancy between the Eulerian and Lagrangian schemes used to derive the external force $`\eta `$ acting on the baryons is a peculiarity of the collisional case. In the DM case, these two techniques – linearizing the RHS of eq. (28) in Eulerian space, and expanding to first order in Lagrangian perturbation theory – lead to identical results.

Approximation schemes for the low-density IGM dynamics which are closely related to our Lagrangian treatment have been studied by Reisenegger & Miralda-Escudé (1995), Gnedin & Hui (1996), and Hui, Gnedin & Zhang (1997). There are, however, important differences that we would like to point out: i) in our model the gas trajectories, and thus the resulting IGM spatial clustering depend on the specific ionization history (different values of $`\alpha `$ in the simplest case); ii) our model is by definition able to exactly reproduce the behavior of the baryon component in the linear regime. In practice, while previous models adopt an IGM filter which is best-fitted to simulations, ours directly derives from baryon dynamics.

#### 4.1.1 Numerical simulations of the IGM distribution in the laminar regime

In order to display the effect of the IGM filter on baryon trajectories, we have produced a set of numerical simulations based on the Lagrangian scheme. $`256^3`$ particles were laid down on a uniform cubic grid and moved according to eq. (40). Realizations of the peculiar gravitational potential were obtained in Fourier-space according to standard algorithms. The model shown in Figure 3 is a spatially flat, vacuum dominated cold dark matter (CDM) model, with $`h=0.65`$ (where $`h`$ is the Hubble constant, $`H_0`$, in units of $`100`$ km s<sup>-1</sup> Mpc<sup>-1</sup>), a cosmological constant with present-day density parameter $`\mathrm{\Omega }_{0\mathrm{\Lambda }}=0.7`$, CDM with $`\mathrm{\Omega }_{0\mathrm{C}\mathrm{D}\mathrm{M}}=0.26`$ and baryons with the remaining $`\mathrm{\Omega }_{0\mathrm{b}}=0.04`$. A Zel’dovich primordial power-spectrum of adiabatic perturbations is assumed, with the standard normalization of $`\mathrm{\Lambda }`$CDM models<sup>2</sup><sup>2</sup>2To correct for the small inaccuracy introduced by the use of $`a`$, instead of the DM growth-factor, as time variable in our treatment, we apply a correction factor $`(1+z)D_+(z)`$ to the $`\sigma _8`$ normalization, where $`D_+(z)`$ is the linear growth-factor of DM density fluctuations, normalized to unity at $`z=0`$., $`\sigma _8=0.99`$ (Viana & Liddle 1999), where $`\sigma _8`$ is the rms mass fluctuation in a sharp-edged sphere of radius $`8h^1`$ Mpc; we adopt the Bardeen et al. (1986) CDM transfer function, corrected to account for the baryon contribution, as in (Sugiyama 1995).

The box-size is $`10.24`$ Mpc. Three different IGM thermal histories are considered: in all cases we assume that reionization occurred at $`z_i=7`$, but we assume that the mean IGM temperature evolves according to different power-laws afterwards, $`T_0(z)(1+z)^\alpha `$, with $`\alpha =1.0,0.8,0.6`$. Figure 3 shows a thin ($`0.02`$ Mpc) slice of our simulation box at $`z=3`$. The corresponding DM distribution is also shown for comparison. Note that different values of $`\alpha `$ lead to a different small-scale distribution of the IGM, because of i) the different time dependence of the Jeans scale and ii) the different trajectories followed by the baryons after reionization. The comoving Jeans wavenumber $`k_J`$ has been set to the same value ($`7`$ Mpc<sup>-1</sup>) in all our models, at redshift $`z=3`$.

#### 4.1.2 The IGM density field in the laminar regime

Once particle trajectories are known, it is immediate to obtain the density field, using mass conservation, $`1+\delta _\mathrm{b}(𝐱,a)=𝐱_\mathrm{b}(𝐪,a)/𝐪^1`$. If we adopt the Lagrangian expression of eq. (40), we get

$$1+\delta _\mathrm{b}(𝐱(𝐪,a),a)=\left(\mathrm{det}\left[\delta _{ij}a\frac{^2\psi _\mathrm{b}(𝐪,a)}{q_iq_j}\right]\right)^1.$$
(43)
The strain tensor $`^2\psi _\mathrm{b}/q_iq_j`$ can be locally diagonalized along principal axes, whose direction will generally depend upon time, unlike the DM case. Calling $`\lambda _\mathrm{i}`$ the corresponding eigenvalues, we write

$$1+\delta (𝐱(𝐪,a),a)=\underset{i=1}{\overset{3}{}}\left[1a\lambda _i(𝐪,a)\right]^1,$$
(44)
which shows that a caustic singularity will form at the finite time $`a_{\mathrm{sc}}(𝐪)=1/\lambda _{}(𝐪,a_{\mathrm{sc}})`$, where $`\lambda _{}(𝐪,a_{\mathrm{sc}})`$ is the largest eigenvalue of the time-dependent strain tensor. The time dependence of the eigenvalues, which on small scales becomes oscillatory, implies that the shell-crossing condition can even be met more than once by a given mass element along the same principal axis.

The extrapolation of our Lagrangian approximation beyond the formation of the first pancakes, leads to an artificial diffusion of these structures. This problem becomes more and more severe at low redshift, making any simplified description of the IGM – and of the Ly$`\alpha `$ forest – based on the Lagrangian trajectories quite unreliable. The stochastic adhesion model presented in the next section aims precisely at overcoming this problem.

One might use eq. (44) for the density field to obtain an analytical expression for the probability distribution function (PDF) of the IGM density field, by a simple extension of the traditional Doroshkevich (1970) formalism (e.g. Kofman et al. 1994). Such a technique has been applied by Reisenegger & Miralda-Escudé (1995), who adopted a smoothed Zel’dovich approximation for the IGM evolution. Our model, however, allows to get a more refined description in which the PDF is sensitive to the IGM thermal history. Similarly, one might also obtain statistical information on the Ly$`\alpha `$ forest in quasar absorption spectra, e.g. in terms of the column density distribution function of the lines (e.g. Hui, Gnedin & Zhang 1997); this too would display a dependence on the assumed $`T_0(z)`$, which might be tested against observational data. These, and related applications of our Lagrangian algorithm will be presented elsewhere.

The emergence of shell-crossing singularities in the velocity pattern of our collisional fluid should be understood as an artifact of the linearized treatment of pressure gradients in the Euler equation. This feature can be also seen as follows. If we evaluate the force $`\eta `$ on large scales, to lowest order in $`k/k_J`$ we find, for $`\alpha =1`$,

$$(k,a)\frac{3}{2a}\left(\frac{a}{a_i}\right)^{1/2}\left(\frac{k}{k_J}\right)^2,$$
(45)
which comes only from the homogeneous part of eq. (19), and

$$(k,a)\frac{6}{a}\frac{(1\alpha )}{(2\alpha 1)(52\alpha )}\left(\frac{k}{k_J}\right)^2,$$
(46)
for $`\alpha 1`$, from the inhomogeneous term of eq. (22). In both cases the result can be expressed in the form $`\eta \nu _{\mathrm{eff}}(a)^2𝐮_0`$, where $`\nu _{\mathrm{eff}}>0`$, for $`0.5<\alpha 1`$. This shows that, as long as the velocity is in the linear regime, the effect of the linearized pressure gradient is similar to that of an effective kinematical viscosity<sup>3</sup><sup>3</sup>3A related point has been recently made by Buchert & Domínguez (1998)., with coefficient $`\nu _{\mathrm{eff}}`$, which smooths the velocity field. As soon as the system enters the mildly non-linear regime, however, $`𝐮`$ deviates from its initial value and this smoothing is no more effective in stopping the formation of caustic singularities in the densest regions.

Two processes will actually prevent the formation of multi-streams in our physical system: the first is due to the binding action of gravity, which tends to keep pancakes and filaments thin, an effect which is experienced both by the DM and IGM components; the second is the nonlinear action of the gas pressure, which is instead characteristic of the baryonic component. Because of the difficulty to deal analytically with the gas pressure beyond any perturbative approximation, we will find convenient in the next section to introduce an artificial viscosity term in the baryon Euler equation, whose effect is to smear these shell-crossing singularities of the velocity field.

## 5 The stochastic adhesion approximation

The model for the formation of structures in the baryon distribution discussed in the last section is based on a suitable approximation for the force exerted on fluid elements by gravity and by the surrounding fluid patches through pressure gradients. The latter are most important on small scales where their effect is to smooth the gas fluctuations relative to the DM ones, thus preventing the occurrence of shell-crossing singularities in the collisional component; nonlinear effects here manifest themselves through the formation of shock-waves. In our scheme, just like in all similar schemes \[e.g. the modified Zel’dovich scheme used by Gnedin & Hui (1996) and Hui, Gnedin & Zhang (1997)\], the gas pressure is only included through a linear approximation. Because of this fact, when the system enters the nonlinear regime, i.e. when two particles come very close in space, shell-crossing can affect their evolution leading to the subsequent occurrence of multi-streams. In order to prevent this unphysical phenomenon we will adopt the same technique which proved so successful in the DM case to extend the Zel’dovich treatment beyond its actual range of validity: we add a kinematical viscosity term to our approximate Euler equation. This method is at the basis of the adhesion approximation for the formation of large-scale structure in the collisionless case (Gurbatov et al. 1985, 1989; Kofman & Shandarin 1988). The adhesion model is based on the three-dimensional generalization of Burgers’ equation of strong turbulence. It is the simplest equation which describes the formation and subsequent merging of shocks. According to the adhesion model, DM particles move according to the Zel’dovich approximation until they fall into pancakes, when, owing to the viscous force, they stick together. Next, pancakes drain into filaments, and filaments into clumps. The thickness of these structures is monitored by the value of the kinematical viscosity coefficient $`\nu `$, and becomes infinitely thin as $`\nu 0`$.

We assume that the equation of motion which governs the IGM dynamics is (from now on we will avoid the subscript ‘b’ on baryon quantities, where unnecessary)

$$\frac{D𝐮}{Da}=\nu ^2𝐮\eta ,$$
(47)
where the kinematical viscosity coefficient $`\nu `$ is here assumed to be small, but non-zero. In our collisional case, moreover, there can be an extra reason to add such a term: the gas is effectively experiencing some shear viscosity, although on scales much smaller than those under consideration. The physical viscosity term should actually depend upon space and time, through the inverse of the local gas density. One might however argue that, in the physically relevant limit, where only a tiny viscous force is present, even a constant $`\nu `$ will produce the correct qualitative trend. It might be interesting to mention an alternative interpretation of the viscosity term in self-gravitating systems. According to Domínguez (2000), a term which resembles the one for kinematical viscosity is the unavoidable consequence of the coarse-graining process inherent in a hydrodynamical description.

The equation above is known as the forced Burgers equation of nonlinear diffusion (e.g. Barabási & Stanley 1995; Frisch & Bec 2000). Its possible application to the cosmological structure formation problem was suggested long ago by Zel’dovich et al. (1985, 1987) and recently studied in greater detail by Jones (1996, 1999). There are a number of differences between our approach and the one by Jones, which lead to a different dynamics of the IGM. First, we used a different time variable (the scale-factor $`a`$ instead of the cosmic time $`t`$) to define the baryon velocity field and its acceleration, thus leading to the form of eq. (28), where the force on the RHS was replaced by its linearized expression. Second, we do include some effect of the baryonic pressure in our external random potential $`\eta `$, which is instead absent in Jones’ treatment, where the role of $`\eta `$ is played by the linear gravitational potential generated by the dark matter component. Thus, our scheme, unlike the one by Jones, is able to exactly reproduce the evolution of the baryons at the linear level. Finally, our random potential has a non-trivial time dependence, unlike the one assumed by Jones.

The forced Burgers equation can be transformed into a Bernoulli-like equation for the velocity potential, namely

$$\frac{\mathrm{\Phi }}{a}+\frac{1}{2}(\mathrm{\Phi })^2=\nu ^2\mathrm{\Phi }\eta .$$
(48)
The problem is fully specified once the initial conditions for the velocity potential and the statistics of the noise are given. Our initial velocity potential is $`\mathrm{\Phi }_0(𝐪)=\phi _0(𝐪)`$. The statistics of the noise follows directly from that of the linear gravitational potential, which we assume to be a Gaussian random field. Then, also $`\eta `$ is a Gaussian process with zero mean and auto-correlation function

$$\eta (𝐱,a)\eta (𝐱+𝐫,a^{})=\frac{1}{2\pi ^2}_0^{\mathrm{}}\frac{dk}{k^2}(k,a)(k,a^{})P(k)j_0(kr),$$
(49)
where the average is over the $`\phi _0`$ ensemble, $`P(k)`$ is the linear power-spectrum of DM density fluctuations, extrapolated to the present time, and $`j_0`$ is the spherical Bessel function of order zero. According to this formula, $`\eta `$ is a stochastic process with the following properties:

i) It is a colored (i.e. non-white) Gaussian random process both in space and time.

ii) It is stationary (i.e. homogeneous and isotropic) in space but not in time (as its auto-correlation function does depend on $`|𝐱𝐱^{}|`$ only, but not on $`|aa^{}|`$: this is typical of the cosmological case.

iii) On small scales, $`\eta `$ oscillates rapidly in time for all $`\alpha 1`$, with a period which generally depends on the wavenumber.

iv) The non-separability of the space and time dependence of its auto-correlation function implies that $`\eta `$ behaves as a stochastic process both in space and time. More precisely, the time evolution of the noise in each given spatial point cannot be predicted on the basis of local initial conditions only, as it depends on the realization of the Gaussian field $`\phi _0`$ over the whole Lagrangian space, through the time evolution of the IGM filter.

Let us finally stress an important peculiarity of our cosmological application of the forced Burgers equation: the same stochastic process $`\phi _0`$ which underlies the external random potential $`\eta =\phi _0`$ also provides the initial condition for the velocity potential, $`\mathrm{\Phi }_0=\phi _0`$.

A plot of the function $`(k,a)`$ vs. redshift is given in Figure 4 for various values of $`\alpha `$, to display the time-dependence of the random force and, in particular, its oscillatory character on small scales. In Figure 5 the dimensionless power-spectra $`\mathrm{\Delta }^2(k)`$ of two terms contributing to the evolution of the velocity potential are shown: that of the random potential $`\eta `$, $`\mathrm{\Delta }_\eta ^2(k,z)(1/2\pi ^2)k^1^2(k,z)P(k)`$, and that of the linear velocity potential $`\mathrm{\Phi }_0`$, multiplied by a factor $`1/a=1+z`$ to give an estimate of the ‘deterministic’ term in eq. (48), namely $`\mathrm{\Delta }_{(1+z)\mathrm{\Phi }_0}^2(k,z)(1/2\pi ^2)(1+z)^2k^1P(k)`$. Note that the two contributions become of the same order at $`kk_J`$ at the relevant redshifts; on smaller scales the random potential dominates the IGM dynamics.

The above evolution equation for $`\mathrm{\Phi }`$, in cases when the random potential $`\eta `$ is white-noise in time, has been extensively studied in condensed matter physics, where it became popular as the Kardar-Parisi-Zhang (KPZ) equation (Kardar, Parisi & Zhang 1986). The KPZ equation describes the growth of an interface under random particle deposition. Here the viscosity coefficient $`\nu `$ plays the role of temperature and the velocity potential $`\mathrm{\Phi }`$ is interpreted as the height above an initially flat surface, which is driven by the random noise and gradually becomes rough.

### 5.1 Solving the random heat equation

Equation (48) can be easily related to a linear, partial-differential equation, by means of the nonlinear Hopf-Cole transformation $`\mathrm{\Phi }=2\nu \mathrm{ln}𝒰`$ (e.g. Burgers 1974), which leads to the linear diffusion equation with a multiplicative potential, also called ‘random heat’ equation

$$\frac{𝒰(𝐱,a)}{a}=\nu ^2𝒰(𝐱,a)+\frac{\eta (𝐱,a)}{2\nu }𝒰(𝐱,a).$$
(50)
Starting from the solution of the latter equation for the Hopf-Cole transformed velocity potential $`𝒰`$ \[which has been dubbed ‘expotential’ by Weinberg & Gunn (1990a)\] one can easily find the velocity by transforming everything back

$$𝐮=2\nu \frac{𝒰}{𝒰}$$
(51)
One first obtains an expression for the ‘transition kernel’ $`𝒦(𝐱,a|𝐪,0)`$, representing the particular solution obtained from a Dirac delta function, $`\delta (𝐱𝐪)`$, at the initial time $`a=0`$. This has a formal solution given in terms of the (Euclidean) Feynman-Kac path-integral <sup>4</sup><sup>4</sup>4The solution of eq. (50) in the absence of a potential, i.e. for the ‘free’ adhesion approximation, is reviewed in Appendix B. formula (e.g. Feynman & Hibbs 1965), namely

$$𝒦(𝐱,a|𝐪,0)=_𝐪^{𝐱(a)}\left[𝒟𝐱(\tau )\right]e^{S/2\nu },$$
(52)
where the action is given by

$$S[𝐱]=_0^a𝑑\tau (𝐱,\dot{𝐱},\tau )=_0^a𝑑\tau \left[\frac{\dot{𝐱}^2}{2}\eta (𝐱,\tau )\right],$$
(53)
with $``$ the Lagrangian of a particle moving in the potential $`\eta `$.

In our diffusion equation the transition kernel is immediately understood as the conditional probability of finding a particle in $`𝐱`$, at time $`a`$, given that it was initially in the Lagrangian position $`𝐪`$. To better understand the path-integral solution it is however convenient to think it in connection with a quantum mechanical problem. It is in fact immediate to realize that an inverse ‘Wick rotation’ of the time variable $`aia`$, together with the formal replacement $`2\nu \mathrm{}`$ transforms eq. (50) into the Schrödinger equation for a particle, of unit mass, subjected to the potential $`\eta `$. According to the path-integral representation of quantum mechanics, its solution is obtained by ‘integrating’ over all possible paths connecting these two end-points, each one weighed by the action $`S[𝐱]`$, calculated along this path.

Once the kernel is known, the solution of the random heat equation with appropriate initial conditions is obtained through the application of the Chapman-Kolmogorov equation (e.g. van Kampen 1992), by integrating the product of the initial function $`𝒰_0`$ with the transition kernel over the whole Lagrangian space, namely $`𝒰(𝐱,a)=d^3q𝒰_0(𝐪)𝒦(𝐱,a|𝐪,0)`$. In our case $`𝒰_0(𝐪)=\mathrm{exp}[\mathrm{\Phi }_0(𝐪)/2\nu ]=\mathrm{exp}[\phi _0(𝐪)/2\nu ]`$. Thus,

$$𝒰(𝐱,a)=d^3qe^{\mathrm{\Phi }_0(𝐪)/2\nu }_𝐪^{𝐱(a)}\left[𝒟𝐱(\tau )\right]e^{S/2\nu }.$$
(54)
In the limit of vanishing viscosity (corresponding to the classical limit, $`\mathrm{}0`$, in our quantum mechanical analog) the dominant contribution to the path-integral comes from the ‘classical’ path, i.e. that which satisfies the Euler-Lagrange equations of motion,

$$\frac{}{𝐱}\frac{d}{da}\frac{}{\dot{𝐱}}=0,$$
(55)
which in our case leads to Newton’s second law, $`\ddot{𝐱}=\eta `$. Thus, the particle trajectories along the classical path read

$$𝐱_{\mathrm{cl}}(𝐪,a)=𝐪+a𝐮_0(𝐪)_0^a𝑑\tau _0^\tau 𝑑\tau ^{}\eta (𝐱(𝐪,\tau ^{}),\tau ^{}),$$
(56)
with general initial velocity $`𝐮_0`$. Note that there are still infinitely many classical trajectories joining the two end-points $`(𝐪,0)`$ and $`(𝐱,a)`$, corresponding to the freedom to choose the initial velocity $`𝐮_0(𝐪)`$.

We can now expand around the classical trajectories in the standard manner of quadratic approximations (e.g. Feynman & Hibbs 1965), that is $`𝐱=𝐱_{\mathrm{cl}}+\xi `$, subject to the constraint that the fluctuations around the classical trajectory vanish at the end points, namely $`\xi (0)=\xi (a)=0`$. We then have

$$S[𝐱]=S_{\mathrm{cl}}+𝑑\tau \xi (\tau )\left(\frac{\delta S[𝐱]}{\delta 𝐱}\right)_{𝐱_{\mathrm{cl}}}+\frac{1}{2}\underset{i,j=1}{\overset{3}{}}𝑑\tau _1𝑑\tau _2\xi _i(\tau _1)\xi _j(\tau _2)\left(\frac{\delta ^2S[𝐱]}{\delta x_i(\tau _1)\delta x_j(\tau _2)}\right)_{𝐱_{\mathrm{cl}}},$$
(57)
with the symbol $`\delta `$ standing for functional differentiation. The classical action, $`S_{\mathrm{cl}}`$, is a function of the end-points $`(𝐪,𝐱)`$ and of the time interval $`a`$, which has to be calculated by replacing the solution of the Euler-Lagrange equations into its expression.

The action is an extremum for the classical trajectory, thus $`(\delta S/\delta 𝐱)_{𝐱_{\mathrm{cl}}}=0`$. Inserting the expansion (57) in the path-integral gives

$$𝒦(𝐱,a|𝐪,0)=e^{S_{\mathrm{cl}}/2\nu }_\mathrm{𝟎}^\mathrm{𝟎}\left[𝒟\xi (\tau )\right]e^{\frac{1}{4\nu }_{i,j=1}^3_0^a𝑑\tau _0^a𝑑\tau _1_0^a𝑑\tau _2\xi _i(\tau _1)\xi _j(\tau _2)\left(\frac{\delta ^2S}{\delta x_i(\tau _1)\delta x_j(\tau _2)}\right)_{𝐱=𝐱_{\mathrm{cl}}}},$$
(58)
where we used the fact that the Jacobian of the transformation from $`𝐱`$ to $`\xi `$ is unity. Note that the path-integral in (58) starts and ends at zero, because the fluctuations vanish at the end-points. Contrary to an ordinary integral, the fact that the upper and lower limits are the same does not imply that the path-integral vanishes. This can be clearly seen by considering the path-integral of a free particle (see Appendix B), which can be evaluated exactly and does not vanish even if one is only interested in the path which starts and ends at the same point (e.g. Feynman & Hibbs 1965). Since the integral over $`\xi `$ starts and ends at $`\mathrm{𝟎}`$, it can only be a function, $`F(a)`$, of the end-times. Therefore, we can rewrite the above path-integral as $`𝒦(𝐱,a|𝐪,0)=F(a)\mathrm{exp}[S_{\mathrm{cl}}/2\nu ]`$. We need not find the explicit expression for the pre-factor $`F(a)`$, since this function cancels for the velocity. The solution of the random heat equation, for small $`\nu `$, reads

$$𝒰(𝐱,a)=F(a)d^3qe^{\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪,a)/2\nu },$$
(59)
with

$$\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪,a)S_{\mathrm{cl}}(𝐱,𝐪,a)+\mathrm{\Phi }_0(𝐪).$$
(60)
Replacement into eq. (51) yields

$$𝐮(𝐱,a)=\frac{d^3q_𝐱S_{\mathrm{cl}}(𝐱,𝐪,a)e^{\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪,a)/2\nu }}{d^3qe^{\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪,a)/2\nu }}.$$
(61)
In the case of vanishing force, eq. (61) reduces to the well-known exact solution of the three-dimensional Burgers equation (see Appendix B), which underlies the adhesion approximation for the dark matter component. In our case, it represents a useful approximation, which could be used for a numerical evaluation of the baryon velocity field, by extending the technique applied by Nusser & Dekel (1990) and Weinberg & Gunn (1990a,b). The main practical difficulty in using the solution (61) is, however, represented by the lack of an explicit expression for the force $`\eta `$, as a function of $`𝐱`$ and $`a`$. A possible shortcut would be to exploit the Lagrangian approximation for the trajectories, introduced in Section 4.1. This method will be discussed in Section 6.

Once the peculiar velocity field of eq. (61) is known at each Eulerian point $`𝐱`$ as a function of time, one can obtain the actual baryon trajectories, by numerically integrating the integral equation (e.g. Nusser & Dekel 1990; Weinberg & Gunn 1990a,b)

$$𝐱(𝐪,a)=𝐪+_0^a𝑑\tau 𝐮(𝐱(𝐪,\tau ),\tau ),$$
(62)
starting from Lagrangian particle positions on a grid. From these trajectories the baryon density can be found, according to the standard formula $`1+\delta (𝐱,a)=𝐱(𝐪,a)/𝐪^1`$.

### 5.2 Steepest-descent approximation

In the limit of small $`\nu `$ we can evaluate the integral over initial positions in eq. (59) using the steepest-descent, or saddle-point, approximation. We have

$`𝒰(𝐱,a)`$ $`=`$ $`F(a)e^{\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪_s,a)/2\nu }{\displaystyle d^3\delta qe^{\frac{1}{4\nu }_{i,j=1}^3\delta q_i\delta q_j\left(\frac{^2\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪,a)}{q_iq_j}\right)_{𝐪=𝐪_s}}}`$ (63)
$`=`$ $`F(a)(4\pi \nu )^{3/2}{\displaystyle \underset{s}{}}j_s(𝐱,𝐪_s,a)e^{\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪_s,a)/2\nu },`$
where

$$j_s(𝐱,𝐪_𝐬,a)=\left(\mathrm{det}\left[\frac{^2\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪,a)}{q_iq_j}\right]_{𝐪=𝐪_s}\right)^{1/2}.$$
(64)
The sum in eq. (63) extends over the saddle points $`𝐪_s`$, i.e. the Lagrangian coordinates corresponding to the absolute minima of the function $`\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪,a)`$, for fixed $`𝐱`$ and $`a`$. These are found by solving the equation

$$_𝐪\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪,a)|_{𝐪=𝐪_s}=_𝐪\left(S_{\mathrm{cl}}(𝐱,𝐪,a)+\mathrm{\Phi }_0(𝐪)\right)|_{𝐪=𝐪_s}=0.$$
(65)
Since, for the classical action, $`_𝐪S_{\mathrm{cl}}=𝐮_0(𝐪)`$, we see that the classical orbit (56) passes through the saddle point $`𝐪_s`$ if the initial velocity satisfies $`𝐮_0(𝐪)=_𝐪\mathrm{\Phi }_0(𝐪)`$ and $`\mathrm{\Phi }(𝐱,𝐪,a)\mathrm{\Phi }(𝐱,𝐪_s,a)`$, for any $`𝐪𝐪_s`$.

We then obtain

$$𝐱(𝐪_s,a)=𝐪_s+a_𝐪\mathrm{\Phi }_0(𝐪_s)_0^a𝑑\tau _0^\tau 𝑑\tau ^{}_𝐱\eta (𝐱_{\mathrm{cl}}(𝐪_s,\tau ^{}),\tau ^{}),$$
(66)
for the particle trajectory, which coincides with eq. (36), as $`\mathrm{\Phi }_0=\phi _0`$.

Replacing eq. (63) in eq. (51) we finally obtain our saddle-point solution of the forced Burgers equation for the baryon velocity field

$$𝐮(𝐱,a)=\underset{s}{}w_s(𝐱,𝐪_s,a)_𝐱S_{\mathrm{cl}}(𝐱,𝐪_s,a),$$
(67)
where the (normalized) weights read

$$w_s(𝐱,𝐪_s,a)=\frac{j_s(𝐱,𝐪_s,a)e^{\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪_s,a)/2\nu }}{_sj_s(𝐱,𝐪_s,a)e^{\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪_s,a)/2\nu }}=\frac{j_s(𝐱,𝐪_s,a)}{_sj_s(𝐱,𝐪_s,a)}$$
(68)
(the last step is justified by the fact that, for a given $`𝐱`$ and $`a`$, all the absolute minima have the same $`\mathrm{\Phi }_{\mathrm{cl}}`$).

Since the Eulerian gradient of the action is the velocity along the classical trajectory, we have from eq. (56)

$$_𝐱S_{\mathrm{cl}}(𝐱,𝐪_s,a))=𝐮_{\mathrm{cl}}(𝐱(𝐪,a),a)=\frac{𝐱𝐪_s}{a}_0^ad\tau _𝐱\eta (𝐱(𝐪_s,\tau ),\tau )+_0^ad\tau _0^\tau d\tau ^{}_𝐱\eta (𝐱(𝐪_s,\tau ^{}),\tau ^{}).$$
(69)
At early times, there is a single contribution to the steepest-descent velocity field, coming from a unique Lagrangian saddle point $`𝐪_s`$ that minimizes $`\mathrm{\Phi }_{\mathrm{cl}}`$ for given $`𝐱`$ and $`a`$, so that $`𝐮(𝐱,a)=𝐮(𝐱_{\mathrm{cl}}(𝐪_s,a),a)`$. With time passing the mapping $`𝐪_s𝐱`$ becomes many-to-one. At this stage, according to the steepest-descent solution, the flow at $`𝐱`$ becomes a weighted average of the velocity of all the mass elements converging to that point at $`a`$ along their classical orbits. This sets the onset of the epoch of pancake formation in our model. This interpretation of our solution will become more clear from the geometrical construction outlined in the next Section.

In concluding this Section, let us come back to our analogy with the KPZ model. A different version of it, which is also often used in condensed matter physics, is obtained by making the Hopf-Cole nonlinear transformation on the KPZ equation, which leads to the linear diffusion equation for the free energy $`𝒵=\mathrm{exp}(\mathrm{\Phi }/2\nu )`$ (Kardar & Zhang 1987). Its path-integral representation is then interpreted as follows. The path traced by a particle starting from an initial point $`𝐱(0)=𝐪`$ and arriving at a final point $`𝐱(a)=𝐱`$, moving in $`D`$ spatial dimensions, can be viewed as a polymer in $`D+1`$ dimensions joining $`(𝐪,0)`$ to $`(𝐱,a)`$, with the time variable interpreted as a spatial one along the direction of main extension of the polymer. One specifically uses the term directed polymer to emphasize that the path is constrained to go only forward in time without crossing itself. The case studied here of vanishing viscosity is frequently referred to, in the language of KPZ and directed polymers, as the ‘strong coupling’, or ‘zero-temperature’, limit. This coincides with the limit of large Reynolds number, in the language of turbulence, or to the classical limit, in our quantum mechanical analog.

## 6 Drawing the skeleton of the IGM distribution

There are many possible ways to look for a numerical solution of the forced Burgers equation. Interesting techniques might be devised, based on the analogy with the KPZ model. In the $`\nu 0`$ limit, in fact, the problem of finding the ‘free energy’ $`𝒵`$, corresponding to our potential $`𝒰`$, reduces to a simple optimization problem. In the KPZ context, one frequently uses a numerical scheme referred to as the transfer matrix method (Kardar & Zhang 1987). First, one generates a random energy landscape (e.g. in 2D, on a square lattice, one puts random numbers with a certain distribution and correlation, depending on the statistical properties of the noise $`\eta `$). The optimal particle starts at a lattice point $`𝐪`$, chooses the next lattice point which would cost it the least energy, jumps to it and so on, continuing its journey to the final point $`𝐱`$. In this way an optimal polymer, joining $`(𝐪,0)`$ to $`(𝐱,a)`$, is formed with the minimum amount of energy. The free energy $`𝒵`$ is then evaluated by weighing the optimal polymer, which ends there, over all the possible starting points. A 1D illustration of this numerical method is given in Figure 6.

The application of a similar scheme to the cosmological structure formation problem is however complicated by the colored and non-stationary time-correlation properties of our random potential, as well as by the fact that the cosmological interesting regime does not usually occur at asymptotically large times, unlike what happens in the case of directed polymers.

Let us next describe an alternative method, based on our steepest-descent solution of the forced Burgers equation. Two main technical difficulties arise, when one tries to implement the saddle-point solution to simulate the baryon distribution. The first concerns the velocity along the individual classical trajectories, which contribute to the sum in eq. (67): this involves a step-by-step integration of the force along the trajectory. The second concerns the weight $`w_s`$ which requires explicit knowledge of $`S_{\mathrm{cl}}`$ along these classical trajectories. The problem is equivalent to that of finding the velocity potential, $`\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪,a)=S_{\mathrm{cl}}(𝐱,𝐪,a)+\mathrm{\Phi }_0(𝐪)`$, which solves the ‘classical’ Hamilton-Jacobi equation (37) with initial conditions $`\mathrm{\Phi }_{\mathrm{cl}}(𝐪,𝐪,0)=\mathrm{\Phi }_0(𝐪)`$. The technique is therefore completely analogous to that used in the free adhesion approximation for DM particles, but is made more involved by the presence of the random potential $`\eta `$.

A substantial simplification is obtained, if we adopt the same approximation of Section 4.1, which basically consists in expanding the particle orbits to first order in the displacement vector. Let us describe the method in detail. For each orbit the system contains a cyclic variable corresponding to the initial velocity, namely

$$𝐮_0=\frac{𝐱𝐪}{a}+_𝐪(\psi \phi _0),$$
(70)
where $`\psi `$ is the baryon potential defined in eq. (41). Replacing $`\dot{𝐱}=𝐮_0+a_𝐪(\phi _0\psi )`$ in the action (53) and then using (70), we obtain, to second order in the displacement vector $`𝐒`$,

$$\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪,a)=\frac{(𝐱𝐪)^2}{2a}a(𝐱𝐪)_𝐪\dot{\psi }(𝐪,a)\psi (𝐪,a)a\dot{\psi }(𝐪,a)\frac{1}{2}_0^a𝑑\tau \tau ^2\left(_𝐪\dot{\psi }(𝐪,\tau )\right)^2+𝒪(𝐒^3),$$
(71)
where, for consistency, we expanded the potential $`\eta `$ to first order in the displacement vector around the Lagrangian coordinate (as $`\eta `$ itself is a first-order quantity), i.e.

$$\eta (𝐱,a)=\eta (𝐪,a)+_𝐪\eta (𝐪,a)\left[a𝐮_0(𝐪)_𝐪\left(\psi (𝐪,a)\phi _0(𝐪)\right)\right]+𝒪(𝐒^3),$$
(72)
It is immediate to check that the velocity potential of eq. (71) gives the required solution of the Hamilton-Jacobi equation (37), if the potential $`\eta `$ in its RHS is consistently expanded as in eq. (72). From eq. (71) we can immediately obtain an explicit expression for $`j_s`$ which appears in our saddle-point solution through the weight of eq. (68). We have

$$j_s(𝐪_𝐬,a)=\left(\mathrm{det}\left[\frac{\delta _{ij}}{a}\frac{^2\psi (𝐪,a)}{q_iq_j}+a\frac{^2\dot{\psi }(𝐪,a)}{q_iq_j}\right]_{𝐪=𝐪_s}\right)^{1/2}+𝒪(𝐒^2).$$
(73)
The Lagrangian approximation for the classical trajectory

$$𝐱_{\mathrm{cl}}(𝐪_s,a)=𝐪_sa_𝐪\psi (𝐪_s,a)$$
(74)
and the corresponding velocity

$$𝐮_{\mathrm{cl}}(𝐱,𝐪_s,a)=\frac{𝐱𝐪_s}{a}a_𝐪\dot{\psi }(𝐪_s,a)$$
(75)
are the remaining ingredients of the saddle-point solution (67). Note that everything is now explicitly given in terms of the linear baryon potential $`\psi `$ and its time derivatives.

We are now ready to introduce a geometrical construction of the saddle-point solution, obtained in analogy with the adhesion approximation, which aims at obtaining the skeleton of the IGM distribution. The graphic technique has been largely applied in connection with the solution of the one-dimensional Burgers equation (e.g. Burgers 1974). In the cosmological context this method has been described in a number of papers (e.g. Gurbatov et al. 1985, 1989; Kofman & Shandarin 1988; Shandarin & Zel’dovich 1989; Vergassola et al. 1994; Sahni & Coles 1995) and applied in one, two and three-dimensional numerical simulations of the large-scale structure of the Universe (Kofman, Pogosyan & Shandarin 1990; Williams et al. 1991; Kofman et al 1992; Sahni, Sathyaprakash & Shandarin 1994; Sathyaprakash et al. 1995), finding good agreement with N-body simulations in models which have sufficient large-scale power (e.g. Sathyaprakash et al. 1995).

Let us start by noting that, for given $`𝐱`$ and $`a`$, the classical velocity potential $`\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪,a)`$ describes a random three-dimensional hypersurface in Lagrangian space. The idea is that absolute minima of this hypersurface correspond to Lagrangian points of first-touching with the hyperplane $`\mathrm{\Sigma }_h(𝐪)=h`$, as the parameter $`h`$ is raised from $`\mathrm{}`$ <sup>5</sup><sup>5</sup>5In the free adhesion case it is actually more convenient to look for the first-touching points of the paraboloid $`𝒫_h(𝐪|𝐱,a)=(𝐪𝐱)^2/2a+h`$, with the random hypersurface $`\mathrm{\Phi }_0(𝐪)`$, as the latter is both $`𝐱`$ and $`a`$ independent..

The geometrical construction can proceed along similar lines as for the adhesion model, namely:

i) At each time $`a`$, and for each point $`𝐱`$ in Eulerian space, find the hyperplane $`\mathrm{\Sigma }_h(𝐪)`$ tangential to $`\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪,a)`$ at the Lagrangian point(s) $`𝐪_{s_i}`$, by increasing the parameter $`h`$ from $`\mathrm{}`$, so that the two hypersurfaces do not cross anywhere.

ii) Once the saddle point(s) $`𝐪_{s_i}(𝐱,a)`$ are found, the corresponding velocity at $`𝐱`$ at time $`a`$ is given by eq. (67), as the weighted sum of single-particle velocities (75), with the weights given by eq. (68), with $`j_s`$ from eq. (73) and $`S_{\mathrm{cl}}`$ from eq. (71).

At early times, the hypersurface $`\mathrm{\Phi }_{\mathrm{cl}}`$ is a narrow, distorted paraboloid whose apex starts in $`𝐪=𝐱`$ and gradually moves away from it. At these early times, the curvature of this random hypersurface is large, so that there is a one-to-one correspondence between Lagrangian and Eulerian space. As time passes, our classical velocity potential becomes shallower, so that it becomes increasingly frequent to find points $`𝐱`$ such that the first-touching condition of $`\mathrm{\Sigma }_h`$ and $`\mathrm{\Phi }_{\mathrm{cl}}`$ occurs at multiple saddle points $`𝐪_{s_i}`$. The Lagrangian to Eulerian correspondence has now become many-to-one, with the degree of the mapping being inversely proportional to the dimensionality of the structure formed. Two simultaneous points of touching, $`𝐪_{s_1}(𝐱,a)`$ and $`𝐪_{s_2}(𝐱,a)`$, signal the formation of a sheet, or pancake, in the corresponding Eulerian location $`𝐱`$, three indicate the presence of a filament, and four that of a knot. As is well known from the adhesion model, filaments occur at the intersection of pancakes, while nodes at the intersection of filaments.

There is only a caveat in our geometrical construction, owing to the approximate form of eq. (71). Because of this, eq. (65) is only approximately satisfied, namely $`_𝐪\mathrm{\Phi }_{\mathrm{cl}}`$ calculated from eq. (71) contains spurious terms which are of second-order in the displacement vector. As a consequence, our accuracy in the determination of the saddle points is limited to first order. Nevertheless, our Lagrangian approximation is likely to be quite accurate not only at early times, when the evolution is linear, but also at later times, as far as the skeleton of the IGM structure is concerned. In fact, the dominant contribution to our saddle-point solution in denser regions, such as pancakes, filaments and knots, is expected to come from the deepest potential wells \[see, e.g., the related interpretation of the path-integral solution by Zel’dovich et al. (1985, 1987)\], i.e. precisely from those Lagrangian sites where the particles experience the smallest displacement from their initial positions (e.g. Matarrese et al. 1992).

As for the standard adhesion model, the above geometrical construction cannot describe the internal structure of dense regions, like pancakes, filaments and knots, as it is only designed to find the sites where these structure form and to follow the subsequent evolution through a continuous merging process (e.g. Shandarin & Zel’dovich 1989). From this point of view, it is clear that a numerical algorithm based on the direct integration of the solution in eq. (61), would be far more suitable to follow the baryon distribution on smaller scales. This kind of numerical integration can be performed by extending the technique used by Weinberg & Gunn (1990a,b) and Melott, Shandarin & Weinberg (1994). For the free adhesion approximation, Weinberg & Gunn notice that the exact solution \[our eq. (95)\] appears as a Gaussian convolution integral, which is easy to handle numerically. For our stochastic adhesion model, this type of calculation is obviously made more complex by the presence of the random force term, which does not generally allow to get an exact explicit expression for the action. This problem, however, is enormously simplified if one adopt the approximate expression of eq. (71), for the classical action. As mentioned previously, the last step in this scheme would consist in the numerical integration of eq. (62) to obtain the actual particle positions. Although the computational time required to perform the integration in eq. (61) plus that of eq. (62) implies only a modest speed up, compared to a particle-mesh N-body code; one should keep in mind that our stochastic adhesion model actually aims at describing a physical situation where the exact treatment of the two-fluid evolution would require a full hydrodynamical code.

## 7 Statistics of the density field

The stochastic adhesion model provides a direct insight on the properties of the velocity field away from the linear regime. A vast literature is devoted to the study of the statistical properties of this field, in cases when the external potential $`\eta `$ is a Gaussian random field which is either time-independent (e.g. Zel’dovich et al. 1987), or with white-noise time-correlation properties (e.g. Bouchaud et al. 1995; E et al. 1997). It is generally found that the system develops intermittency at late times. Based on this property, Jones (1999) argues that the velocity field is normally distributed even during the nonlinear stage. He also suggests that this may provide a dynamical motivation for a Lognormal PDF of the density field in the baryonic component.

More complex is to analyze the statistical behavior of the velocity and density field in our model, owing to the non-trivial time-correlation properties of the noise, as well as to the fact that we are interested in the behavior of the system in the intermediate regime, when most of the matter still resides in pancakes and filaments, rather than having relaxed in the most massive clumps, corresponding to the nodes of the cellular structure. We will only give here a preliminary introduction to this problem, focusing on the density field, rather than the peculiar velocity.

A convenient form of the comoving local matter density can be obtained by integrating the continuity equation along the particle trajectories. One gets (e.g. Matarrese et al. 1992)

$$1+\delta (𝐱(𝐪,a),a)=\mathrm{exp}\left[_0^a𝑑\tau \theta (𝐱(𝐪,\tau ),\tau )\right],$$
(76)
where $`\theta (𝐱,a)_𝐱𝐮(𝐱,a)`$ and the time integration is along the particle trajectory. The validity of this integral expression for the density is restricted to the single-stream regime. We are then allowed to apply it in our case, provided the viscosity coefficient $`\nu `$ takes a finite value. For each Eulerian position $`𝐱`$ at time $`a`$ there is, in fact, a unique Lagrangian element $`𝐪`$ implicitly defined through eq. (62).

A relevant statistical information on the density field could then be obtained by studying its (disconnected) moments, $`(1+\delta )^p`$. We have

$$\left(1+\delta (𝐱,a)\right)^p_\mathrm{E}=\mathrm{exp}\left[(p1)_0^a𝑑\tau \theta (𝐱(𝐪,\tau ),\tau )\right]_\mathrm{L},$$
(77)
where the subscripts ‘E’ and ‘L’ denote the statistical ensemble over which the average is performed: the one on the LHS involves the Eulerian PDF, whereas that on the RHS is most conveniently obtained with the Lagrangian PDF. To account for this difference we multiplied the RHS by the Jacobian $`𝐱/𝐪=(1+\delta )^1`$ (e.g. Kofman et al. 1994). Note that the Lagrangian-to-Eulerian mapping is always one-to-one, as the velocity field solving the forced Burgers equation is non-singular everywhere for any finite value of $`\nu `$; thus the actual fluid elements, which move according to eq. (62), do not experience orbit-crossing during their evolution.

The main technical difficulty in dealing with eq. (76), or any equivalent expression for the density field, resides in the inversion of the Lagrangian map $`𝐪𝐱(𝐪,a)`$ (e.g. Vergassola et al. 1994). The same problem, of course, makes it generally difficult to deal analytically with its statistical properties (such as the disconnected moments above). An alternative technique to reconstruct the matter density field from the velocity, solving the Burgers equation, is discussed by Gurbatov (1996; see also Gurbatov, Malakhov & Saichev 1991), who propose a modification of the continuity equation to simplify this inversion problem. Another interesting approach to the same problem is adopted by Vergassola et al. (1994), who implement a Fast Legendre Transform algorithm. Here we will sketch the first steps of a different method, which aims at relating the mass density in $`𝐱`$ at time $`a`$ to the properties of the classical velocity field at the saddle-points $`𝐪_s(𝐱,a)`$.

### 7.1 The velocity divergence

The first step to obtain the density field is to write a suitable expression for the velocity divergence $`\theta `$. Using eq. (51), we can write the velocity divergence in the general form

$$\theta =\frac{𝐮^2}{2\nu }2\nu \frac{^2𝒰}{𝒰},$$
(78)
Next, we can replace for $`𝐮`$ and $`𝒰`$ the expressions of eq. (61) and (59) (which would be exact in the free adhesion case), which gives

$$\theta (𝐱,a)=\frac{3}{a}\frac{1}{2\nu }[[\left(𝐮_{\mathrm{cl}}(𝐱,𝐪,a)[[𝐮_{\mathrm{cl}}(𝐱,𝐪,a)]]\right)^2]],$$
(79)
where we introduced the symbol

$$[[𝐀(𝐱,𝐪,a)]]\frac{d^3q𝐀(𝐱,𝐪,a)e^{\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪,a)/2\nu }}{d^3qe^{\mathrm{\Phi }_{\mathrm{cl}}(𝐱,𝐪,a)/2\nu }},$$
(80)
for any scalar, vector or tensor $`𝐀`$. To get the first term on the RHS of eq. (79) we used the approximate form for the classical action deriving from eq. (71).

To better understand the meaning of the expression (79) for the velocity divergence, we can consider the small $`\nu `$ limit, using the steepest descent approximation. We get

$$\theta (𝐱,a)=\underset{s}{}w_s(𝐪_s,\tau )\theta _{\mathrm{cl}}(𝐱(𝐪_s,a),a)\frac{1}{2\nu }\delta 𝐮^2(𝐱,a),$$
(81)
where $`\theta _{\mathrm{cl}}_𝐱𝐮_{\mathrm{cl}}=_𝐱^2\mathrm{\Phi }_{\mathrm{cl}}`$, and

$$\delta 𝐮^2(𝐱,a)=\underset{s}{}w_s(𝐪_s,a)𝐮_{\mathrm{cl}}^2(𝐪_s,a)\underset{s}{}\underset{s^{}}{}w_s(𝐪_s,a)w_s^{}(𝐪_s^{},a)𝐮_{\mathrm{cl}}(𝐪_s,a)𝐮_{\mathrm{cl}}(𝐪_s^{},a),$$
(82)
from which it becomes evident its role as mean square velocity dispersion caused by the convergence of several classical streams in $`𝐱`$. Note that, by definition, $`\delta 𝐮^20`$. In order to obtain the expression of eq. (81), we expanded the classical velocity to first order around the saddle point,

$$u_{\mathrm{cl}}^i(𝐱,𝐪,a)=u_{\mathrm{cl}}^i(𝐱,𝐪_s,a)\frac{1}{a}\left(\delta _j^i+a^2\frac{^2\dot{\psi }(𝐪_s,a)}{q_iq^j}\right)\left(q^jq_s^j\right)$$
(83)
and wrote the classical velocity divergence in the form

$$\theta _{\mathrm{cl}}(𝐱(𝐪,a),a)=\frac{3}{a}\frac{1}{a}\left[\delta _{ij}+a^2\frac{^2\dot{\psi }(𝐪,a)}{q^iq^j}\right]\left[\delta ^{ij}a\frac{^2\psi (𝐪,a)}{q_iq_j}\right]^1,$$
(84)
which follows from taking the Eulerian divergence of eq. (75). Note that, while in the free adhesion case, eq. (81) is the direct result of the saddle-point approximation, in our stochastic adhesion model, that expression has been derived only to first order in the displacement vector. Nonetheless, we believe that the validity of eq. (81) is much more general than its approximate derivation might suggest.

### 7.2 Evolution of the density field

We are now in a position to discuss the evolution of the density field in the various stages of the structure formation process.

#### 7.2.1 Linear stage

At early times, when the system is well approximated by Eulerian first-order perturbation theory, $`_0^a𝑑\tau \theta a^2\psi 1`$, and the density fluctuation is linearly related to that of the initial gravitational potential, which we assumed to be a uniform Gaussian process. Density fluctuations at this time are also Gaussian distributed.

#### 7.2.2 Weakly nonlinear stage

Later on, the system enters a weakly nonlinear regime, when density fluctuations are no longer small, but caustic formation is still a sporadic event. In this laminar stage the evolution of the system is generally well described by Lagrangian first-order perturbation theory (see Section 4.1). Our solution (81) yields

$$1+\delta (𝐱(𝐪,a),a)=\mathrm{exp}\left(_0^a𝑑\tau \theta _{\mathrm{cl}}(𝐱(𝐪_s,\tau ),\tau )\right),$$
(85)
as, for each $`𝐱`$, there is a unique steepest-descent solution, corresponding to the particle, with initial coordinate $`𝐪_s`$, which reaches that point at time $`a`$. This expression is of course equivalent to the more conventional one of Section 4.1.2 \[e.g. Appendix A in (Matarrese et al. 1992)\]. At this early nonlinear stage, the density distribution is no longer Gaussian. Obtaining the PDF of the mass density at this stage is a classical problem which has been solved long ago by Doroshkevich (1970), and analyzed in more detail by several authors (e.g. Kofman et al. 1994; Bernardeau & Kofman 1994), for the collisionless case. As anticipated in Section 4.1.2, the presence of non-zero acceleration in the motion of our collisional fluid elements does not imply relevant modifications compared with the DM case, so, for instance, the functional form of the PDF is unchanged.

#### 7.2.3 Mildly nonlinear stage

Let us finally come to the most interesting issue: the analysis of the mildly nonlinear stage, when dense, spatially extended structures have grown in the system. It is this process which is accompanied by the occurrence of strong phase coherence both in the velocity and density fields. According to our stochastic adhesion model, the formation of a shock singularity in $`𝐱`$ at time $`a`$ is accompanied by the occurrence of multiple saddle-point solutions $`𝐪_s(𝐱,a)`$. As a consequence of these multiple solutions, the velocity dispersion term of eq. (82) will suddenly become non-zero. Because of the ‘diverging’ (as $`\nu 0`$) factor $`1/2\nu `$, however, any non-zero $`\delta 𝐮^2`$ would signal a ‘singularity’ (for $`\nu 0`$) in the local density. One has, therefore, the following picture: as $`\nu 0`$ the mass tends to get concentrated in the cellular structure, which becomes a thin cobweb of delta-like sheets, filaments and nodes. On the other hand, as soon as the Lagrangian-to-Eulerian mapping becomes many-to-one our expression (76) for the mass density will lose its validity, and a more clever technique should be used to account for the actual mass converging into these nonlinear structures. A possible hint in this direction might be given by the calculation of the velocity field on a pancake, made by Kofman (1989), within the free adhesion model. A more detailed analysis of this problem will be presented elsewhere.

## 8 Conclusions

In this paper we have proposed a new dynamical model for the formation of mildly nonlinear structure in the IGM, based on a ‘stochastic adhesion’ approximation to the Navier-Stokes equation, which governs the evolution of the collisional baryon component. The main idea of our model is that a random external force can be used to approximate the composition of Hubble drag, local gravitational force and pressure gradients. Such a force, consistently calculated from first-order perturbation theory, has two main effects. First, it produces a slight distortion of the large-scale filamentary cosmic web that characterizes the IGM distribution, compared with the DM one. Second, it adds a small-scale noise to the otherwise deterministic evolution of the fluid, thereby imprinting characteristic ripples on top of the cellular network. This small-scale roughness, which is the memory of acoustic oscillations in the baryon fluid, can be used as a probe of the IGM equation of state and thermal history after hydrogen reionization.

This paper was mostly devoted to the presentation of the stochastic adhesion model. A few applications were shown, only for the simplified inviscid ($`\nu =0`$) case, where the trajectories followed by patches of the collisional fluid are obtained through a first-order Lagrangian approximation scheme, somehow similar to that adopted by Gnedin & Hui (1996) and Hui, Gnedin & Zhang (1997). Already at this level, however, the precise IGM thermal history is found to affect the resulting clustering properties of the gas, through the detailed structure of denser regions. The second part of our paper was devoted to the solution of the forced Burgers equation, which governs the evolution of the IGM peculiar velocity field in our model. Among our main results is the saddle-point solution, eq. (67), and the Lagrangian approximation of Section 6: both are completely new and might have much wider applications than the present context. We also sketched a geometrical representation of our saddle-point solution, analogous to that used for the standard adhesion model, which allows to draw the skeleton of the IGM distribution.

We did not address here the important issue of whether pancakes or filaments provide the dominant structure in the large-scale matter distribution. According to the cosmic web description of Bond et al. (1996) the nonlinear evolution of the DM brings a filamentary network into relief. Bond & Wadsley (1997) applied similar ideas to the IGM distribution probed by the Ly$`\alpha `$ forest. Quite recently, Colombi, Pogosyan & Souradeep (2000) adopted topological descriptors to conclude that the structure of the Universe at the percolating threshold is predominantly filamentary, which confirms the results of a percolation analysis of N-body simulations by Sathyaprakash, Sahni & Shandarin (1996). In both the free and stochastic adhesion models, this aspect of the emerging structures on large scales is expected to be strictly related to the stage of dynamical evolution of the system, as well as to the value of the viscosity coefficient $`\nu `$ and to the resolution scale at which the structure is analy

The formation of a cellular structure on large scales, appearing as a cobweb of interconnected sheets, filaments and nodes, both for the DM and IGM components, is accompanied by the growth of specific phase coherence in the velocity and density field. We have shown how this phenomenon, which has been dubbed ‘intermediate intermittency’ (Zel’dovich et al. 1985, 1987), is most naturally described in terms of the forced Burgers equation, which underlies our model. A wavelet analysis of the Ly$`\alpha `$ forest, such as that recently initiated by Meiksin (2000) and Theuns & Zaroubi (2000), would be the ideal tool for a statistical analysis of this phenomenon, even on those small scales where the effects of the thermal gas pressure are expected to imprint characteristic features on the baryon distribution, which should tell us about the equation of state and detailed thermal history of the IGM.

In concluding this paper let us emphasize that the stochastic adhesion model presented here might represent the ideal tool to study the statistical properties of the low-column density Ly-alpha forest, because of its ability to accurately describe the IGM distribution from the largest scales, poorly probed by present-day hydro-simulations, down to those regions which have undergone mildly non-linear evolution (baryon density contrast up to $`510`$), which cannot be properly represented by smoothed versions of the Zel’dovich approximation.

## Acknowledgments.

We would like to thank U. Frisch, L. Moscardini, A. Nusser, A. Stella and M. Viel for enlightening discussions and for technical help. RM acknowledges the TMR Network FMRXCT980183 for financial support.

## Appendix A First-order Lagrangian perturbation theory for the baryons

We will briefly derive here an expression for the baryon trajectories within first-order Lagrangian perturbation theory. Lagrangian approximation schemes have been successfully applied to the weakly nonlinear evolution of the DM fluid (e.g. Buchert 1989, 1992; Moutarde et al. 1991; Bouchet et al. 1992; Catelan 1995). Adler & Buchert (1999) have recently applied the Lagrangian perturbation technique to the case of a collisional self-gravitating fluid, corresponding to $`f_\mathrm{b}=1`$ in our equation (7). The case discussed here is instead that of a collisional fluid moving in the gravitational field caused by a DM component, corresponding to $`f_\mathrm{b}=0`$ in our equation (7).

The particle trajectories of both components can be represented, before shell-crossing, in the general form

$$𝐱(𝐪,a)=𝐪+𝐒(𝐪,a)$$
(86)
Mass conservation (starting from a uniform distribution) implies

$$1+\delta (𝐱(𝐪,a),a)=\frac{1}{\mathrm{det}(x_\alpha /q_\beta )}=\frac{1}{J}=1_q𝐒+O(S^2),$$
(87)
where $`J`$ is the Jacobian determinant of the transformation from Lagrangian to Eulerian coordinates.

Taking the divergence of the baryon Euler equation we obtain

$$_𝐱\ddot{𝐒}_\mathrm{b}=\frac{3}{2a}[_𝐱\dot{𝐒}_\mathrm{b}+\frac{\delta _{\mathrm{DM}}}{a}+\frac{1}{(\gamma 1)ak_J^2}_𝐱^2(1+\delta _\mathrm{b})^{\gamma 1}]$$
(88)
Substituting for the overdensities from eq. (87) we have

$$_𝐱\ddot{𝐒}_\mathrm{b}+\frac{3}{2a}_𝐱\dot{𝐒}_\mathrm{b}=\frac{3}{2a^2}\frac{1J_{\mathrm{DM}}}{J_{\mathrm{DM}}}\frac{3}{2(\gamma 1)a^2k_J^2}_𝐱^2\left(\frac{1J_\mathrm{b}}{J_\mathrm{b}}\right)^{\gamma 1}$$
(89)
Next, we expand everything to first order in the displacement vector $`𝐒`$ and obtain

$$_𝐪\ddot{𝐒}_\mathrm{b}+\frac{3}{2a}_𝐪\dot{𝐒}_\mathrm{b}\frac{3}{2a^2}_𝐪𝐒_{\mathrm{DM}}=\frac{3}{2a^2k_J^2}_𝐪^2_𝐪𝐒_\mathrm{b}.$$
(90)
Note that at a given time, the same Eulerian position $`𝐱`$ is generally reached by the two components (DM and baryons) starting from different Lagrangian positions $`𝐪_{\mathrm{DM}}`$ and $`𝐪_\mathrm{b}`$. To lowest order in the displacement vector, however, we are allowed to ignore this difference and set simply $`𝐪=𝐪_{\mathrm{DM}}=𝐪_\mathrm{b}`$.

We now assume that the flow is irrotational $`_𝐱\times \dot{𝐒}=0`$, which, to lowest order implies $`𝐒=_𝐪\mathrm{\Psi }`$, and reduces (90) to the scalar differential equation

$$\ddot{\mathrm{\Psi }}_\mathrm{b}+\frac{3}{2a}\dot{\mathrm{\Psi }}_\mathrm{b}\frac{3}{2a^2}\frac{1}{k_J^2}_q^2\mathrm{\Psi }_\mathrm{b}=\frac{3}{2a^2}\mathrm{\Psi }_{\mathrm{DM}},$$
(91)
which can be solved in Fourier space. Using the well-known DM solution $`\mathrm{\Psi }_{\mathrm{DM}}=a\phi _0`$, we obtain

$$\ddot{\mathrm{\Psi }}_\mathrm{b}+\frac{3}{2a}\dot{\mathrm{\Psi }}_\mathrm{b}+\frac{3}{2a^2}\frac{k^2}{k_J^2}\mathrm{\Psi }_\mathrm{b}=\frac{3}{2a}\phi _0,$$
(92)
which can be solved in exactly the same manner as we previously solved (16) to obtain the baryon overdensity. The solution can be expressed in the form $`\mathrm{\Psi }_\mathrm{b}=a\psi _\mathrm{b}`$, where $`\psi _\mathrm{b}(𝐤,a)=W_\mathrm{b}(k,a)\phi _0(𝐤)`$ with $`W_\mathrm{b}`$ the IGM linear filter function obtained in the main text.

## Appendix B Path-integral solution for the free adhesion model

In the case of zero potential $`\eta `$, eq. (47) reduces to the Burgers equation. By dropping the potential $`\eta `$ we obtain the free particle trajectory

$$𝐱(𝐪,a)=𝐪+a𝐮_0(𝐪)$$
(93)
and the classical action reduces to

$$S_{\mathrm{cl}}(𝐱,𝐪,a)=_0^a𝑑\tau \frac{\dot{𝐱}^2(𝐪,\tau )}{2}=\frac{a𝐮_0^2(𝐪)}{2}=\frac{(𝐱𝐪)^2}{2a}.$$
(94)
We therefore have

$$𝒰(𝐱,a)=F(a)d^3qe^{(S_{\mathrm{cl}}(𝐱,𝐪,a)+\mathrm{\Phi }_0(𝐪))/2\nu }=\frac{1}{(4\pi \nu a)^{3/2}}d^3qe^{\frac{(𝐱𝐪)^2}{4\nu a}\mathrm{\Phi }_0(𝐪)/2\nu },$$
(95)
where the pre-factor is evaluated by performing the Gaussian integration over $`\xi `$ in eq. (58).

The above expression (95) is equivalent to the standard expression for the adhesion approximation (e.g. Shandarin & Zel’dovich 1989).

In the steepest descent approximation, valid in the limit of vanishing viscosity, the particle trajectories are given by the solution of

$$_𝐪\left(\mathrm{\Phi }_0+S_{cl}\right)=0,$$
(96)
that is

$$_𝐪\mathrm{\Phi }_0(𝐪)+\frac{𝐱𝐪}{a}=0,$$
(97)
which coincides with eq. (93), for our set of initial conditions. The results of this Appendix also show that the quadratic approximation is exact for free particles.