text:
# Modelling the IGM and the Ly𝛼 forest at high redshift from the dark matter distribution

## 1 Introduction

The low-density Intergalactic Medium (IGM) offers a unique and powerful probe of the high redshift Universe. Numerous weak absorption lines seen in the spectra of distant quasars (the Ly$`\alpha `$ forest) are produced by the small residue of neutral hydrogen in filamentary structures intersected by the line of sight (Bahcall & Salpeter 1965; Gunn & Peterson 1965; see Rauch 1998 for a recent review). Such structures arise naturally in hierarchical cold dark matter dominated models of structure formation, in which the IGM is highly ionised by the UV-background produced by stars and galaxies. Simulations show that the low-column density ($`\mathrm{N}_{\mathrm{HI}}10^{14.5}`$ cm-2) absorption lines are produced by small fluctuations in the warm ($`T10^4`$ K) photo-heated IGM, which smoothly traces the mildly non-linear dark matter filaments and sheets on scales larger the Jeans scale (Cen et al. 1994, Petitjean et al. 1995, Miralda-Escudé et al. 1996, Zhang et al. 1998).

This picture is supported by analytical studies based on simple models for the IGM dynamics. Such models are based on either a local non-linear mapping of the linear density contrast, obtained for example by applying a lognormal transformation (Coles & Jones 1991) to the IGM (Bi, Börner & Chu 1992; Bi 1993; Bi, Ge & Fang 1995; Bi & Davidsen 1997, hereafter BD97), or on suitable modifications of the Zel’dovich approximation (Zel’dovich 1970), to account for the smoothing caused by the gas pressure on the Jeans scale (Reisenegger & Miralda-Escudé 1995; Gnedin & Hui 1998; Hui, Gnedin & Zhang 1997; Matarrese & Mohayaee 2002). Many properties of the Ly$`\alpha `$ absorbers can be understood with straightforward physical arguments (Schaye 2001, Zhang et al. 1998).

The most convincing support for this picture comes, however, from the comparison of observed spectra with mock spectra computed from hydrodynamical numerical simulations (Cen et al. 1994; Zhang, Anninos & Norman 1995, 1997; Miralda-Escudé et al. 1996; Hernquist et al. 1996; Charlton et al. 1997; Theuns et al. 1998). These simulated spectra accurately reproduce many observed properties of the Ly$`\alpha `$ forest. Although there are still some discrepancies between observed and simulated spectra, especially in the Doppler parameters of the absorption lines (Theuns et al. 1998; Bryan et al. 1999; Meiksin et al. 2001).

Observationally, the unprecedented high resolution observations of the HIRES and UVES spectrographs on the Keck and VLT telescopes, as well as observations with the Hubble Space Telescope, have lead to great advances in the understanding of the Ly$`\alpha `$ forest. HIRES allowed to detect lines with column densities as low as $`\mathrm{N}_{\mathrm{HI}}10^{12}`$ cm-2, while HST made a detailed analysis of the low-redshift Ly$`\alpha `$ forest at $`z<1.6`$ possible. Recent results include limits on the baryon density (Rauch et al. 1997), the temperature and equation of state of the IGM (Schaye et al. 2000; Theuns et al. 2001; Ricotti, Gnedin & Shull 2000; Bryan & Machachek 2000, McDonald et al. 2001), the power spectrum of the density fluctuations (Croft et al. 1998, 2001, Gnedin & Hamilton 2001, Zaldarriaga et al. 2001), the geometry of the Universe (Hui et al. 1999; McDonald et al. 2001; Viel et al. 2002) and direct inversions of the density field (Nusser & Haehnelt 1999, 2000; Pichon et al. 2001). Kim et al. (2001) used high resolution VLT/UVES data, to make an extensive analysis of the Ly$`\alpha `$ forest in the redshift range $`1.5<z<4`$. Dobrzycki et al. (2001) presented results on the clustering and evolution of the lines at $`z<1.7`$ using the HST/FOS spectrograph.

Many aspects of the warm photo-ionised Intergalactic medium can be well modelled by hydrodynamical simulations. Hydro simulations are, however, still rather limited in dynamic range. The box size of hydro simulations which resolve the Jeans mass of the photo-ionised IGM with a temperature of $`10^4K`$ probe the large scale fluctuations of the density field rather poorly. This leads to uncertainties due to cosmic variance in the fluctuations on scales approaching the box size and missing fluctuations on scales larger than the simulation box. Because of limited computational resources it is also hardly possible to perform extensive parameter studies. In order to overcome these problems approximate methods for simulating the Ly$`\alpha `$ forest in QSO absorption spectra are often used (e.g. Gnedin & Hui 1998; Croft et al. 1998; Croft et al. 1999; Meiksin & White 2001). We test here a wide range of such approximate methods against a large high-resolution hydro-dynamical simulation.

The plan of the paper is as follows. In Section 2 we briefly describe the hydro-simulation used. Section 3 presents the lognormal model for the IGM and an improved model based on the implementation of the gas probability distribution function taken from the hydro-simulation. We improve our modelling in Section 4, using two different ways of filtering the linear dark matter density field to model pressure effects. In this Section we present a further improvement, based on modelling the gas distribution starting from the non-linear dark matter density field obtained from the hydro-simulation and we show how the smoothing on the Jeans length of the evolved dark-matter density field poorly reproduces the real gas distribution. In Section 5 we compare the different methods proposed in terms of the 1-point and 2-points PDF of the flux. Section 6 contains a general discussion and our conclusions.

## 2 The hydro-dynamical simulation

In the following Sections we will test our approximate schemes to model the Ly$`\alpha `$ forest against a large high-resolution hydro-dynamical simulation. We therefore summarise here some parameters of the hydro-dynamical simulation used. The simulation techniques are described in more detail in Theuns et al. (1998). We analyse a total of 7 outputs at redshifts $`z=49`$, $`z=10`$, $`z=4`$, $`z=3.5`$, $`z=3`$, $`z=2.25`$, $`z=2`$ of a periodic, cubic region in a $`\mathrm{\Lambda }`$CDM Universe. The cosmological parameters are: $`\mathrm{\Omega }_{0m}=0.3`$, $`\mathrm{\Omega }_{0\mathrm{\Lambda }}=0.7`$, $`H_0=100h`$ km s-1Mpc with $`h=0.65`$ and $`\mathrm{\Omega }_{0b}h^2=0.019`$. The comoving size of the box is $`12.0/h`$ Mpc. There are $`256^3`$ DM particles and $`256^3`$ gas particles, whose masses are $`m_{DM}=1.13\times 10^7M_{}`$ and $`m_{IGM}=1.91\times 10^6M_{}`$. The input linear power spectrum was computed with CMBFAST (Seljak & Zaldarriaga 1996), and normalised to the abundance of galaxy clusters using $`\sigma _8=0.9`$ (Eke, Cole & Frenk 1996), where $`\sigma _8`$ denotes the mass fluctuations in spheres of radius $`8h^1`$ Mpc.

The simulation code is based on HYDRA (Couchman et al. 1995) and combines Smoothed Particle Hydrodynamics (see e.g. Monaghan 1992) with P3M for self-gravity. Spline interpolation over gas particles allows the computation of smooth estimates for density, temperature and velocity and their gradients. The width of the spline kernel is matched to the local particle number density, in this way high density regions have higher numerical resolution than lower density ones. Photo-ionisation and photo-heating rates are computed using the fits in Theuns et al. (1998).

Along many randomly chosen sight lines parallel to one of the axis of the simulation box, we compute the gas and DM (over) density and peculiar velocity, the gas temperature, and the neutral hydrogen (density weighted) density, temperature and peculiar velocity. This allows us to compute the absorption spectra. The smooth estimate for the dark matter fields are computed using SPH interpolation as well, with a smoothing length chosen to given of order $`32`$ neighbour contributions per particle (Appendix A).

## 3 Models using analytic descriptions of the dark matter distribution

### 3.1 The lognormal model

We have started with the model introduced by Bi and collaborators (Bi et al. 1992, 1995; Bi 1993; BD97) for generating a Ly$`\alpha `$ absorption spectrum along a LOS. This simple model predicts many properties of the absorption lines, including the column density distribution and the distribution of line widths ($`b`$ parameters), which can be directly compared with observations (BD97). Recently, the BD97 model has been used by Roy Choudhury et al. (2000, 2001) to study neutral hydrogen correlation functions along and transverse to the line-of-sight. Feng & Fang (2000) adopted the BD97 method to analyse non-Gaussian effects on the transmitted flux stressing their importance for the reconstruction of the initial mass density field. Viel et al. (2002) implemented a variant of the BD97 model to simulate multiple systems of QSOs and found correlations in the transverse direction in agreement with observations (see also Petry et al. 2002 for an analysis of the correlations in the transverse direction).

The BD97 model is based on the assumption that the low-column density Ly$`\alpha `$ forest is produced by smooth fluctuations in the intergalactic medium which arise as a result of gravitational growth of perturbations. Since the fluctuations are only mildly non-linear, density perturbations in the intergalactic medium $`\delta _0^{\mathrm{IGM}}(𝐱,z)`$ can be related to the underlying DM perturbations by a convolution, which models the effects of gas pressure. In Fourier space one has:

$$\delta _0^{\mathrm{IGM}}(𝐤,z)=W_{\mathrm{IGM}}(k,z)D_+(z)\delta _0^{\mathrm{DM}}(𝐤)$$
(1)
where $`D_+(z)`$ is the growing mode of density perturbations (normalised so that $`D_+(0)=1`$) and $`\delta _0^{\mathrm{DM}}(𝐤)`$ is the Fourier transformed DM linear over density at $`z=0`$. The low-pass filter $`W_{\mathrm{IGM}}(k,z)=(1+k^2/k_J^2)^1`$ depends on the comoving Jeans length

$$k_J^1(z)H_0^1\left[\frac{2\gamma k_BT_m(z)}{3\mu m_p\mathrm{\Omega }_{0m}(1+z)}\right]^{1/2},$$
(2)
where $`k_B`$ is Boltzmann’s constant, $`T_m`$ the gas temperature, $`\mu `$ the molecular weight and $`\gamma `$ the ratio of specific heats. $`\mathrm{\Omega }_{0m}`$ is the present-day matter density. The Jeans length depends on density and temperature. Jeans smoothing is therefore an adaptive smoothing of the density field. However, for simplicity in most practical implementations the comoving Jeans length is assumed to be constant and computed for the mean density and the mean temperature $`T_0`$ at mean density $`\delta =0`$. As we will see later this is a severe restriction which significantly affects the results.

Gnedin & Hui (1998) adopt a different and more accurate expression for the IGM filter $`W_{\mathrm{IGM}}(k,z)`$, which, however, does not allow a simple matching with the non-linear regime (see also the discussion in Section 4.1). More accurate window-functions have also been proposed by Nusser (2000) and Matarrese & Mohayaee (2002). In what follows, we take $`T_0(z)1+z`$, which leads to a constant comoving Jeans scale. This assumption should not be critical as the redshift intervals considered here are small.

Bi & Davidsen adopt a simple lognormal (LN) transformation (Coles & Jones 1991) to obtain the IGM density in the mildly non-linear regime from the linear density field,

$$1+\delta _{\mathrm{IGM}}(𝐱,z)=\mathrm{exp}\left[\delta _0^{\mathrm{IGM}}(𝐱,z)\frac{(\delta _0^{\mathrm{IGM}})^2D_+^2(z)}{2}\right]$$
(3)
where $`1+\delta _{IGM}(𝐱,z)=n_{\mathrm{IGM}}(𝐱,z)/\overline{n}_{\mathrm{IGM}}(z)`$ and $`\overline{n}_{\mathrm{IGM}}(z)1.12\times 10^5\mathrm{\Omega }_{0b}h^2(1+z)^3`$ cm-3. The IGM peculiar velocity $`𝐯^{\mathrm{IGM}}`$ is related to the linear IGM density contrast via the continuity equation. As in BD97, we assume that the peculiar velocity is still linear even on scales where the density contrast gets non-linear; this yields

$$𝐯^{\mathrm{IGM}}(𝐤,z)=E_+(z)\frac{i𝐤}{k^2}W_{\mathrm{IGM}}(k,z)\delta _0^{\mathrm{DM}}(𝐤)$$
(4)
with $`E_+(z)=H(z)f(\mathrm{\Omega }_m,\mathrm{\Omega }_\mathrm{\Lambda })D_+(z)/(1+z)`$. Here $`f(\mathrm{\Omega }_m,\mathrm{\Omega }_\mathrm{\Lambda })d\mathrm{ln}D_+(z)/d\mathrm{ln}(1+z)`$ (e.g. Lahav et al. 1991, for its explicit and general expression) and $`H(z)`$ is the Hubble parameter at redshift $`z`$,

$$H(z)=H_0\sqrt{\mathrm{\Omega }_{0m}(1+z)^3+\mathrm{\Omega }_0(1+z)^2+\mathrm{\Omega }_{0\mathrm{\Lambda }}}$$
(5)
where $`\mathrm{\Omega }_{0\mathrm{\Lambda }}`$ is the vacuum-energy contribution to the cosmic density and $`\mathrm{\Omega }_0=1\mathrm{\Omega }_{0m}\mathrm{\Omega }_{0\mathrm{\Lambda }}`$ ($`\mathrm{\Omega }_0=0`$ for a flat universe). The procedure to obtain 1D LOS random fields from 3D random fields is described in BD97 and Viel et al. (2002). We will explore below a more accurate mapping between the linear and non-linear gas density which uses the PDF of the gas density obtained from hydro-dynamical simulations.

The neutral hydrogen density can be computed from the total density, assuming photo-ionisation and taking the optically thin limit, $`n_{\mathrm{HI}}(𝐱,z)=f_{\mathrm{HI}}(T,\mathrm{\Gamma },Y)n_\mathrm{H}(𝐱,z)`$, where $`n_\mathrm{H}=(1Y)n_{\mathrm{IGM}}`$ and $`Y`$ is the helium abundance by mass. The photo-ionisation rate $`\mathrm{\Gamma }\mathrm{\Gamma }_{12}\times 10^{12}`$s-1 is related to the spectrum of ionising photons as $`\mathrm{\Gamma }=_{\nu _{\mathrm{th}}}^{\mathrm{}}4\pi J(\nu )/h_p\nu \sigma (\nu )𝑑\nu `$, where $`\sigma (\nu )`$ is the photo-ionisation cross-section and $`h_P\nu _{\mathrm{th}}`$ the hydrogen ionisation threshold ($`h_P`$ denotes Planck’s constant). $`J(\nu )=J_{21}(\nu _0/\nu )^m\times 10^{21}\mathrm{erg}\mathrm{s}^1\mathrm{Hz}^1\mathrm{cm}^2\mathrm{sr}^1`$ with $`\nu _{\mathrm{th}}`$ the frequency of the HI ionisation threshold, and $`m`$ is usually assumed to lie between 1.5 and 1.8. In the highly ionised case ($`n_{\mathrm{HI}}n_{\mathrm{IGM}}`$) of interest here, one can approximate the local density of neutral hydrogen as (e.g. Hui, Gnedin & Zhang 1997)

$`{\displaystyle \frac{n_{\mathrm{HI}}(𝐱,z)}{\overline{n}_{\mathrm{IGM}}(z)}}`$ $``$ $`10^5\left({\displaystyle \frac{\mathrm{\Omega }_{0b}h^2}{0.019}}\right)\left({\displaystyle \frac{\mathrm{\Gamma }_{12}}{0.5}}\right)^1\left({\displaystyle \frac{1+z}{4}}\right)^3`$ (6)
$`\times \left({\displaystyle \frac{T(𝐱,z)}{10^4\mathrm{K}}}\right)^{0.7}\left(1+\delta _{\mathrm{IGM}}(𝐱,z)\right)^2.`$
The temperature of the low-density IGM is determined by the balance between adiabatic cooling and photo-heating by the UV background, which establishes a local power-law relation between temperature and density, $`T(𝐱,z)=T_0(z)(1+\delta ^{\mathrm{IGM}}(𝐱,z))^{\gamma (z)1}`$, where both the temperature at mean density $`T_0`$ and the adiabatic index $`\gamma `$ depend on the IGM ionisation history (Meiksin 1994; Miralda-Escudé & Rees 1994; Hui & Gnedin 1997; Schaye et al. 2000; Theuns et al. 2001).

Given the neutral density, the optical depth in redshift-space at velocity $`u`$ (in km s-1) is

$$\tau (u)=\frac{\sigma _{0,\alpha }c}{H(z)}_{\mathrm{}}^{\mathrm{}}𝑑yn_{\mathrm{HI}}(y)𝒱[uyv_{}^{\mathrm{IGM}}(y),b(y)]$$
(7)
where $`\sigma _{0,\alpha }=4.45\times 10^{18}`$ cm2 is the hydrogen Ly$`\alpha `$ cross-section, $`y`$ is the real-space coordinate (in km s-1), $`𝒱`$ is the standard Voigt profile normalised in real-space, and $`b=(2k_BT/mc^2)^{1/2}`$ is the thermal width. Velocity $`v`$ and redshift $`z`$ are related through $`d\lambda /\lambda =dv/c`$, where $`\lambda =\lambda _0(1+z)`$. For the low column-density systems considered here, the Voigt profile is well approximated by a Gaussian: $`𝒱=(\sqrt{\pi }b)^1\mathrm{exp}[(uyv_{}^{\mathrm{IGM}}(y))^2/b^2]`$. As stressed by BD97 peculiar velocities affect the optical depth in two different ways: the lines are shifted to a slightly different location and their profiles are altered by velocity gradients. In our modelling, we treat $`\mathrm{\Gamma }_{12}`$ as a free parameter, which is varied the observed effective opacity $`\tau _{\mathrm{eff}}(z)=\mathrm{ln}\mathrm{exp}(\tau )`$ is matched (e.g. McDonald et al. 1999; Efstathiou et al. 2000) at the median redshift of the considered range ($`\tau _{\mathrm{eff}}=0.12`$ and $`\tau _{\mathrm{eff}}=0.27`$ at $`z=2.15`$ and $`z=3`$, respectively, in our case). We determine $`\mathrm{\Gamma }_{12}`$ by requiring that the ensemble averaged effective optical depth is equal to the observed effective optical depth. The transmitted flux is then simply $`=\mathrm{exp}(\tau )`$.

### 3.2 Improving the mapping from the linear dark matter density to the non-linear gas density

In this subsection we describe how the mapping from the linear DM density to the non-linear gas density can be improved. As described in the previous section the semi-analytical modelling of Ly$`\alpha `$ forest spectra involves two main steps:

* smoothing of the linear DM density field to obtain a linear gas density field
* a local mapping from the linear to the non-linear gas density.

In this section we will substitute the log-normal mapping used by BD97 for the second step by a rank-ordered mapping from the linear gas density field to the PDF of the gas density in the hydro-simulation. Figure 1 compares the PDF of the DM density field obtained using SPH interpolation and that of the gas density field. Pressure forces push the gas out into the surrounding voids (‘Jeans smoothing’) and for this reason the PDF of the gas drops below that of the DM density at low densities. At high density the gas is converted into stars and the PDF of the gas density drops again below the PDF of the DM density. As a result the PDF of the gas density is more peaked (Theuns, Schaye & Haehnelt 2000). Figure 2 compares the resulting probability distribution of the flux for simulated spectra using the lognormal model our ‘improved’ model (hereafter PDF model) and the numerical simulations. For this comparison we have produced simulated spectra for the lognormal model and the PDF model with the same cosmological parameters as in our hydrodynamical simulation. We further imposed a power law temperature-density relation with $`T_0=10^{4.3}`$ K and $`\gamma =1.2`$, which fits the simulation well (Figure 3). The simulated spectra have larger length ($`16000\mathrm{km}\mathrm{s}^1`$) than the spectra extracted from hydro-simulations ($`1400\mathrm{km}\mathrm{s}^1`$), but are obtained with the same spectral resolution, approximately of $`2\mathrm{km}\mathrm{s}^1`$. We take into account that the hydro-simulations is missing large scale power by applying a similar cut-off in the power spectrum used to calculate the semi-analytical spectra. This should allow a fair comparison. The peculiar velocity field is again assumed to be given by linear theory (Eq. 4). All spectra are scaled to $`\tau _{eff}=0.27`$ at $`z=3`$ and no noise is added.

The probability distribution of the flux for our PDF model agrees significantly better with that of the spectra obtained from our hydro-dynamical simulation. There are still small but significant differences. While the agreement for this 1 point statistic seems reasonable the agreement becomes very poor for some 2-point statistics of the flux. The function $`P(F_1,F_2,\mathrm{\Delta }v)dF_1dF_2`$ is the probability that two pixels separated in the spectrum by a velocity difference $`\mathrm{\Delta }v`$, have transmitted fluxes in intervals $`dF_1`$, $`dF_2`$ around $`F_1`$ and $`F_2`$, respectively. In the right panel of Figure 2 we show the mean flux difference,

$$\mathrm{\Delta }F(F_1,\delta v)=P(F_1,F_2,\delta v)(F_1F_2)𝑑F_2,$$
(8)
for 2 values of $`F_1`$, as a function of $`\mathrm{\Delta }v`$ (see Miralda-Escudé et al. 1997, Theuns et al. 2000, Kim et al. 2001). The results for the lognormal model and for the PDF model are almost identical, which is not surprising since we checked that the two functions that map the linear density field into the non-linear density field are similar.

For strong absorbers ($`0F0.1`$) the agreement between the numerical simulation and our models which rely on an analytical description of the DM density field is very poor. What could be the reason for this discrepancy? We have made the following major simplifications: (i) the peculiar velocity field is assumed to be linear; (ii) we do not simulate the intrinsic scatter present in the relation $`T_{IGM}`$ vs. $`\delta _{IGM}`$; (iii) we recover directly the neutral hydrogen density along the LOS by assuming eq. (6), which is valid only under some assumptions (see Subsection 3.1); (iv) we input only the 1-point probability distribution function of the IGM, without taking into account the higher order moments of the distribution; (v) the correlations are assumed to be those predicted by linear theory, modified by the non-linear mapping, while the correlations of the hydro-simulation are different.

Concerning point $`(i)`$ it has been shown by Hui et al. (1997) that the peculiar velocity field can affect the shape of the absorption features, while it has very small effect on the column density distribution functions of the lines. We have run the calculated spectra varying the peculiar velocity field and have not found a strong dependence of the PDF of the flux on the peculiar velocity quantity. This means that this statistics is mainly influenced by the underlying density field. We have compared the PDF of the peculiar velocity field predicted by linear theory with the peculiar velocity found in the $`z=3`$ output of hydro-simulation for the IGM and we have found that the differences are not big. The assumption of a linear velocity field should therefore be a good approximation.

If we try to simulate the intrinsic scatter in the relation $`T\delta `$, left panel of Figure 3, for example by using as input a higher order polynomial fit instead of the power-law relation, the improvement in the PDF of the flux is negligible. Provided we set the same $`T_0`$ of the simulations and a ‘reasonable’ value of $`\gamma `$, the power-law equation of state is a good approximation and the scatter can be neglected.

The same considerations are valid for point $`(iii)`$. In fact, we have found that eq. (6) is in good agreement with hydro-simulations (Figure 3, right panel), showing that the approximations we made (optical thin limit, high ionisation state) are not affecting the results significantly.

We thus believe that point $`(iv)`$ and $`(v)`$ are the main simplification involved in the model. This is somewhat expected as in the definition of the optical depth, eq. (7), the summation is done over nearby pixels, so it is influenced by the spatial correlations, which are assumed to be those predicted by the linear theory, modified by the local non-linear mapping. In the following Sections we will thus turn to models that are based on the DM density field obtained directly from a $`N`$-body simulation (Fig. 13 below). These models reproduce the flux correlation for strong absorbers much better. This demonstrates that correlations in the density field introduced by the non-linear evolution are indeed responsible for the shape of strong lines. These are not not reproduced accurately by our local mapping from the linear to non-linear gas density. The marked differences for strong absorbers are not unexpected as the $`z=3`$ output of hydro-simulation contains very strong absorption systems from fully collapsed objects which are not related to the Ly$`\alpha `$ forest. The spectra of the LOGN model have been produced using the Jeans length at the mean temperature and with a fixed value of $`\gamma `$. If we use other values for the temperature $`T_0`$ and for $`\gamma `$ we can obtain a better agreement in other flux intervals. In Figure 2 the range of flux values in which the agreement with the flux PDF of the hydro-dynamical simulation and in $`\mathrm{\Delta }F`$ is good corresponds to intermediate absorbers ($`0.2<F<0.4`$), this is in part determined by the requirement of having a fixed $`\tau _{eff}=0.27`$ for the ensemble of simulated spectra in each model, which scales the neutral hydrogen fraction.

## 4 Models using Numerical simulations of the dark matter distribution

### 4.1 Zel’dovich modelling of the gas distribution

In this Section we present models of the gas distribution based on a modified filtering of the initial conditions of the dark matter density field. The first method is based on the truncated Zel’dovich approximation (TZA). Among the possible approximations, TZA has been found to produce the best agreement with N-body results (see Coles et al. 1993, Melott et al. 1994, Sathyaprakash et al. 1995). Hui et al. (1997, hereafter HGZ) showed that the TZA, with an appropriate smoothing and with a recipe which allows to convert density peaks into absorption lines, successfully reproduces the observed column density distribution of the Ly$`\alpha `$ forest lines over a wide range of $`N_{HI}`$. Gnedin & Hui (1998) presented a more accurate semi-analytical model of the Ly$`\alpha `$ forest by combining a particle mesh solver modified to compute also an effective potential due to gas pressure. They showed that a particle mesh solver, with an appropriate filtering of the initial condition, can be used to model the low density IGM. All these different methods have the advantage of being much faster than a hydro-simulation.

We are now going to make a comparison, LOS by LOS, between the simulated gas distribution and the effective gas distribution of the $`z=3`$ output of the hydro-simulation. Our simulation is significantly larger and with higher resolution than those used in previous studies.

Following Hui et al. (1997), we define a non-linear wavenumber $`k_{nl}`$ (see also Melott et al. 1994):

$$D_+^2(z)_0^{k_{nl}}P(k)d^3𝐤=1,$$
(9)
where $`D_+(z)`$ is the linear growth factor for the density perturbations and $`P(k)`$ the linear power-spectrum. We filter the linear density field with a Gaussian window $`W(k,k_s)=\mathrm{exp}(k^2/2k_s^2)`$, with $`k_s1.5k_{nl}`$. At this point we displace particles from the initial Lagrangian coordinates $`𝐪`$ according to the truncated Zel’dovich approximation:

$$𝐱(𝐪,z)=𝐪+D_+(z)_𝐪\varphi _f(𝐪),$$
(10)
with $`\varphi _f(q)`$ the initial filtered velocity potential (see e.g. Coles et al. 1995), with $`\varphi (𝐤)`$ the actual potential used in the initial conditions of the simulation. The filtering is intended to prevent shell-crossing. To mimic baryonic pressure, we smooth the initial density field with a Gaussian window $`\mathrm{exp}(k^2/2k_J^2)`$, with $`k_J`$ given by eq. (2), as in HGZ. The final filtering scale is effectively the smaller of $`k_s`$ and $`k_J`$. In our case, $`k_s5.5\mathrm{Mpc}^1`$ and $`k_J7\mathrm{M}\mathrm{p}\mathrm{c}^1`$, these wave numbers correspond to scales of $`\lambda _s1.14\mathrm{Mpc}`$ and $`\lambda _J0.9\mathrm{Mpc}`$, respectively. This means that the amount of filtering, before displacing particles, is given by the condition of eq. (9).

We also outline here a further approach based on the assumption that the displacement between DM and IGM particles is ‘Zel’dovich-like’. We refer to this method as the Zel’dovich Displacement (ZD method), while the method in which we assume a perfect tracing between the DM and IGM density field will be referred as the DM method. The difference, at redshift $`z`$, between the Eulerian position of a gas and a DM particle which have the same Lagrangian coordinate $`𝐪`$, according to the Zel’dovich approximation, is:

$$\mathrm{\Delta }𝐱(𝐪,z)=D_+(z)[_𝐪\psi _{IGM}(𝐪,z)_𝐪\varphi _{DM}(𝐪)],$$
(11)
where $`\psi _{IGM}`$ and $`\varphi _{DM}`$ are the IGM and DM velocity potentials. In Fourier space we have (Matarrese & Mohayaee 2002) $`\psi _{IGM}(𝐤,z)=W_{IGM}(𝐤,z)\varphi _{DM}(𝐤)`$ and from eq. (11) we get:

$$\mathrm{\Delta }𝐱(𝐤,z)=D_+(z)[W_{IGM}1]i𝐤\varphi _{DM}(𝐤).$$
(12)
The above equation shows that for $`kk_J`$ $`W_{IGM}1`$, $`\mathrm{\Delta }𝐱(𝐤,z)0`$, while for $`kk_J`$, $`W_{IGM}0`$, $`\mathrm{\Delta }𝐱(𝐤,z)i𝐤\varphi _{DM}(𝐤)`$. In the first limit, Jeans smoothing is not effective and the IGM traces the DM, in the second limit the effect of gas pressure prevents large displacements of baryons from their initial location. According to this method, the displacement between DM and IGM particles with the same Lagrangian coordinate is given by a filtering of the initial DM density field. We calculate this displacement at $`z=3`$ with the filter proposed by Gnedin & Hui (1998), i.e. instead of $`W_{IGM}`$ we use $`W_{ZD}=\mathrm{exp}(k^2/2k_{ZD}^2)`$, with $`k_{ZD}2.2k_J`$. This choice is motivated by the fact that a Gaussian filter gives an excellent fit to baryon fluctuations for a wide range of wave numbers. The filtering scale $`k_{ZD}`$ is the value predicted by linear theory at $`z=3`$, assuming that reionization takes places at $`z7`$ (Gnedin & Hui 1998). The location of the gas particle is then found be adding this displacement to the actual position at $`z=3`$ of the DM particle with the same Lagrangian coordinate.

In Figure 4 we show four panels, which are slices of thickness $`0.1h^1`$ comoving Mpc at the same position along the z-axis, for the IGM and the DM distribution of the hydro-simulations (top panels), the TZA field (bottom left panel) and the IGM field obtained with ZD method (bottom right).

By comparing the top panels, one can see that the gas is more diffuse than the dark matter. TZA reproduces the main filaments but we know that the agreement will be better in the low-density regions. ZD (with a filtering at $`k_{ZD}16\mathrm{M}\mathrm{p}\mathrm{c}^1`$) seems promising, at least by eye. This method allows some diffusion around the dark matter to reproduce the gas distribution. A more quantitative comparison is needed to see how good our IGM density field is compared to the hydro one. We compare LOS by LOS the SPH interpolated IGM density fields with the ‘true’ IGM field of the hydro-simulation. This test is pretty severe as we are not going to filter the density fields and the comparison is made ‘pixel-by-pixel’.

The results are shown in Figure 5 where we compare the amount of scatter predicted by the TZA approximation (left panel) and the ZD method (right panel). As expected, the test without any filtering on the final density field shows that the TZA tracks the simulation, but the scatter is very large. The ZD method agrees better with the simulation, but the scatter is still rather significant. A more detailed comparison is shown in Figure 6, where we plot the mean, and scatter around the mean, as a function of density.

Clearly ZD works significantly better than TZA, the average value is better reproduced and also the scatter is significantly smaller. The inadequacy of the TZA was also shown by Bond & Wadsley (1997). The amount of scatter and the average value found with the ZD method are consistent with the results of Gnedin & Hui (1998; see their Figure 4).

### 4.2 Using the mean $`\delta _{DM}\delta _{IGM}`$ relation of the hydro simulation to predict the gas distribution from the DM distribution.

The results of the previous section have demonstrated that ‘filtering’ techniques on the initial conditions do not result in density fields that agree well with hydro-simulations, at least for a ‘pixel-by-pixel’ comparison with a high-resolution simulation.

In this section we present an alternative approach that starts from the actual DM density distribution of the numerical simulation and “predicts” the gas distribution using a fit to the mean relation between gas and DM density. In this way the displacement of the gas with respect to the dark matter is modelled statistically in real space. The techniques of the previous Section are all based on different filtering schemes which smooth the IGM density field over a constant scale set by the Jeans length at the mean density. This filtering is done in Fourier space and leads to equal smoothing of all dark matter structures independent of their density. As discussed above this strong simplification is responsible for the rather strong discrepancies with the gas distribution in our hydro-dynamical simulation.

We compute IGM and DM over densities, IGM peculiar velocity and IGM temperature obtained with SPH interpolation (see Appendix A for details). From the computed IGM and DM density profiles along different LOS we argue that a simple relation between the DM and the gas IGM density fields does not exist (see also Cen et al. 1994, Croft et al. 1998; Zhang et al. 1998). For a wide range of moderate over densities, i.e. $`0.75\mathit{}\mathit{<}\text{ }\delta \mathit{}\mathit{<}\text{ }5`$, the PDF of the gas distribution lies above that of the DM (Fig. 1): this of course does not determine the relation between gas and dark matter, but shows that we should expect a significant number of regions in which $`\delta _{IGM}`$ is larger than $`\delta _{DM}`$.

This is what we find in the hydro-simulations in which there are regions at high $`\delta `$ where the gas is more concentrated than the dark matter. At low $`\delta `$ there is an almost perfect agreement between the gas and the dark matter distribution, while at moderate and large over densities the relation between gas and dark matter is not unique. The highest density peaks are usually related to peaks in the temperature of the gas of and strong gradients in its peculiar velocity, suggesting the presence of a shock. However, it can happen that a region shows a gas density more peaked than the DM, while the velocity field and the temperature show no particular behaviour. This means that modelling of the gas distribution using only the information along the LOS will not be very accurate on a point-to-point basis, but we might nevertheless hope to obtain a model for the gas distribution which has the correct properties in a statistical sense.

We start by plotting the values of $`\delta _{IGM}`$ vs. $`\delta _{DM}`$ obtained from the simulation at different redshifts. Figure 7 demonstrates that the scatter between the IGM and DM densities increases with decreasing redshift, as expected. At $`z=10`$, there is almost perfect agreement between the two fields, but when cosmic structures get non-linear the physics becomes more complex. In addition, reionization occurs in our simulation at $`z6`$, so before that the simulations do not resolve the then much smaller Jeans length. At later redshifts, the effect of the baryonic pressure can be seen. For under dense regions, the IGM density tends to be higher than the DM density, because pressure is pushing gas into the low density voids. For larger values of $`\delta _{DM}`$ the scatter increases and gas can be either denser or less concentrated than the dark matter, depending on star formation which sets-in in the simulations at over densities $`80`$ (see Aguirre, Schaye and Theuns 2002 for a more detailed description of the star formation recipe adopted).

We fit the relation between gas and DM density with a 3rd order polynomial; the results are shown in Figure 8, where the range of $`\delta _{DM}`$ plotted is the one of interest for the Ly$`\alpha `$ forest, $`1<\text{ }\delta <\text{ }6`$. If we set $`y=\mathrm{log}(1+\delta _{IGM})`$ and $`x=\mathrm{log}(1+\delta _{DM})`$ the fitting function is:

$`y`$ $`=`$ $`0.02^{\pm 0.03}+0.91^{\pm 0.05}x0.15^{\pm 0.05}x^2`$ (13)
$``$ $`0.08^{\pm 0.04}x^3,1.5x1.0,`$
where the errors reported are the $`1\sigma `$ uncertainties of each coefficient. The same function has been found to be a reasonable good fit at $`z=4`$ and $`z=2`$ as well, so it can be used to simulate the gas distribution in this redshift range for a $`\mathrm{\Lambda }`$CDM model. From Figure 8 one can see that the fitting function predicts that the gas is on average more concentrated than the dark matter for a wide range of $`\delta _{DM}`$, although significant scatter is present.

Is it possible to reduce this scatter using other information?

In Figure 9 (left panel) we plot the peculiar velocity gradient of the dark matter density fields vs. the difference $`\mathrm{\Delta }`$ between the ‘true’ value of $`\mathrm{log}(1+\delta _{IGM})`$ and the ‘fit’ value obtained with eq. (13). Both these quantities have been smoothed over a scale of $`100\mathrm{km}\mathrm{s}^1`$ to reduce noise. $`\mathrm{\Delta }`$ is weakly anti-correlated with the dark matter peculiar velocity gradient, with regions with a negative gradient lying above the mean fit and hence have $`\mathrm{\Delta }>0`$. On average, a negative velocity gradient indicates that the gas is being compressed and so may be undergoing moderate or strong shocks. As a consequence, the gas is also being heated, and this introduces the correlation between $`\mathrm{\Delta }`$ and the temperature of the gas (Fig. 9, right panel). This correlation is stronger, because the gas temperature is a more direct indicator of a shock. Conversely, negative $`\mathrm{\Delta }`$, i.e. points below the fitting function, are related to colder regions of gas and small positive gradients of the dark matter peculiar velocity field. As expected, the bulk of the points is in quiet regions with $`dv_{pec}/dx0`$ and temperatures between $`10^{4.2}`$ K and $`10^{4.3}`$ K.

These results suggest that we can indeed model the gas distribution that corresponds to a given dark matter density and peculiar velocity distribution statistically. We use the fit of eq. (13) to predict the gas density field from the dark matter density field. The modelling can be further improved by taking the peculiar velocity of the dark matter into account. For this purpose we fit the relation shown in the left panel of Figure 9 with a 2nd order polynomial because the correlation is weak. Models obtained using the density-density fit and models using a two parameter density-density plus the density-peculiar velocity fit will be referred to as FIT-1 or FIT-2, respectively.

The fitting function of eq. (13) gives a good approximation to the gas distribution for a $`\mathrm{\Lambda }`$CDM model in the redshift range $`2<\text{ }z<\text{ }4`$. The same technique can be applied to other cosmological models and different redshift ranges. We will now make a LOS by LOS comparison for the the gas density distribution of the two models FIT-1 and FIT-2 to that of our hydro-dynamical simulation. In Figure 11 we show the scatter plots for the gas (over)densities for the two models. Figure 10 shows the scatter in terms of the standard and mean deviation for both models. The simple DM model where the gas is assumed to trace the dark matter exactly is also shown (thick, thin and dotted lines, respectively).

Figure 11 illustrates the improvement of using the fitting procedure of model FIT-2. The $`\chi ^2`$ for these fits, where $`\chi ^2N_{bin}^1_{i=1}^{N_{bin}}\frac{(\delta _i\delta )^2}{\sigma _i^2}`$ ($`N_{bin}=25`$ is the number of bins) are 0.24, 0.21 and 0.19 for the FIT-1, FIT-2 and DM, respectively, showing that the scatter is effectively reduced with the second method but still in the total interval $`1.5<\mathrm{log}(1+\delta _{IGM})<1`$.

When the density range is constrained to be in the interval $`1<\mathrm{log}(1+\delta _{IGM})<0.6`$, relevant for the Lyman$`\alpha `$ forest, we obtain $`\chi ^2`$ of 0.06, 0.03 and 0.10 for the FIT-1, FIT-2 and DM model respectively, illustrating the improvement of the fitting procedure. Comparing with Figure 6, we see that the FIT-2 method is significantly better than ZD in reproducing the mean values for $`\mathrm{log}(1+\delta _{IGM})<\text{ }0.5`$ (however, for larger over densities ZD shows a better agreement than FIT-2 or FIT-1), while the rms values are basically equivalent. If we apply the fitting procedure in the range $`\delta _{IGM}<\text{ }1`$ which represents the bulk of the IGM, the mean and rms deviation between the fitted and true IGM densities are smaller than 10 and 30 per cent, respectively.

### 4.3 Jeans smoothing of the evolved DM density field

Here we investigate the effect of smoothing the evolved rather than the initial DM density field on a constant scale set by the Jeans length at mean density. We use the DM distribution of the $`z=3`$ output of the hydro-simulation to compute the dark matter density field on a cube of $`128^3`$ mesh points using SPH interpolation as described in the Appendix. The cell size of this mesh of $`0.14`$ co-moving Mpc approximately resolves the Jeans length $`\lambda _J1`$Mpc. We have to choose a mesh because the filtering is done in 3D using fast Fourier transforms. However, we note that also 1D smoothing along LOS is in rough agreement with the 3D one. To model the Jeans smoothing we convolve the DM density field with a Gaussian filter $`W=\mathrm{exp}(k^2/2k_f^2)`$, with $`k_f=k_J7\mathrm{M}\mathrm{p}\mathrm{c}^1`$. By comparing Figure 8 with Figure 12 one can see that the Jeans smoothed dark matter density field is very different from IGM distribution in the hydro-simulation. Given this large discrepancy between hydro-simulations and Jeans smoothed dark-matter density field we choose not to analyse into the details the density statistics as we did for the models of the previous subsection.

The main reason for the large difference between the true IGM density distribution, and the one obtained from Jeans smoothing the DM density field, is the simplification of a constant smoothing length. In reality, the Jeans length depends on temperature, and therefore should be adaptive. Unfortunately, Fourier space filtering techniques smooth all structures in the same way, since they are global operations. In low-density regions, the amount of smoothing in hydrodynamical simulations is small, but the Jeans smoothed density field is very different from the original DM field. On the other hand, for large over densities of the dark matter, Jeans filtering at the mean gas density underestimates the amount of smoothing, the region at $`\delta _{IGM}<\delta _{DM}`$, for large values of $`\delta _{DM}`$, is more populated in Figure 12 than in Figure 8. As we will see in the next section the resulting flux statistics are very different as well.

## 5 A comparison of flux statistics for the improved methods

In this section we compare the flux distribution obtained with the different methods for modelling the gas density field. Figure 13 shows the 1-point probability distribution function for the flux (left panel) and the mean flux difference $`\mathrm{\Delta }F`$ (Eq.8, right panel), for three models: (i) DM (dotted line) is the PDF we get by setting $`\delta _{IGM}=\delta _{DM}`$; (ii) ZD (dashed lines) is the PDF obtained with the ‘Zel’dovich Displacement’ method of Section 4.1 with a $`k_f16\mathrm{M}\mathrm{p}\mathrm{c}^1`$; (iii) FIT-2 (continuous line) is the PDF obtained with the fitting technique of the previous Section and adding the scatter inferred from the dark matter peculiar velocity field. We also show the flux PDF of the hydro-dynamical simulation (triangles) which is the true PDF extracted from the $`z=3`$ output of the $`\mathrm{\Lambda }`$CDM model.

All these methods try to predict the gas distribution from the dark matter distribution of the numerical simulation. We thereby assume that $`v_{HI}v_{IGM}v_{DM}`$ which we have checked to be a good approximation. The local abundance of neutral hydrogen is computed from the ionisation equilibrium equation, eq. (6). We further need to assume the temperature at mean density, $`T_0`$, and the slope of the temperature-density relation, $`\gamma `$. We take $`T_0=10^{4.3}`$K and $`\gamma =1.2`$. All simulated spectra have been scaled to the same effective optical depth, $`\tau _{eff}0.27`$ (which is a good fit to the observed effective optical depth at redshift $`z=3`$).

The PDF of the flux is very similar for the three models. We performed a quantitative comparison based on the Kolmogorov-Smirnov (KS) test which characterises the difference between two models from the maximum absolute deviation, $`d_{\mathrm{KS}}`$, between the two cumulative flux distributions (see Meiksin et al. 2001 for further details on KS test applied to the PDF of the flux). We calculate $`d_{KS}`$ for the lognormal model (LOGN) and the ‘improved’ model (PDF) obtained by implementing the 1-point probability distribution function of the IGM as well. Fluxes obtained with TZA have not been computed as the large scatter found (left panel, Figure 5) suggests that this approximation is the least accurate. The KS test was performed for a total of $`310^5`$ pixels (Table LABEL:tab:KS).

The best agreement with the hydro-simulation is obtained with the FIT-2 method. The KS indicates a slightly better agreement than for the FIT-1 method. The ZD method gives better agreement than assuming that the gas traces the dark matter. The PDF and LOGN methods are the least accurate; this means that the proposed fitting methods determine a better agreement in terms of the pdf of the flux and can be considered as an effective improvement compared to the PDF and LOGN methods.

The methods based on the actual DM density distribution (FIT-1, FIT-2, DM) also reproduce the 2-point flux distribution of the hydro-simulation dramatically better than those based on linear theory or the lognormal model as can be seen by comparing Fig. 13 (right panel) with Fig. 2. Note especially the significant improvement in the shape of strong lines. As discussed above this implies that this statistic is mainly influenced by the correlations in the underlying dark matter density fields which are now the same for all the models.

We choose to make a final plot to quantify the differences between the Jeans smoothed dark matter density field and the hydro-simulations in terms of the PDF of the flux. Results are shown in Figure 14. The number of pixels per LOS is 128, significantly lower than in the comparisons between the other improved methods 4.2. To produce ‘Jeans smoothed’ spectra we assume again that $`v_{HI}v_{DM}`$ and a power-law equation of state for the dark matter over density. This time there are also significant differences between the the two models also in terms of the PDF of the flux. The hydro PDF predicts more dense regions than the Jeans smoothed one, and correspondingly less regions with small densities. This is somewhat counterintuitive, since the fitting technique smooths less in low density regions compared with Jeans smoothing. However, we require both types of spectra to have the same $`\tau _{eff}`$. Effectively, this means that the two PDFs for $`\delta `$ are multiplied by a correction factor, which can be seen as a scaling of the neutral hydrogen, in such a way that the hydro PDF has a larger number of regions at large $`\delta `$ than the Jeans smoothed field. The $`d_{KS}`$ of the two distribution is $`0.1`$ slightly better than the LOGN model value.

The main conclusions of this Section are: (i) the use of the numerically simulated dark matter density distribution for the prediction of the gas density distribution results in a significant improvement in the comparison with the hydro-simulations; (ii) the modelling of the gas distribution using the fit to the mean relation between dark matter and gas density in the hydro-simulations results in much better agreement than the methods where a filtering scheme is applied to the initial or evolved density field.

## 6 Discussion and conclusions

Many aspects of the warm photo-ionised Intergalactic medium can be well modelled by hydrodynamical simulations. These are, however, still rather limited in dynamic range and lack the possibility of extensive parameter studies due to limited computational resources. In order to overcome these problems we have tested several approximate methods for simulating the Ly$`\alpha `$ forest in QSO absorption spectra. The modelling consists of two main steps (i) modelling the DM distribution and (ii) mapping the DM distribution into a gas distribution.

Methods which use an analytic description of the DM distribution like the lognormal model give a rather poor description of the gas distribution compared to numerical simulations. In Section 3 we have shown that this results in flux PDFs which are in rough agreement with the PDF extracted from hydro-simulations. The agreement can be improved with a variant of the lognormal model where a mapping between linear and non-linear IGM density fields calibrated by hydro-dynamical simulations is used. However, in both cases, the 2-point PDF of the flux differs significantly from that obtained from numerical simulations. This is because the 2-point PDF is strongly affected by the underlying correlations in the DM density field, which are not well reproduced in models based on an extrapolation of linear theory. To take properly account of these correlations on a point-to-point basis numerical DM simulations are required.

We thus tested a variety of schemes to relate the DM distribution of numerical simulations to the gas distribution. These schemes are supposed to take into account that the gas is smoothed on a Jeans scale relative to the dark matter. We first investigated two approximations based on the Zel’dovich approximation, the truncated Zel’dovich approximation (TZA) and a scheme which we call Zel’dovich displacement (ZD). The latter is based on the assumption that the displacement between DM and IGM at the same Lagrangian coordinate depends on the initial DM density field, filtered on a suitable scale. The TZA reproduces the gas density field very poorly in a LOS by LOS comparison. It is actually worse than if we assume that the gas traces the DM faithfully. The ZD method which allows diffusion on a scale smaller than half of the Jeans length to mimic baryonic pressure fares somewhat better. The scatter in plots of predicted vs. simulated densities is nevertheless only slightly smaller than in a model where gas traces dark matter. We have also tested the ZD method with smoothing on a global Jeans length and have again found poor agreement with spectra extracted from hydro-simulations. It may be that filtering techniques based on the Zel’dovich approximation can reproduce some statistical properties of the Ly$`\alpha `$ forest like the column density distribution reasonably well, but they fail in reproducing the flux distribution in detail.

To make progress we have thus investigated the relation of the gas density and the DM density in the hydro-dynamical simulation in more detail. The relation between $`\delta _{IGM}`$ and $`\delta _{DM}`$ can be well fit with a 3rd order polynomial. There is considerable scatter around the mean relation and we have examined the correlations between deviations from the fit with other physical quantities along the LOS. There is a weak correlation of these deviations with the filtered dark matter peculiar velocity gradient and a somewhat stronger correlations with the gas temperature. This indicates that the deviations are due to moderate or strong shocks in the gas component.

Combining the DM simulations with the fitted relation between DM density and gas density gives good results for both the one and two-point distribution of the flux. If we introduce an additional correlation of the gas density with the DM peculiar velocity gradient as found in the hydro simulations the agreement is further improved (a method which we called FIT-2). These fitting methods give significantly better results than the other methods we have discussed.

Smoothing of the dark matter density field with a constant global Jeans scale calculated for mean density and temperature results in a gas distribution very different from that found in hydro-simulations. This is not too astonishing as it does not take into account that the Jeans length depends on temperature and density. Jeans smoothing at the mean temperature overestimates the smoothing in low density regions and underestimates it at higher density. This leads to significant differences of the flux statistics compared to hydro-simulations. Smoothing of the $`\mathrm{𝑒𝑣𝑜𝑙𝑣𝑒𝑑}`$ dark matter density field on a global Jeans scale is therefore not a promising technique. The reason for the success of our FIT-1 and FIT-2 schemes is their adaptive nature which takes into account – at least to some extent – the density/temperature dependence of the Jeans scale.

We conclude that large high-resolution DM simulations combined with a two-parameter fit of the DM density gas density relation obtained from hydro-dynamical simulations are the best compromise between computational expense and accuracy when a large dynamic range and/or an extensive parameter study are required.

## Acknowledgments

We thank Bepi Tormen, Francesco Miniati, Lauro Moscardini, Joop Schaye and Simon White for useful discussions and technical help. MV acknowledges partial financial support from an EARA Marie Curie Fellowship under contract HPMT-CT-2000-00132. TT thanks PPARC for the award of an Advanced Fellowship. This research was conducted in collaboration with Cray/SGI, utilising the COSMOS super computer at the Department for Applied Mathematics and Theoretical Physics in Cambridge. This work was supported by the European Community Research and Training Network ‘The Physics of the Intergalactic Medium’.

## Appendix A SPH interpolation along LOS

In this Appendix we describe the SPH computation of physical quantities along a given line of sight through the box. We follow the same procedure described in (Theuns et al. 1998). We divide the sight line into $`N2^{10}`$ bins of width $`\mathrm{\Delta }`$ in distance $`x`$ along the sight line. For a bin $`j`$ at position $`x(j)`$ we compute the density and the density weighted temperature and velocity for the gas and density and weighted velocity for the dark matter from:

$`\rho _X(j)`$ $`=`$ $`{\displaystyle \underset{i}{}}𝒲_{ij}`$ (14)
$`(\rho v)_X(j)`$ $`=`$ $`{\displaystyle \underset{i}{}}v_{X(i)}𝒲_{ij}`$ (15)
$`(\rho T)_X(j)`$ $`=`$ $`{\displaystyle \underset{i}{}}T(i)𝒲_{ij}.`$ (16)
$`X(i)`$ is a label indicating the abundance of species $`X`$ of particle $`i`$ ($`X=\mathrm{H}\mathrm{i}`$; $`X=`$IGM and $`X=`$DM denotes neutral hydrogen, total gas and dark matter density, respectively).

Here, $`𝒲_{ij}=mW(q_{ij})/h_i^3`$ and $`m`$ is the SPH particle mass which is the same for all SPH particles (but different for DM and gas particles). For $`W`$ we use the M4 spline (Monaghan 1992) given by

$`W(q)`$ $`=`$ $`{\displaystyle \frac{1}{\pi }}(1.+q^2(1.5+0.75q))\mathrm{if}q1`$ (17)
$`=`$ $`{\displaystyle \frac{1}{\pi }}(0.25(2.q)^3)\mathrm{if}\mathrm{\hspace{0.17em}1}q2`$
$`=`$ $`0\mathrm{elsewhere}.`$
We have defined:

$$q_{ij}=\frac{|𝐱(i)𝐱(j)|}{h_i},$$
(18)
where $`𝐱(i)`$ and $`h_i`$ are the position and SPH-smoothing length of particle i. Note that $`h`$ is defined in such a way that on average 32 particles are within $`2h(i)`$ from particle i. In this way, for each pixel along the LOS, we compute the contribution of all the particles which influence this region with a weight given by eq. (17); this is the ‘scatter’ interpretation (see, for example, Hernquist & Katz 1989). For the computation of the spectra we label bins according to velocity and we adopt the procedure described in Section 3.


format: text

text:
# 1 Introduction

## 1 Introduction

Exclusive double-diffractive-like processes of the type

$$ppp+M+p$$
(1)
can significantly increase the physics potential of high energy proton colliders. Here $`M`$ represents a system of invariant mass $`M`$, and the + signs denote the presence of rapidity gaps which separate the system $`M`$ from the protons. Such processes allow, on the one hand, novel studies of QCD and of the diffractive amplitude at very high energies, while, on the other hand, allow an exceptionally clean experimental environment to identify New Physics signals (such as the Higgs boson, SUSY particles, etc., see  and ref. therein). Moreover tagging two forward protons offers an attractive extension of the proton collider physics programme to studies of high-energy $`\gamma \gamma `$ collision physics; see, for example, \[1–3\].

In such events we produce a colour-singlet state $`M`$ which is practically free from soft secondary particles. Moreover, if forward going protons are tagged we can reconstruct the ‘missing’ mass $`M`$ with good resolution, and so have an ideal means to search for new resonances and to study threshold behaviour phenomena. We have to pay a price for ensuring such a clean diffractive signal. In particular, the diffractive event rate is suppressed by the small probability, $`\widehat{S}^2`$, that the rapidity gaps survive soft rescattering effects between the interacting hadrons, which can generate secondary particles which populate the gaps \[4–12\].

In general, we may write the survival factor $`\widehat{S}^2`$ in a multi-channel eikonal framework in the form

$$\widehat{S}^2=\frac{_i\left|_i(s,b_t^2)\right|^2\mathrm{exp}\left(\mathrm{\Omega }_i(s,b_t^2)\right)d^2b_t}{_i\left|_i(s,b_t^2)\right|^2d^2b_t}$$
(2)
where the incoming proton is decomposed into diffractive eigenstates, each with its own opacity111Really we deal with a matrix $`\mathrm{\Omega }_{jj^{}}^{ii^{}}`$, where the indices refer to the eigenstates of the two incoming and two outgoing hadrons . $`\mathrm{\Omega }_i`$. The amplitudes $`_i(s,b_t^2)`$ of the process of interest may be different in the different diffractive eigenstates. They are expressed in impact parameter $`b_t`$ space at centre-of-mass energy $`\sqrt{s}`$. It is important to recall that the suppression factor $`\widehat{S}^2`$ is not universal, but depends on the particular hard subprocess, as well as on the kinematical configurations of the parent reaction, see, for example, .

Double-diffractive Higgs production,

$$ppp+H+p,$$
(3)
at the LHC, is a good example of illustrating the pros and cons of such exclusive processes. Let us assume a Higgs boson of mass $`M_H=120\mathrm{GeV}`$ and consider detection in the $`b\overline{b}`$ decay channel. The disadvantage is that to ensure the survival of the rapidity gaps in (3), the predicted cross section is low, $`\sigma 2\mathrm{fb}`$, corresponding to a survival factor $`\widehat{S}^2=0.02`$. The advantage is that, by tagging the outgoing protons, the signal-to-background ratio222The ratio quoted in (4) corresponds to a cut on the $`b`$ jet transverse momenta of $`p_t(b)>0.4M_H`$. is extremely favourable relative to other Higgs signals ,

$$\frac{\mathrm{signal}(Hb\overline{b})}{b\overline{b}\mathrm{QCD}\mathrm{background}}\stackrel{>}{}\mathrm{\hspace{0.33em}4},$$
(4)
if the missing mass resolution $`\mathrm{\Delta }M`$ obtained by the proton taggers is $`1\mathrm{GeV}`$. Indeed, with dedicated forward proton spectrometers at the LHC, the process $`ppp+H+p`$ may even be the light Higgs discovery channel.

For completeness, we note that for Higgs production via photon-photon fusion, the survival factor is much larger, $`\widehat{S}^20.86`$ , and the corresponding cross section $`\sigma _{\gamma \gamma }(Hb\overline{b})0.1\mathrm{fb}`$. In this case the signal-to-background ratio is

$$\frac{\mathrm{signal}(\gamma \gamma Hb\overline{b})}{b\overline{b}\mathrm{QED}\mathrm{background}}\frac{7\mathrm{GeV}}{\mathrm{\Delta }M_{b\overline{b}}}.$$
(5)
There are other physics reasons why it is desirable to tag the recoil protons in double diffractive processes. First, it offers a valuable experimental probe of the opacities $`\mathrm{\Omega }_i(s,b_t^2)`$ of the proton. The relevant Feynman diagrams for process (1) are shown in Fig. 1. There is appreciable interference between the amplitudes without and with the soft rescattering corrections, which are shown in Figs. 1(a) and 1(b) respectively. We will show that the interference depends on the transverse momenta $`\stackrel{}{p}_{1t},\stackrel{}{p}_{2t}`$ of the recoil protons and on the azimuthal angle $`\varphi `$ between these momenta. This dependence can be used to probe333Another possibility, to probe the opacity of the proton, is to study the process with rapidity gaps mediated by photon exchange . By varying the momentum transfer of the photon we sample different impact parameters $`b_t`$ and hence scan the opacity $`\mathrm{\Omega }(s,b_t^2)`$. the different soft rescattering models and the behaviour of the opacities $`\mathrm{\Omega }_i(s,b_t^2)`$.

Secondly, we need to understand, and if possible to predict, the $`\stackrel{}{p}_{it}`$ behaviour of the diffractive cross sections in order to plan experiments and to evaluate the acceptance and efficiency of the leading proton detectors at the LHC.

The content of the paper is as follows. In Section 2 we recall how the bare perturbative amplitude $`A^{(a)}`$ of Fig. 1(a) may be calculated, and then in Section 3 we illustrate the effect of rescattering corrections in terms of a simple model. Throughout the paper, it is safe to neglect the rescattering of the centrally produced system $`M`$ on either proton. This system (for example, a Higgs boson, high $`E_T`$ dijet, etc.) is massive and has small size, so, due to colour transparency , the cross section $`\sigma (Mp)`$ is very small. Rescattering is computed realistically in Section 4 and the predictions for a general double-diffractive process, $`ppp+M+p`$, are presented in Section 5. A detour is made in Section 6 to discuss double-diffractive light meson production, for which data already exist. Of course, this is beyond the region of validity of a perturbative QCD approach, but, surprisingly, the perturbative predictions agree qualitatively with interesting features of these data. Our conclusions are presented in Section 7.

## 2 The bare amplitude

The amplitude $`A^{(a)}`$ of Fig. 1(a), describing the high energy double-diffractive production of a heavy system $`M`$, can be expressed in terms of the generalised (skewed) unintegrated gluon densities $`f_g(x,x^{},t,Q_t,\mu )`$. Here $`\mu M/2`$ is the scale of the hard $`ggM`$ subprocess, and $`t`$ is the transverse momentum squared transferred through the ‘hard’ QCD Pomeron (that is the two-gluon system). Essentially the gluon distribution $`f_g`$ opens up and describes the internal structure of the ‘hard’ QCD Pomeron, whose exchange mediates the diffractive process (1).

For the exclusive reaction (1) the bare amplitude of Fig. 1(a) is, to single log accuracy, given by

$$A^{(a)}=\frac{1}{N_c^21}\frac{d^2Q_t}{Q_t^4}f_g(x_1,\mathrm{}Q_t,\mu )f_g(x_2,\mathrm{}Q_t,\mu )$$
(6)
where $``$ is the matrix element of the hard $`ggM`$ subprocess. For example, the cross section for the $`gggg`$ subprocess, relevant to high $`E_T`$ dijet production , is

$$\frac{d\widehat{\sigma }}{dt}=||^2=\frac{9}{4}\frac{\pi \alpha _s^2}{E_T^4}.$$
(7)
For small $`|x_ix_i^{}|`$, which is appropriate for high energy double diffraction, and $`t=0`$, the skewed unintegrated density $`f_g`$ can be calculated from knowledge of the conventional integrated gluon . The precise form of the $`t`$ dependence of $`f_g`$ is not well known. Recall, however, that $`f_g(\mathrm{}Q_t,\mu )`$ contains a Sudakov-like factor $`T(Q_t,\mu )`$ which reflects the chance that a gluon with transverse momentum $`Q_t`$ remains untouched in the evolution up to the hard scale $`\mu `$—a necessary condition for the survival of the rapidity gap, see, for example, . It is this $`T`$ factor which provides the infrared stability of the $`Q_t`$ integral of (6)444Moreover, the effective anomalous dimension of the gluon distribution additionally suppresses the contribution from the low $`Q_t`$ domain .. For example, for the production of a system of mass $`M\stackrel{>}{}\mathrm{\hspace{0.33em}100}\mathrm{GeV}`$ at the LHC, the saddle point of the $`Q_t^2`$ integration occurs at $`Q_t^2\stackrel{>}{}\mathrm{\hspace{0.33em}3}`$$`4\mathrm{GeV}^2`$. On the other hand, the transferred momenta satisfy $`|t_i|\stackrel{<}{}\mathrm{\hspace{0.33em}0.5}\mathrm{GeV}^2`$, which are small in comparison with $`Q_t^2`$. Therefore it is natural to assume the factorization

$$f_g(x,x^{},t,Q_t,\mu )=\beta (t)f(x,x^{},t=0,Q_t,\mu ),$$
(8)
where the whole $`t`$ dependence is given by the effective form-factor $`\beta (t)`$ of the QCD Pomeron-proton vertex. In other words, we separate the dependence of the $`ppp+M+p`$ cross section on the transverse variables ($`t_1,t_2`$ or $`\stackrel{}{p}_{1t},\stackrel{}{p}_{2t}`$) from the dependence on the ‘longitudinal’ variables (the initial $`pp`$ energy $`\sqrt{s}`$, the mass $`M`$ and rapidity $`y`$ of the system $`M`$). That is the bare amplitude is given by

$$A^{(a)}(\stackrel{}{p}_{1t},\stackrel{}{p}_{2t})=\beta (t_1)\beta (t_2)A_M$$
(9)
where $`t_ip_{it}^2`$, and where $`A_MA^{(a)}\left(p_{1t}=p_{2t}=0\right)`$ may be calculated from (6).

## 3 Absorption correction: a first look

To calculate the absorptive or soft rescattering amplitude $`A^{(b)}`$ of Fig. 1(b), it is convenient to use the momentum representation. We may neglect the spin-flip component in the proton-Pomeron vertex555This component is expected to be small and consistent with zero. If we note the similarity between the photon and Pomeron vertices then the magnitude of the isosinglet spin-flip amplitude is proportional to $`\left|\frac{1}{2}(\mu _p^a+\mu _n^a)\right|\stackrel{<}{}0.06`$, where the anomalous magnetic moments $`\mu ^a`$ of the neutron and proton cancel each other almost exactly.. We perform the detailed calculation in Section 4, but, first, we estimate the qualitative features of the rescattering effect by assuming, in this Section, that the amplitude for elastic proton-proton scattering, at energy $`\sqrt{s}`$ and momentum transfer $`k_t`$, has the simplified form

$$A_{pp}(s,k_t^2)=A_0(s)\mathrm{exp}(Bk_t^2/2).$$
(10)
From the optical theorem we have $`\mathrm{Im}A_0(s)=s\sigma _{pp}^{\mathrm{tot}}(s)`$, and for the small contribution of the real part it is sufficient to use $`\mathrm{Re}A_0/\mathrm{Im}A_00.12`$ in the energy regime of interest. $`B`$ is the slope of the elastic $`pp`$ differential cross section, $`d\sigma _{pp}/dt\mathrm{exp}(Bt)`$.

Using the above elastic $`pp`$ amplitude we may write the rescattering contribution, Fig. 1(b), to the $`ppp+M+p`$ amplitude as

$$A^{(b)}=i\frac{d^2k_t}{8\pi ^2}\beta (t_1)\beta (t_2)A_M\frac{A_0}{s}e^{Bk_t^2/2}$$
(11)
where now $`t_1(\stackrel{}{k}_t\stackrel{}{p}_{1t})^2`$ and $`t_2(\stackrel{}{k}_t+\stackrel{}{p}_{2t})^2`$. If we take an exponential form for the QCD Pomeron vertices,

$$\beta (t)=e^{bt/2},$$
(12)
then the integral in (11) can be evaluated, to give

$$A^{(b)}(\stackrel{}{p}_{1t},\stackrel{}{p}_{2t})=\frac{iA_0}{4\pi s(B+2b)}\mathrm{exp}\left(\frac{b^2|\stackrel{}{p}_{1t}\stackrel{}{p}_{2t}|^2}{2(B+2b)}\right)A^{(a)}(\stackrel{}{p}_{1t},\stackrel{}{p}_{2t}).$$
(13)
To gain insight it is useful to compute the numerical value of $`A^{(b)}`$ at the LHC energy using reasonable values of the parameters. We take $`b=4\mathrm{GeV}^2`$, $`B=20\mathrm{GeV}^2`$, $`\sigma _{pp}^{\mathrm{tot}}=100\mathrm{mb}`$ and, for the moment, neglect the real part, $`\mathrm{Re}A_0`$, of the $`pp`$ elastic amplitude. We obtain

$$A^{(b)}=0.73\mathrm{exp}\left(C|\stackrel{}{p}_{1t}\stackrel{}{p}_{2t}|^2\right)A^{(a)},$$
(14)
where $`C=0.29\mathrm{GeV}^2`$. Thus in the back-to back configuration with $`\stackrel{}{p}_{1t}\stackrel{}{p}_{2t}0.5\mathrm{GeV}`$, the absorptive correction $`A^{(b)}`$ completely cancels the bare amplitude $`A^{(a)}`$, and we predict a deep diffractive dip. Moreover we see that the position of the dip depends on the azimuthal angle $`\varphi `$ between the transverse momenta $`\stackrel{}{p}_{1t}`$ and $`\stackrel{}{p}_{2t}`$ of the tagged protons. For $`\varphi =180^{}`$ the momentum transfer occurs mainly through the elastic amplitude $`A_{pp}`$, with $`|t_1|`$ and $`|t_2|`$ minimized simultaneously, and hence the amplitude $`A^{(b)}`$ becomes larger.

## 4 Detailed treatment of rescattering corrections

To make a realistic calculation of the rescattering corrections we must improve the description of $`pp`$ soft interaction. We use the model of Ref. . It embodies (i) pion-loop insertions to the Pomeron trajectory, (ii) two-channel eikonal description of proton-proton rescattering and (iii) high mass diffractive dissociation. The parameters of the model were tuned to describe all the main features of the soft $`pp`$ data throughout the CERN-ISR to the Tevatron energy interval. In terms of the two channel eikonal the incoming proton is described by two diffractive eigenstates $`|\varphi _i`$, each with its own absorptive cross section.

The eigenstates were taken to have the same profile in impact parameter space, and absorptive cross sections

$$\sigma _i=a_i\sigma _0\mathrm{with}a_i=\mathrm{\hspace{0.33em}1}\pm \gamma ,$$
(15)
where $`\gamma =0.4`$ and $`\sigma _0`$ is defined in Ref. . That is the two channel opacity is

$$\mathrm{\Omega }_{jj^{}}^{ii^{}}=\delta _{ii^{}}\delta _{jj^{}}a_ia_j\mathrm{\Omega }.$$
(16)
The impact parameter representation of the elastic amplitude is thus

$$\mathrm{Im}\stackrel{~}{A}_{pp}(s,b_t)=s\left(1\frac{1}{4}\left[e^{(1+\gamma )^2\mathrm{\Omega }/2}+\mathrm{\hspace{0.25em}2}e^{(1\gamma ^2)\mathrm{\Omega }/2}+e^{(1\gamma )^2\mathrm{\Omega }/2}\right]\right).$$
(17)
When we allow for the extra $`(1\pm \gamma )^2`$ factors, which reflect the different Pomeron couplings to the two diffractive eigenstates in the $`ppp+M+p`$ production amplitude, that is in the right-hand part of Fig. 1(b), we obtain the effective amplitude of $`pp`$ rescattering,

$$\mathrm{Im}\stackrel{~}{A}_{pp}(s,b_t)=s\left(1\frac{1}{4}\left[(1+\gamma )^2e^{(1+\gamma )^2\mathrm{\Omega }/2}+\mathrm{\hspace{0.25em}2}(1\gamma ^2)e^{(1\gamma ^2)\mathrm{\Omega }/2}+(1\gamma )^2e^{(1\gamma )^2\mathrm{\Omega }/2}\right]\right).$$
(18)
The optical density $`\mathrm{\Omega }(s,b_t^2)`$ was given in Ref.  for Tevatron $`(\sqrt{s}=2\mathrm{TeV})`$ and LHC $`(\sqrt{s}=14\mathrm{TeV})`$ energies.

As before, we work in momentum space, and replace (10) by

$$A_{pp}(s,k_t^2)=\frac{1}{2\pi }d^2b_t\mathrm{\hspace{0.17em}4}\pi \stackrel{~}{A}_{pp}(s,b_t)e^{i\stackrel{}{k}_t\stackrel{}{b}_t}.$$
(19)
In this way we obtain a more realistic evaluation of the rescattering amplitude $`A^{(b)}`$ of Fig. 1(b). However, from the naïve evaluation of Section 3, we anticipate that there will still be a strong cancellation between the bare amplitude $`A^{(a)}`$ and the absorptive correction $`A^{(b)}`$, originating from the imaginary part of the elastic rescattering $`A_{pp}`$.

First, we must introduce the real part of $`A_{pp}(s,k_t^2)`$. However, it is not justified to use a constant ratio $`\mathrm{Re}A_{pp}/\mathrm{Im}A_{pp}\mathrm{\hspace{0.25em}0.12}`$ for $`k_t`$ away from zero. To account for the $`k_t`$ dependence we use the dispersion relation result, recalling that the total $`pp`$ cross section increases logarithmically with energy. Then the even-signature amplitude satisfies

$$\mathrm{Re}\stackrel{~}{A}_{pp}(s,b_t)=\frac{\pi }{2}s\frac{\left(\mathrm{Im}\stackrel{~}{A}_{pp}(s,b_t)/s\right)}{\mathrm{ln}s}.$$
(20)
To convert this relation from the impact parameter space amplitude $`\stackrel{~}{A}_{pp}(s,b_t)`$ to the momentum space amplitude $`A_{pp}(s,k_t^2)`$ we use (19).

Second, we should specify the form of the QCD Pomeron-proton vertex, $`\beta (t)`$. The most consistent choice is to take the dipole form used in Ref.

$$\beta (t)=\frac{1}{(1t/a_1)}\frac{1}{(1t/a_2)}.$$
(21)
For comparison we also evaluated the $`ppp+M+p`$ cross section using the alternative form

$$\beta (t)=e^{bt/2}$$
(22)
with $`b=4\mathrm{GeV}^2`$, which is consistent with the $`\gamma pJ/\psi +p`$ HERA data666The latest HERA data for $`J/\psi `$ elastic photoproduction prefer $`b=4.5\mathrm{GeV}^2`$ . However, for a heavier system, the smaller slope $`b=4\mathrm{GeV}^2`$ looks more reasonable. .

## 5 Predictions for tagged protons at the LHC

We are now in a position to predict the transverse momentum dependence of the outgoing protons in the double-diffractive production of a heavy system of mass $`M`$, that is in the process $`ppp+M+p`$. We show the results in the form

$$M^2\frac{^4}{yM^2^2p_{1t}^2p_{2t}}=M^2\frac{^2}{yM^2}F(\stackrel{}{p}_{1t},\stackrel{}{p}_{2t})$$
(23)
where $`M^2d^2/dydM^2`$ is the Pomeron-Pomeron luminosity given in Ref.  and where the factor $`F`$ contains the explicit $`\stackrel{}{p}_{1t}`$ and $`\stackrel{}{p}_{2t}`$ dependence. The luminosity, given in , was integrated over $`d^2p_{1t}d^2p_{2t}`$ with the assumption that the bare amplitude had an exponential $`t`$ behaviour,

$$A^{(a)}\mathrm{exp}\left(b_0(t_1+t_2)/2\right)$$
(24)
with $`b_0=4\mathrm{GeV}^2`$. In addition, the soft survival probability $`S^2`$ was averaged over the available $`t_1,t_2`$ domain. Here we unfold the luminosity to expose the $`\stackrel{}{p}_{1t}`$ and $`\stackrel{}{p}_{2t}`$ dependence. In order to be able to use the published $`M^2^2/yM^2`$ luminosity, we therefore compute

$$F(\stackrel{}{p}_{1t},\stackrel{}{p}_{2t})=\frac{\beta ^2(t_1)\beta ^2(t_2)}{S^2\pi ^2/b_0^2}\frac{^2S^2(\stackrel{}{p}_{1t},\stackrel{}{p}_{2t})}{^2p_{1t}^2p_{2t}}.$$
(25)
From the product of $`F`$, computed in this way, and the luminosity given in Ref. , we obtain the luminosity as a function of $`\stackrel{}{p}_{1t}`$ and $`\stackrel{}{p}_{2t}`$, as well as of $`y`$ and $`M`$. This resultant luminosity, (23), needs only be multiplied by the appropriate hard subprocess $`gg^{PP}M`$ cross section777The notation $`gg^{PP}`$ is to indicate that the hard gluons, which interact to form the system $`M`$, originate within overall (colourless) hard Pomeron exchanges. to obtain the differential cross section for any $`ppp+M+p`$ diffractive process. Various hard subprocess cross sections were listed, and discussed, in Ref. .

The factor $`F`$ is plotted in Figs. 2(a,b,c) as a function of $`p_{1t}`$ for three values of $`p_{2t}=0.2,0.4`$ and $`0.7\mathrm{GeV}`$ respectively. In each case the factor is shown for $`\varphi =0^{},90^{}`$ and $`180^{}`$, where $`\varphi `$ is the azimuthal angle between $`\stackrel{}{p}_{1t}`$ and $`\stackrel{}{p}_{2t}`$. The continuous and dashed curves are obtained using (21) and (22), respectively, for the QCD Pomeron-proton vertex. Recall that the model of Ref.  was fitted to ‘soft’ diffractive $`pp`$ data in the region $`|t|0.5\mathrm{GeV}^2`$, and so, strictly speaking, $`\beta (t)`$ of (21) should only be applied for $`p_t\stackrel{<}{}\mathrm{\hspace{0.33em}0.7}\mathrm{GeV}`$. However, we hope we can reliably evaluate rescattering corrections up to $`p_t1\mathrm{GeV}`$, as the typical values of $`k_t`$ (the momentum transferred through the elastic amplitude $`A_{pp}`$ in Fig. 1(b)), which are controlled mainly by the elastic slope $`B/2`$, are much less than $`1\mathrm{GeV}`$. As discussed before, the absorptive corrections are stronger in the back-to-back configuration; already for $`p_t<0.7\mathrm{GeV}`$ the $`\varphi =180^{}`$ curves reveal a rich dip structure.

In the final plot, Fig. 2(d), we compare the prediction for $`F`$ obtained using the elastic amplitude determined in the two-channel eikonal model of Ref. , with a naïve estimate based on a simple one-channel approach where the elastic $`pp`$ amplitude is given by the Gaussian formula of (10), that is the amplitude is described by single Pomeron exchange. However, we keep the parameters found in , that is $`\sigma _{pp}^{\mathrm{tot}}=102\mathrm{mb}`$, $`B=20.7\mathrm{GeV}^2`$ and $`\mathrm{Re}A_{pp}/\mathrm{Im}A_{pp}=0.12`$ at $`k_t=0`$. In both cases we use (22) for $`\beta (t)`$. The large difference between the realistic and naïve predictions demonstrates their sensitivity to the model used for soft rescattering.

In Fig. 3 we show the azimuthal dependence of the ‘soft’ survival factor

$$S^2(\stackrel{}{p}_{1t},\stackrel{}{p}_{2t})=\frac{|A^{(a)}+A^{(b)}|^2}{|A^{(a)}|^2},$$
(26)
as a function of the azimuthal angle $`\varphi `$, for different choices of $`p_{1t}`$ and $`p_{2t}`$. The rich structure of $`S^2`$ is apparent, which feeds through into the double-diffractive cross section. As anticipated, we observe a flatter behaviour in $`\varphi `$ for small $`p_{1t}`$ and $`p_{2t}`$, while for larger $`p_t0.7\mathrm{GeV}`$ a diffractive dip already occurs for $`\varphi 90^{}`$.

## 6 Application to double diffractive meson production

An interesting $`\varphi `$ behaviour has been observed by the WA102 collaboration at CERN  at lower energies ($`\sqrt{s}30\mathrm{GeV}`$) in fixed target central double-diffractive meson production,

$$ppp+X+p,$$
(27)
where a partial wave analysis of the $`X`$ channel allows the identification of a wide range of meson resonances.

It has been emphasized in Ref.  that within the Regge framework the Reggeon-Reggeon$``$meson vertex, $`V(RRX)`$, embodied in the amplitude $`A^{(a)}`$, may contain an azimuthal dependence like

$$|V|^2=1+a\mathrm{cos}\varphi +b\mathrm{cos}^2\varphi .$$
(28)
In fact, with an appropriate choice of parameters, a phenomenological description of the $`\varphi `$ dependences observed in the meson-production data can be achieved without any rescattering corrections .

For a hard subprocess at scale $`\mu =M/2`$ with $`Q_t^2t_i`$ we have no such dependence in the vertex $`V(IPIPM)`$; the coefficients in (28) satisfy $`a,b\stackrel{<}{}|t_i|/4Q_t^2`$. On the other hand, at least part of the azimuthal effects observed in the meson data  may originate in the rescattering corrections discussed in the present paper.

To study this further, we adopt the perturbative QCD viewpoint, which offers a dynamical basis in which to understand the structure of process (27). Of course, it is questionable to use perturbative QCD to describe the production of rather light mesons, via (27), where we have no hard scale. On the other hand it is natural to expect a smooth matching between the ‘soft’ and perturbative regimes. In this way we may obtain a qualitative interpretation of observed features the data. As we shall see below, this indeed turns out to be the case.

Recall that, in general, for forward going protons ($`p_{it}Q_t`$), two ‘hard’ QCD Pomerons can produce only a P-even state with the longitudinal projection of its spin $`J_z=0`$ . Also note that, as the QCD Pomerons are built from gluons, the underlying fusion subprocess $`gg^{PP}X`$ may provide a favourable environment for the production of exotic meson states containing gluons (such as glueballs, hybrids, etc.); the cross section of the $`gg^{PP}q\overline{q}`$ subprocess is much smaller than that of $`gg^{PP}gg`$, especially in the $`J_z=0`$ channel. Next for $`J_z=0`$ the vertex $`gg^{PP}X(2^{++})`$ is strongly suppressed if the $`2^{++}`$ state is a normal non-relativistic $`q\overline{q}`$ meson888The origin of this result can actually be traced to the absence of the $`\gamma \gamma `$ decay mode of $`2^{++}`$ positronium in the $`J_z=0`$ state , and then to the absence of the $`gg`$ decay of a $`J_z=0`$ $`q\overline{q}`$ system .. This is the result of gauge invariance. Indeed, for $`X`$ made from a non-relativistic $`q\overline{q}`$ pair, the fusion process $`gg^{PP}(q\overline{q})`$ looks like a single local vertex. The distance between the two gluon vertices is of the order of the inverse constituent quark mass ($`1/m_q`$), which is much smaller than the size of the ($`q\overline{q}`$) bound state. Now, the structure of the local $`gg^{PP}X(2^{++})`$ interaction is fixed by gluon gauge invariance. Then the requirement that the polarization tensor $`T_{\mu \nu }`$ of the $`2^{++}`$ meson satisfies $`T_{\mu \mu }=0`$ means that the vertex vanishes for $`J_z=0`$ . The consequence is that the forward double-diffractive production of normal quarkonium $`q\overline{q}(2^{++})`$ states should be suppressed , and so Pomeron-Pomeron fusion produces relatively more exotic (non-$`q\overline{q}`$) mesons, such as glueballs, $`(q\overline{q}g)`$ states, etc. In other words the process $`ppp+X(2^{++})+p`$ indeed acts as a filter for separating out exotic mesons from normal mesons .

The next observation is that mesons produced in the process $`ppp+X+p`$ by Pomeron-Pomeron fusion have larger transverse momenta, $`P_t`$, where $`P_t=|\stackrel{}{p}_{1t}+\stackrel{}{p}_{2t}|`$. The reasons are that (i) the $`J_z=0`$ selection rule is absent at larger $`p_{it}`$, (ii) for a $`1^{++}`$ meson, the vertex contains a factor $`ϵ^{\mu \nu \alpha \beta }p_{1t\mu }p_{2t\nu }`$ and so prefers larger $`p_{it}`$999By analogous arguments, the forward production of the unnatural spin-parity states, $`0^+`$ and $`2^+`$, should be strongly suppressed and also favours large $`p_{it}`$. This does not depend on whether or not the mesons are quarkonium states. Also, $`1^+`$ production would tend to occur at large $`p_{it}`$., and, finally, (iii) heavier (exotic) mesons101010This, of course, is also valid for $`0^{++}`$ mesons. tend to have larger $`P_t`$ due simply to kinematics.

Inspection of Fig. 3 shows that for small $`p_{it}`$ the cross section decreases with increasing $`\varphi `$, while for large $`p_{it}`$ it increases. In the latter domain the dominant contribution comes from the back-to-back configuration.

All the qualitative features described above are indeed observed in the available data  for process (27):

* Pomeron-Pomeron fusion in double-diffractive meson production may act as a glueball filter: the final state is enriched by non-$`q\overline{q}`$ mesons, which have smaller $`P_T`$ and are produced mainly with tagged protons in the $`\varphi =0`$ configuration;
* the normal $`q\overline{q}`$ light mesons have larger $`P_T`$ and their cross sections peak at $`\varphi =180^{}`$111111The preference for double-diffractive $`f_2`$(1270)-meson production in the $`\varphi >90^{}`$ domain has been observed at higher energies at the ISR .;
* the $`2^{++}`$ channel is produced mainly in the $`J_z=0`$ state .

These expectations can be confirmed by observing the double-diffractive production of heavier quarkonia, like $`\chi _c`$ and $`\chi _b`$, at the Tevatron and at RHIC. The heavier mesons sample smaller distances and so their production should be better described by perturbative QCD. Of course, $`\chi _c`$ is probably still not heavy enough, but nevertheless it would be interesting to compare $`\chi _c(2^{++})`$ production with the enhanced $`2^{++}`$ glueball production rate.

## 7 Conclusions

It is well known that, in general, absorptive effects in inelastic diffractive processes are much stronger than in the elastic amplitude (see, for example, ). Such rescattering clearly violates Regge factorization and leads to non-trivial correlations between the transverse momenta $`\stackrel{}{p}_{1t}`$ and $`\stackrel{}{p}_{2t}`$ of the forward going protons in processes of the type $`ppp+M+p`$. Measuring the $`p_{it}`$ and the azimuthal angle $`\varphi `$ distributions can provide an interesting possibility to probe the opacity $`\mathrm{\Omega }(s,b_t)`$ of the incoming proton and, moreover, to test the dynamics of soft rescattering.

One of the best examples to study these effects is exclusive high-$`E_T`$ dijet production, $`ppp+\mathrm{dijet}+p`$, where the cross section for the hard subprocess is large and well known.

Although questionable, the above perturbative formalism was applied to central double-diffractive meson production at lower energies, at which data exist. Surprisingly, the qualitative features of these data were reproduced.

## Acknowledgements

We thank Frank Close, J. Lamsa, Risto Orava and Albert de Roeck for interesting discussions. One of us (VAK) thanks the Leverhulme Trust for a Fellowship. This work was partially supported by the UK Particle Physics and Astronomy Research Council, by the Russian Fund for Fundamental Research (grants 01-02-17095 and 00-15-96610) and by the EU Framework TMR programme, contract FMRX-CT98-0194 (DG 12-MIHT).