text:
# Anisotropic relaxations introduced by Cd impurities in Rutile TiO2: first-principles calculations and experimental support

Also at \]LURE, Université Paris-Sud, Orsay, France.

## Abstract

We present an ab initio study of the relaxations introduced in TiO2 when a Cd impurity substitutes a Ti atom and the experimental test of this calculation by a perturbed-angular-correlation (PAC) measurement of the orientation of the electric-field gradient (EFG) tensor at the Cd site. The ab-initio calculation predicts strong anisotropic relaxations of the nearest oxygen neighbors of the impurity and a change of the orientation of the largest EFG tensor component, V33, from the $`[001]`$ to the $`[110]`$ direction upon substitution of a Ti atom by a Cd impurity. The last prediction is confirmed by the PAC experiment that shows that V33 at the Cd site is parallel to either the $`[110]`$ or the $`[1\overline{1}0]`$ crystal axis.

Metal impurities in oxides are a challenging problem not only from a fundamental point of view. A reliable calculation of impurity properties such as, e.g., their energy levels in semiconductor oxides is also of great technological importance. Impurities introduce structural atomic relaxations in the host and modify the electronic structure of the system and it is the interplay of these two effects which makes an adequate theoretical description very difficult. Recently, there have been several attempts to deal with this problem at the ab initio level impurities ; cdte\_blaha .

Experimental tests of the ab initio predictions are of fundamental importance for the evaluation of models of impurity-induced relaxations. The structural relaxations modify the local electron distribution at the impurity site. A quantity reflecting this distribution is the electric-field gradient (EFG). As the EFG decreases with 1/r3 from the producing charge density, it is very ”short-sighted” and therefore particularly sensitive to slight local changes. Consequently, the EFG at the impurity site carries information on the relaxations induced by the impurity and its experimental value may serve as test of the ab initio calculations. The EFG can be determined by a measurement of its electric hyperfine interaction with a quadrupole moment of the impurity nucleus.

In this paper we present an ab initio calculation of the relaxations at Cd impurities on substitutional Ti sites in TiO2 and the experimental test of this calculation by a measurement of the EFG tensor experienced by the Cd impurities. This particular host (TiO2)-impurity (Cd) combination was suggested mainly by two aspects. First, strength, symmetry, and orientation of the EFG tensor at Ti nuclei in the rutile structure of pure TiO2 are known from previous NMR measurements of the 49Ti quadrupole interaction in single crystals nmr\_tio2 . Ab initio calculations blaha\_tio2  are in agreement with this experimental result. Second, 111Cd is one of the most favorable probe nuclei for perturbed-angular-correlation (PAC) experiments and when introduced as an impurity into TiO2 is known to reside on substitutional Ti sites. PAC is a hyperfine spectroscopic technique that can be applied at very low probe concentrations and is therefore particularly well suited for the investigation of the EFG at isolated impurities. Strength (in terms of the absolute value of the largest principal-axis component, V33) and symmetry of the EFG tensor of 111Cd on Ti sites in TiO2 have been measured previously by 111Cd PAC in polycrystalline TiO2 lieb ; catchen .

Recently, the EFG of Cd in TiO2 has been calculated by Sato et al. akai  imposing the constraint of isotropic relaxations of the nearest oxygen neighbors and by some of us without this constraint our\_nqi . Although the absolute value of V33 provided by experiment is fairly well reproduced, the results of both calculations differ considerably: $`(i)`$ opposite signs are predicted and $`(ii)`$ according to our calculation the substitution of Ti by Cd should change the EFG orientation from the $`[001]`$ direction of pure TiO2 to the $`[110]`$ direction while according to Sato et al. the $`[001]`$ direction is not affected. These two approaches to the problem of a Cd-impurity in TiO2 have two shortcomings: $`(i)`$ the size of the supercell (SC) considered in both calculations (12 atoms) was very small for a dilute impurity calculation and $`(ii)`$ the charge state of the impurity was not taken into account in a self-consistent way. In the present work we have removed both shortcomings by using much bigger SCs and considering different impurity charge states in a fully self-consistent way.

The first-principles calculations were performed with the WIEN97 implementation of the Full-Potential Linearized-Augmented-Plane-Wave (FLAPW) method wien97 , which has been successfully used to determine the EFG in many compounds succes\_EFG ; cdte\_blaha . We worked in LDA approximation LDA . We have considered a $`2\times 2\times 3`$ supercell of 72 atoms consisting of twelve unit cells of TiO2 (see Fig. 1),

where one of the 24 titanium atoms is replaced by a Cd atom. Each Cd atom is therefore almost equidistant from its Cd images at a distance of around 9 Å. Substitution of a Ti by a Cd atom makes the resulting system metallic because of the lack of two electrons to fill up the $`p`$-O band. It is a question of central importance whether the real system we want to describe provides the lacking electrons (via an oxygen vacancy, for example) or not. In this letter we present the results of calculations assuming $`(a)`$ that extra electrons are not available (neutral impurity state: Cd0), and $`(b)`$ that the system somehow provides the lacking two electrons (charged impurity state: Cd-). In the situation $`(a)`$ we use the described SC performing calculations on a metallic system. To describe situation $`(b)`$ we add two electrons to the SC that we compensate with an homogeneous positive background in order to have a neutral cell to compute total energy and forces (Cd-(2e)). We have also simulated situation $`(b)`$ by replacing the two most distant oxygen atoms from Cd by two fluorine atoms in order to provide the lacking two electrons without introducing any artificial background (Cd-(fluorine)).

The value of the parameter RKMAX Soldner , which controls the size of the basis-set, was set to 6 and we introduced local orbitals (LO) to include Ti-3$`s`$ and 3$`p`$, O-2$`s`$ and Cd-4$`p`$ orbitals. Integrations in reciprocal space were performed using the tetrahedron method taking 8 k-points in the irreducible Brillouin zone for the metallic system $`(a)`$, and 2 k-points for the non-metallic ones $`(b)`$. Once self-consistency of the potential was achieved, the forces were computed forces , the ions were displaced newton , and new positions for Cd neighbors were obtained. The procedure was repeated until the forces on the ions were below 0.025 eV/Å. At the relaxed structure, the $`V_{ij}`$ elements of the EFG tensor were obtained from the $`V_{2M}`$ components of the lattice harmonic expansion of the self-consistent potential efg . We have assumed that relaxations preserve the point-group symmetry of the SC, which restricts O1 and O2 displacements to YZ plane and X axis direction, respectively (see Fig. 1). In order to check the stability of the solution obtained, at the end of the relaxation process we performed new calculations with O1 and O2 atoms displaced from their symmetry positions and verified that this solution is a minimum.

In Table 1 we show the results for the relaxation of the 6 nearest oxygen neighbors of the Cd impurity. We see that for both impurity charge-states the relaxations are quite anisotropic, with the Cd-O1 distance bigger than the Cd-O2 distance, opposite to the initial unrelaxed structure. This conclusion contradicts the assumptions of previous studies of this system akai ; lieb  and confirms the tendency predicted in the first calculation performed by some of us with a much smaller SC our\_nqi . Anisotropy in the relaxations of the nearest oxygen neighbors of the Cd-impurity can be understood by inspection of Fig. 2.

Stretching of Cd-O2 bond implies a considerable shortening in Ti-O2 bonds. However, stretching of Cd-O1 bond affects the Ti-O1 bonds not so much, since the structure is more open in this direction. So, at the end of the relaxation, the Cd-O1 bond stretches almost twice as much as the Cd-O2 bond. Relaxations and EFGs are very similar for the two ways of simulating the charged state of the impurity indicating that both approaches are well suited to deal with this problem. These results for EFG agree reasonably well with available experimental information lieb ; catchen  and are very different from the ones obtained for the neutral state of the impurity, Cd0. The huge difference obtained in the $`V_{ii}`$ components for the two impurity charge-states is mainly due to the filling of the impurity state at the Fermi level which has an important component of Cd$`d_{yz}`$ character PRB\_nos . The high $`\eta `$ value obtained for Cd0 is too far from experimental findings suggesting that Cd as an impurity in TiO2 is in a charged state. The results obtained for Cd-(2e) changes very slightly when the coordinates of atoms beyond the Cd-nearest neighbors are relaxed (see Table 2).

In addition we also performed FLAPW calculations for the system Cd-(2e) imposing the constraint of isotropic relaxations of O1 and O2 atoms and obtained V33 pointing in -direction and a high $`\eta `$ value of 0.91 confirming for a converged-size SC that an isotropic relaxation is not consistent with the available experimental data.

To check the accuracy and convergence of our results we have performed additional calculations at the relaxed structure of case Cd-(2e) with $`R_C`$=2.5 Å considering $`RK_{MAX}`$=7, including LO for Cd-$`4d`$ orbitals, increasing the mesh in $`\stackrel{}{k}`$-space, using GGA approximation gga  for exchange-correlation potential and also considering $`2\times 2\times 2`$ and $`2\times 2\times 4`$ SCs. These checks show that d(Cd-O1) and d(Cd-O2) are converged within 0.01 Åand V33 and $`\eta `$ within 0.5 1021V/m2 and 0.1, respectively. Details will be presented in a forthcoming paper PRB\_nos .

In order to check the change in orientation of V33 from the Z to the Y direction upon substitution of a Ti atom by a Cd impurity, we performed PAC experiments in TiO2 single crystals with the 111Cd probe. The PAC technique measures the probability W($`𝐤_1`$,$`𝐤_2`$,t) to detect the second $`\gamma `$-ray of a nuclear-decay $`\gamma `$-$`\gamma `$ cascade at the direction $`𝐤_2`$($`\theta _2`$,$`\varphi _2`$) after a time $`t`$ the first $`\gamma `$-ray was detected at the direction $`𝐤_1`$($`\theta _1`$,$`\varphi _1`$), where $`\theta _i`$ and $`\varphi _i`$ describe the propagation direction of $`\gamma _1`$ and $`\gamma _2`$ with respect to a chosen quantization axis. The angular correlation may be modulated in time by hyperfine interactions which induce transitions between the m-states of the intermediate level of the cascade, changing their population distribution (relative to $`\gamma _1`$) and thus the angular distribution of $`\gamma _2`$. The correlation function W($`𝐤_1`$,$`𝐤_2`$,t) therefore carries information on the strength, symmetry and time dependence of the hyperfine field and is sensitive to the orientation of the interaction with respect to $`𝐤_1`$ and $`𝐤_2`$.

The orientation dependence is easily understood by considering two simple cases. Let us assume an axially symmetric EFG and choose the direction of $`𝐤_1`$ as quantization axis. If we orient the EFG parallel to $`𝐤_1`$, the spin precession does not change the projection of I on $`𝐤_1`$. As a consequence the population of the m-states with respect to $`𝐤_1`$ and thus the correlation function remain constant in time. In case the symmetry axis of the EFG is perpendicular to the direction of the first ($`𝐤_1`$) detector, the spin precession causes a periodic oscillation of the projections of I on $`𝐤_1`$, corresponding to transitions between the different $`m`$-states relative to $`𝐤_1`$, and a periodic time-modulation of the intensity seen by the second detector results, which is invariant under rotation around the symmetry axis of the interaction. Details of the PAC theory can be found in Ref. TDPAC .

For the PAC experiment, we used a commercially available TiO2 single crystal with the dimensions 10x10x0.5 mm. The crystal was cut with the $`[001]`$ axis perpendicular to the 10x10 mm faces and the $`[110]`$ axis perpendicular to the 10x0.5 mm-sides. The mother isotope 111In of the 111Cd PAC probe was implanted at room temperature into the 10x10 mm faces with an energy of 160 keV to a dose of $``$1013 ions/cm2. Following the implantation the single crystal was annealed at 800C for 4 h in order to remove any radiation damage. After this treatment at least 90$`\%`$ of the probes were found to reside on the regular Ti site.

The time-differential perturbed-angular correlation of the 175 and 247 keV $`\gamma `$-$`\gamma `$ transition of the cascade 111Cd was measured with a standard four BaF2 detectors PAC spectrometer for several orientations of the crystal relative to the detectors. The spectra shown in Figs. 3(a) and 3(b) were obtained with the $`[001]`$-axis normal to the detector plane. In the case of Fig. 3(a) the $`[110]`$-axis was parallel to the bisector of two detectors, in fig. 3b it was parallel to one detector. The visual inspection of these spectra already reveals the EFG orientation:

Clearly, the rotation around the $`[001]`$ direction causes substantial changes of the PAC pattern. This observation immediately discards the $`[001]`$ direction as the direction of V33 because, as pointed out above, in the case of an axially symmetric EFG the PAC pattern should be invariant under rotation around the symmetry axis of the interaction. In the present case the EFG presents a slight axial asymmetry ($`\eta =0.18`$). But even then, the changes upon rotation around the V33 axis normal to the detector plane should be imperceptibly small. This is easily checked by calculating the perturbation factor for this geometry.

If in a crystal of tetragonal symmetry such as TiO2 the EFG is not parallel to the $`[001]`$ axis, one must expect at least two identical EFG tensors with equal weight but different orientations of the principal axis system. The PAC pattern obtained with the $`[110]`$-axis at 45 to the detectors (Fig. 3(a)) can be described with a single site, i.e., all 111Cd probes are subject to the same EFG implying that the two expected EFG orientations must be equivalent. With the $`[110]`$-axis along the bisector of the detectors, two equivalent orientations can be obtained only with V33 in the basal plane parallel to either the $`[110]`$\- or $`[1\overline{1}0]`$-axis.

A fit of the exact theory barradas  to the spectrum of Fig. 3(a) gave the values $`\theta 8^{}`$ and $`\varphi `$ = 43(2) for the angle between V33 and the detector plane and the angle between $`𝐤_1`$ and the projection of V33 onto the detector plane, respectively. We estimate that we adjusted the 10x10mm face of the TiO2 single crystal parallel to the detector plane within $`\pm 5^{}`$. So, the value of $`\theta =0^{}`$ is compatible with the experimental accuracy.

The conclusion of V33 lying in the basal plane parallel to either the $`[110]`$\- or $`[1\overline{1}0]`$-axis is confirmed by the spectrum in Fig. 3(b). In this case the two EFG orientations are not equivalent. As discussed above, with the EFG parallel to the $`𝐤_1`$ detector, the angular correlation remains unperturbed while the EFG normal to $`𝐤_1`$ produces a periodic modulation. The superposition of these two configurations leads to a PAC pattern consisting of a small-amplitude oscillation on a constant background. The oscillation is mainly caused by the second configuration, but the exact theory shows that a small oscillatory contribution also comes from the slight asymmetry of the EFG. The least-squares fit to the spectrum in Fig. 3(b) gave the angles $`\varphi _1=3(1)^{}`$ and $`\varphi _2=90(1)^{}`$. The value of the angle $`\theta `$ is consistent with the one obtained from the spectrum in Fig. 3(a). For both spectra the fitted values $`\nu _Q=107.1(3)`$ MHz and $`\eta =0.18(1)`$ are in good agreement with those previously reported lieb ; catchen .

In summary, we have shown that the inclusion of Cd-impurities in TiO2 distorts the oxygen octahedron surrounding the impurity, producing a change in the EFG tensor orientation. Such a change does not result if an isotropic relaxation of the oxygen octahedron is assumed. According to our calculation, the asymmetry of the EFG tensor suggests that Cd in TiO2 is not neutral, but in a charged state. We therefore conclude that a proper theoretical approach of electronic properties of metal impurities in oxide semiconductors should consider the charge state of the impurity and the impurity-induced relaxations without constraints in a fully self-consistent way.

This work was supported by CONICET, ANPCyT, F. Antorchas (Argentina), and DFG (Germany). MR is indebted to J-P. Dallas (CECM, Vitry) for the single crystals orientation and to A. Traverse (LURE) for valuable support on this project. We thank R. Vianden (ISKP) and A.F. Pasquevich (UNLP) for critical comments.


format: text

text:
# Untitled Document

Radiophysics and Quantum Electronics, Vol. 41, No. 8, 1998, 669–673

A DIFFERENTIAL METHOD OF MAXIMUM ENTROPY

> A.T.Bajkova
>
> We consider a differential method of maximum entropy that is based on the linearity of Fourier transform and involves reconstruction of images from the differences of the visibility function. The efficiency of the method is demonstrated with respect to the recovery of source images with bright components against the background of a sufficiently weak extended base. The simulation results are given along with the maps of an extragalactic radio source 0059+581, which were obtained using the standard and differential methods of maximum entropy for three observation dates and show that the principle of differential mapping allows us to increase considerably the dynamic interval of images.

1. THE PRINCIPLE OF DIFFERENTIAL MAPPING

The idea of differential mapping is not new. At present, it is realized most comprehensively in the program package ”DifMap”  at the California Institute of Technology, where the ”cleaning” procedure (CLEAN) is used as deconvolution operation.

The method of differential mapping is based on the fundamental property of the Fourier transformation. According to this method, the bright components of the source, which were reconstructed at the first stage, are subtracted from the initial visibility function, the subsequent reconstruction is performed using the residual visibility function, and the reconstruction results are added at the final stage.

In the case of CLEAN, such a mapping method mainly influences the convergence rate. However, if the maximum entropy method (MEM), which has pronounced nonlinear properties, is used as the deconvolution operation, the principle of differential mapping can improve the reconstruction quality, particularly when the source has bright compact components against the background of a sufficiently weak extended base. An example of such improvement due to subtraction of the bright component from the visibility function is shown in .

What is the cause of the reconstruction quality improvement? The fact is that after the subtraction of bright components, which were reconstructed at the first stage using the MEM, from the initial visibility function we obtain the residual visibility function, in which the proportion of the weak extended component becomes larger. (The subtraction result is called the first-order residual visibility function.) Therefore, we artificially decrease the dynamic interval of the map, which corresponds to the residual visibility function, and, thus simplify the image reconstruction at the second stage.

If the bright components of the image reconstructed at the second stage are subtracted from the first-order residual visibility function, we obtain a second-order residual visibility function, to which a map with a still lower dynamic interval corresponds. Using this procedure, we obtain residual visibility functions of an increasingly high order.

To obtain the desired source map, the map reconstructed at the last stage must be summed up with all the components subtracted from the visibility function at the previous stages.

Formally, the algorithm of two-stage differential mapping is formulated as

$$F_{vis}^{\left(2\right)}=F_{vis}^{\left(0\right)}F_{vis}^{\left(1\right)},$$

where $`F_{vis}^{(0)}`$ is the initial visibility function of the desired brightness distribution over the source $`I(x)`$; $`F_{vis}^{(1)}`$ is the visibility function corresponding to the bright component $`I^{(1)}(x)`$ reconstructed at the first stage; $`F_{vis}^{(2)}`$ is the residual visibility function corresponding to the source map $`I^{(2)}(x)`$ reconstructed at the second stage.

The resulting map has the form

$$I\left(x\right)=I^{\left(1\right)}\left(x\right)+I^{\left(2\right)}\left(x\right).$$

2. SIMULATION RESULTS

In Figs. 1 and 2 we show the results of simulation of the principle of differential mapping using the generalized maximum entropy method (GMEM)  for real images with both positive and negative values (the advantages of the GMEM compared with the MEM are discussed below in Section 3).

The experiment shown in Fig. 1 was carried out for three model sources in the form of one or several bright point components against the background of a weaker extended base that has a Gaussian form. In the figures we give the following source names: ”Model-1”, ”Model-2”, and ”Model-3”. For ”Model-1” the base amplitude is 20% of the amplitude of the bright point-like component, while for ”Model-2” it is only 2%. ”Model-3”, which has four bright components against the background of a 20% Gaussian, is more complicated. To form the visibility function, we used the symmetric coating of the $`UV`$ surface with almost uniform filling over the entire aperture, which included 98 points. To study the reconstruction improvement due to use of the differential mapping principle in its pure form, no noise was added to the visibility function.

In Figs. 1 and 2, for each vertically located source, we show the maps that appear horizontally in the following order (from left to right): the map of the model source, a ”noisy” image resulting from the inverse Fourier transform of the initial visibility function, an image reconstructed at the first stage of differential mapping using the GMEM, and the resulting map of differential mapping. (It must be noted that the bottom level of the contour line corresponds to 1% of the peak value of the flux in all the above images.)

Obviously, the maps obtained at the first stage can be treated as the limiting performance of the GMEM applied to the initial visibility function. The subtraction of the bright components reconstructed at the first stage from the visibility function allows us to improve significantly the reconstruction quality at the second stage. In Table 1 we show the quantitative properties of the image reconstruction quality from Fig. 1 after the first and second stages of the differential mapping.

As is obvious from the table, compared with the standard GMEM, the use of the differential GMEM allows us to increase the signal-to-noise ratio (SNR) of the output image by a factor of 3 to 5, such that the maximum effect is achieved for the source ”Model-2”, which has maximum dynamic interval.

In Fig. 2 we show the results of reconstruction of images of the sources ”Model-1” and ”Model-2” for the case of another filling of the $`UV`$ plane, which includes 98 points distributed uniformly inside the central region. The above region is approximately a quarter of the aperture. This example of mapping is more complicated compared with the previous case because of the total absence of the high spatial harmonics of the source in the visibility function. As is obvious from the figure, in this case the differential GMEM also gives higher-quality results than the standard GMEM.

Table 1

Characteristics of reconstructed images

|  | Maximum | SNR | Maximum | SNR |
| --- | --- | --- | --- | --- |
| Source | error |  | error |  |
|  | at the 1st stage | at the 1st stage | at the 2nd stage | at the 2nd stage |
| Model-1 | 0.4701 | 3.84 | 0.0572 | 15.86 |
| Model-2 | 0.0467 | 11.43 | 0.0088 | 52.00 |
| Model-3 | 0.1890 | 4.10 | 0.0594 | 13.20 |

3. MAPPING FOR THE RADIO SOURCE 0059+581

To demonstrate the potentialities of the differential generalized method of maximum entropy used for the real data processing, we show the results of mapping the extragalactic radio source 0059+581, which is often used as a reference source in the VLBI geodetic programs . This source is of astrophysical interest because it shows fast variation of both the total flux and the structure. The maps that are constructed using the GMEM for a number of dates in the interval from June 1994 to December 1995 are shown in Fig. 3. Variation of the total flux in the same interval is presented in Fig. 4 in .

In the time interval from the middle to the end of 1994 we observed an almost linear decrease in the total flux from $``$ 4 Ja to 1.5 to 2 Ja. As is obvious from Fig.3, we failed to obtain a map with dynamic range sufficient for detecting an extended element from the data of June 26, 1994 when the flux was maximal over the interval in question. Using the standard GMEM, we detected extended elements only starting from October 4, 1994 when their share with respect to the total flux became sufficient to be detected by the selected reconstruction method. In Fig. 4, for comparison, we show maps obtained using both the standard and differential GMEMs for the three dates on which the total flux was sufficiently large and the flux share corresponding to extended elements was insignificant. As is obvious from the figure, the use of differential mapping allowed us to increase the dynamic range of the maps such that extended structures became very distinguished. It must be noted that the minimal level of the contour line corresponds to 1% of the peak value of the flux everywhere in the images. The parameters of the maps obtained using the differential GMEM are given below in Table 2.

Table 2

Map parameters for the source 0059+581

| Date | Total | Peak |
| --- | --- | --- |
|  | Flux(Ja) | flux(Ja) |
| 25.06.94 | 4.23 | 2.40 |
| 25.08.94 | 3.40 | 1.76 |
| 04.10.94 | 3.28 | 1.86 |

The main criterion for choosing the method of maximum entropy (either classical or generalized, the MEM or the GMEM) is the type of the desired distribution. If the desired image is real bipolar (with positive and negative values) or complex with both real and imaginary parts, only the GMEM  can be used.

If it is known a priori that the image is real and nonnegative, the classical MEM can be used. However, one should keep in mind that in case of large errors in the data, the GMEM allows us to obtain  higher-quality maps with a much lower level of nonlinear distortions, as compared with the MEM. At present, the researchers of the IAA RAS carry out research in making maps of the extragalactic sources using the observation data from the astrometric and geodetic VLBI programs. These data were not obtained directly for astrophysical mapping, and, as a rule, have poor calibration. Therefore, to obtain satisfactory maps, the GMEM is recommended instead of the MEM.

Another advantage of the GMEM originates from the fact that an image with negative values can correspond to the residual visibility function obtained after the subtraction of the bright component that is reconstructed at the first stage of the differential algorithm from the initial visibility function. This is possible if the bright component was reconstructed at the first stage with an overestimated amplitude, which is quite typical of any nonlinear method. To rule out the unwanted image distortions when a nonnegative solution is sought using the data assuming the presence of negative components, we should use the generalized method of maximum entropy for real bipolar images rather than the classical method.

4. CONCLUSIONS

Compared with the traditional methods, the nonlinear methods of differential mapping using the maximum entropy method as the deconvolution operation allow us to improve substantially the reconstruction quality of source images containing bright point-like components against the background of weak base. The maps obtained using the differential method of maximum entropy are characterized by a higher dynamic range. To eliminate possible nonlinear distortions in the case of differential mapping, the generalized method of maximum entropy is preferred to the classical method.

This work was supported by the Russian Foundation for Basic Research, grant No 96-02-19177 and the Ministry of Science’s program ”Astronomy. Fundamental Space Research”, grant No 2.1.1.3.

REFERENCES

1. G.Taylor, The Difmap Cookbook, California Institute of Technology, Pasadena (1994).

2. T.Cornwell, ”Very Long Baseline Interferometry and the VLBA”, in: ASP Conference Series (J.A. Zensus, P.J. Diamond, and P.J. Napier, eds) vol.82, 227 (1995).

3. A.T.Bajkova, Commun. IAA RAS, No. 58 (1993).

4. A.T.Bajkova et al, Trans. IAA RAS, No. 1, 22 (1997).

5. T.B.Pyatunina, A.T.Bajkova et al, Trans. IAA RAS, No. 1, 64 (1997).

6. A.T.Bajkova, Izv. Vyssh. Uchebn. Zaved., Radiofiz., 34, No. 8, 919 (1991).

7. A.T.Bajkova, Izv. Vyssh. Uchebn. Zaved., Radiofiz., 38, No. 12, 1267 (1995).


text:
# Ultracold neutrons, quantum effects of gravity and the Weak Equivalence Principle

## Abstract

We consider an extension of the recent experiment with ultracold neutrons and the quantization of its vertical motion in order to test the Weak Equivalence Principle. We show that an improvement on the energy resolution of the experiment may allow to establish a modest limit to the Weak Equivalence Principle and on the gravitational screening constant. We also discuss the influence of a possible new interaction of Nature.

preprint: DF/IST-5.2002

Quite recently, Nesvizhevsky and collaborators have reported on the evidence of quantization of the vertical motion of neutrons bound by the gravitational field of the Earth Nesvizhevsky . Even though the conceptual aspects of the experiment were discussed back in the seventies Luschikov , concrete steps towards the goal of the final experiment were realized more recently Nesvizhevsky1 . The experiment consists in allowing ultracold neutrons generated by a source at the Institute Laue-Langevin reactor (Grenoble) to fall towards a horizontal mirror due to the influence of the Earth’s gravitational field. This potential confines the motion of the neutrons which no longer move continuously in the vertical direction, but rather jump from one height to another, as predicted by quantum mechanics. This impressive experiment complements the Collela, Overhauser and Werner experiment Collela  where neutrons were split and let interfere with the gravitational field, even though in that situation the neutrons were not quantum states. In the Grenoble experiment Nesvizhevsky , the minimum energy of the detected neutrons is $`1.4\times 10^{12}eV`$, corresponding to a vertical velocity of $`1.7cm/s`$. In that work it is stated that a more intense beam and an enclosure mirrored on all sides could lead to an improved energy resolution, down to $`10^{18}eV`$, if the neutron can be kept confined all its lifetime. In Ref. Nesvizhevsky  besides the ground state, three excited states were determined, although with reduced accuracy.

An earlier theoretical work Bertolami1  suggested the use of a version of the Collela, Overhauser and Werner experiment for establishing bounds on the parameters of a possible new interaction of Nature. In this work however, we shall show that an upgraded version of the experiment of Nesvizhevsky and collaborators does not yield any bound on the strength of new interaction of Nature and explore various implications of a highly resolved neutron spectrum.

The Schrödinger problem of a neutron submitted to the local gravitational potential, $`V=mgx`$ has well known solutions (here $`m`$ is the neutron’s gravitational mass and $`x`$ the vertical height). The eigen-solutions of $`\widehat{H}\mathrm{\Psi }_n=E_n\mathrm{\Psi }_n`$, can be expressed in terms of Airy functions, $`\mathrm{\Phi }`$, Landau :

$$\mathrm{\Psi }_n(z)=A_n\mathrm{\Phi }(z),$$
(1)
with eigenvalues determined by the roots of the Airy functions:

$$\alpha _n=\left(\frac{2}{mg^2\mathrm{}^2}\right)^{1/3}E_n.$$
(2)
The variable $`z`$ has a simple relation to the physical height $`x`$ which for the Nesvizhevsky et al. experiment is of order $`\mu m`$:

$$z=\left(\frac{2m^2g}{\mathrm{}^2}\right)^{1/3}(x\frac{E_n}{mg}).$$
(3)
The normalization of the wavefunction $`A_n`$ can be explicitly determined in terms of the integral

$$A_n^2=\left(\frac{\mathrm{}^2}{2m^2g}\right)^{1/3}_{\alpha _n}^{\mathrm{}}\mathrm{\Phi }(z)^2𝑑z.$$
(4)
First of all, directly from Eq. (2), one can extract a precision for the measurement of the local gravity of $`\mathrm{\Delta }g=6\times 10^9`$ ms-1, under the assumption that $`\mathrm{\Delta }E=10^{18}`$ eV. We will show that other effects superpose this uncertainty.

Let us then establish to which extent the Weak Equivalence Principle can be tested once the energy resolution of the Nesvizhevsky et al. experiment is improved. For that one has to realize that in the neutron’s Hamiltonian two masses are present. The neutron’s inertial mass, $`m_i`$, in the kinetic term and the gravitational mass in the gravitational potential. Distinguishing these two masses implies the relationship (2) is slightly altered as $`m`$ must be replaced by $`m^2/m_i`$. It then follows that the corresponding uncertainty in the energy due to the difference between inertial and gravitational masses is given by

$$\frac{m_im}{m}=\frac{3}{2}mg^2\mathrm{}^2\alpha _n^3\frac{\mathrm{\Delta }E}{E_n^4}.$$
(5)
Hence, for the ground state, $`E_0=1.44peV`$ Nesvizhevsky1  and $`\alpha _0=2.338`$, and if $`\mathrm{\Delta }E=10^{18}eV`$ is attained, one can establish a somewhat modest bound:

$$\frac{m_im}{m}<1.9\times 10^6.$$
(6)
The situation is somewhat different in what concerns a new force of Nature. The exciting prospect of a new fundamental interaction, beyond the already known four interactions, sparked in 1986 from the claim that the original Ëotvös experiment, designed to verify the equality of gravitational and inertial masses, revealed evidence of a new force with sub-gravitational strength (see Ref. Fishbach1  and references therein for an exhaustive discussion). Despite the fact that Ëotvös data turned out not be conclusive, the claim has stimulated a great deal of theoretical discussion (see e.g. Nieto  for a complete set of references) as well as the repetition of old experiments using new technology. Actually, the most stringent limit on the equality of inertial and gravitational masses (of $`Cu`$ and $`Pb`$) has arisen in this context and with an accuracy of $`5\times 10^{13}`$ Adelberger .

The simplest way a new interaction or a fifth force could arise would be through the exchange of a light boson coupled to matter with sub-gravitational strength. This could originate from various physical models at the Planck scale such as the extended supergravity theories after dimensional reduction Nieto ; Scherk , the compactification of $`5`$-dimensional generalized Kaluza-Klein theories that include gauge interactions at higher dimensions Bars  and also from string/M-theory. A common feature of these schemes is the appearance of a new Yukawa type modification in the interaction energy, $`V(r)`$, between two point masses $`m_1`$ and $`m_2`$:

$$V(r)=\frac{G_{\mathrm{}}m_1m_2}{r}[1+\alpha _5\mathrm{exp}(r/\lambda _5)],$$
(7)
where $`r=|\stackrel{}{r}_2\stackrel{}{r}_1|`$ is the distance between the masses, $`G_{\mathrm{}}`$ is the gravitational coupling for $`r\mathrm{}`$, $`\alpha _5`$ and $`\lambda _5`$ are the strength and the range of the new interaction. Of course, $`G_{\mathrm{}}`$ is related with the Newtonian gravitational constant. Indeed, the force associated with Eq. (7) is given by:

$$\stackrel{}{F}(r)=V(r)=\frac{G(r)m_1m_2}{r^2}\widehat{𝐫},$$
(8)
where

$$G(r)=G_{\mathrm{}}[1+\alpha _5(1+r/\lambda _5)\mathrm{exp}(r/\lambda _5)].$$
(9)
A particularly interesting implication of the mentioned approaches to the Planck scale physics is that the coupling $`\alpha _5`$ is not an universal constant, but instead a parameter dependent on the chemical composition of the test masses as first pointed out in Ref. Lee . This dependence comes about if one assumes that the new bosonic field couples to the baryon number $`B=Z+N`$ (the sum of protons and neutrons) and would imply a clear violation of the Weak Equivalence Principle.

Several experiments have been performed in order to establish the parameters of a new interaction based on the idea of a composition-dependence differential free fall of bodies (see Refs. Fishbach1 ; Nieto  for discussions on those leading to the most stringent limits). The current data is entirely compatible with predictions of Newtonian gravity from both, composition-independent or composition-dependent, experiments. The bounds on parameters $`\alpha _5`$ and $`\lambda _5`$ can be summarized as follows:

* Laboratory experiments devised to measure deviations from the inverse-square law are sensitive to the range $`10^2m<\lambda _5\text{ }<1m`$ and constrain $`\alpha _5`$ to be smaller than $`10^4`$;
* Nucleosynthesis bounds imply that $`\alpha _5\text{ }<4\times 10^1`$ for $`\lambda _5\text{ }<1m`$ Bertolami2 ;
* Gravimetric experiments sensitive in the range of $`10m\text{ }<\lambda _5\text{ }<10^3m`$ suggest $`\alpha _5<10^3`$;
* Satellite tests probing ranges about $`10^5m\text{ }<\lambda _5\text{ }<10^7m`$ indicate that $`\alpha _5<10^5`$;
* Radiometric data of the Pioneer 10/11, Galileo and Ulysses spacecrafts suggest the presence of a new force with parameters $`\alpha _5=10^3`$ and $`\lambda _54\times 10^{13}m`$ Anderson .

It is striking that, for $`\lambda _5<10^3m`$ and $`\lambda _5>10^{13}m`$, $`\alpha _5`$ is essentially unconstrained. The former range arises in higher dimensional superstring motivated cosmological solutions, where matter fields that are related to open string modes, lie on a lower dimensional brane, while gravity propagates in the bulk Polchinski . In these scenarios, the $`d`$ extra dimensions are not restricted to be small Antoniadis  and the fundamental $`D`$-dimensional scale, $`M_D`$, with $`D=4+d`$, can be considerably smaller than the 4-dimensional Planck scale. Assuming that these quantum gravity effects are just beyond current detection capability in accelerators, then $`M_D\text{ }>few\text{TeV}`$ and hence modifications to Newtonian gravity will occur in the short range region, $`\lambda _5<10^3m`$. This range also emerges if one assumes that scalar Beane  or vector/tensor Bertolami3  excitations are associated with the observed vacuum energy density.

As a matter of fact, the millimeter range has been recently available for experimental verification and as a consequence the effects associated to $`d`$ extra dimensions, that imply that $`\alpha _5=d+1(2d)`$ for a circle (torus) topology Floratos , have been ruled out down to $`\lambda _5\text{ }<0.2mm`$ yielding, for two extra dimensions, the bound on the typical energy scale, $`M_63.5\text{TeV}`$ Hoyle .

It should be mentioned that cosmological considerations tend to yield higher bounds for $`M_D`$. For instance, requiring that nucleosynthesis yields are not affected by brane related effects, imply that $`M_510\text{TeV}`$ Cline  and hence $`\lambda _5<10^4m`$. On the other hand, matching inflationary observable quantities in the context of supergravity models on the brane, yield much higher bounds for $`M_5`$, typically, $`M_5\text{ }>10^{13}\text{GeV}`$ for supergravity mass term chaotic inflation Bento1  and $`M_5\text{ }<10^{16}\text{GeV}`$ for supergravity inflationary models where the potential has the form $`V=V_0[1+c_n(\varphi /M_P)^n]`$, the first term being dominant Bento2 .

Let us examine to which extent the Grenoble experiment can allow imposing limits on the parameters of a putative new interaction of Nature once the energy resolution of the experiment is increased to $`\mathrm{\Delta }E=10^{18}eV`$.

In order to study possible limits on a fifth force which is sensitive to baryon number, we assume that its effects can be treated as a perturbation of the local gravitational potential. In this way, first order perturbation theory provides a correction to the energy of the nth-level defined by:

$$\mathrm{\Delta }E_n^{(1)}=_0^{\mathrm{}}\mathrm{\Psi }_n^{}V(x)\mathrm{\Psi }_n𝑑x,$$
(10)
where $`\mathrm{\Psi }_n`$ is the unperturbed eigenstate and the perturbation potential is given by $`V(x)=\alpha _5fmgx\mathrm{exp}(x/\lambda _5)`$, where $`f`$ is the fraction of the source (Earth) that contributes to the force acting on the neutrons. Since the neutrons are constrained vertically in the length scale $`x\mu m`$, the bounds on $`\alpha _5f`$ become interesting for say $`\lambda \text{ }>10\mu m`$. In this range, the exponential in $`V(x)`$ can be expanded around $`x=0`$. Retaining the first two terms in the expansion, we derive an upper bound for $`\alpha _5f`$ in terms of the energy resolution $`\mathrm{\Delta }E`$ of the experiment:

$$\alpha _5f\left[3.24+5.05\left(\frac{\lambda _5}{1\mu m}\right)^1\right]\text{ }<4.91\times 10^2\left(\frac{\mathrm{\Delta }E}{10^{18}eV}\right).$$
(11)
Here, we have used the numerical values for the following integrals

$$I_0=_{\alpha _0}^{\mathrm{}}\mathrm{\Phi }(z)^2𝑑z=1.24\times 10^4,$$
(12)
$$I_1=_{\alpha _0}^{\mathrm{}}\mathrm{\Phi }(z)z\mathrm{\Phi }(z)𝑑z=3.24\times 10^3,$$
(13)
$$I_2=_{\alpha _0}^{\mathrm{}}\mathrm{\Phi }(z)z^2\mathrm{\Phi }(z)𝑑z=8.64\times 10^4.$$
(14)
Directly from Eq. (11), we conclude that, within the range where no constraints exist $`10\mu m\text{ }<\lambda 1cm`$, the upper limit on $`\alpha _5f`$ would be of the order of $`1.3\times 10^2`$ as long as the energy resolution of $`\mathrm{\Delta }E=10^{18}eV`$ is achieved. Since the tested region for the range is quite small, $`f<<1`$ implying that there is no limit for $`\alpha _5`$ as the effect of a possible new force is far too small to affect neutrons in the experiment. If on the other hand the new force were long range so that $`f=1`$, still no bound on $`\alpha _5`$ would arise given the relationship between Newton’s constant, $`G_N`$ and $`G_{\mathrm{}}`$, namely $`G_N=G_{\mathrm{}}(1+\alpha _5)`$. We mention that the same conclusions can be drawn for the so-called massive version of the Brans-Dicke theory Acharya , where $`\lambda _5=m_\varphi ^1`$ and $`\alpha _5=1/2\omega +3`$, $`m_\varphi `$ being the scalar field mass and $`\omega `$ the Brans-Dicke coupling parameter.

We turn now to the issue of the gravitational screening. A possible absorption of the gravitational force between two bodies, when a medium is screened by another, has been the subject of a series of tests. This effect is reminiscent of the magnetic permeability of materials, and a screening or extinction coefficient, $`h`$, was proposed by Quirino Majorana Majorana  in 1920 in order to measure the ability of an object of dimension $`L`$ with density $`\rho (r)`$ to shield the gravitational force between masses $`m_1`$ and $`m_2`$:

$$F^{}=\frac{Gm_1m_2}{r^2}\mathrm{exp}\left[h_0^L\rho (r)𝑑r\right].$$
(15)
Naturally, $`h`$ must be quite small. Several attempts to measure this constant from general principles have been made. For instance, Weber Weber  has argued that quasi-static shielding could be predicted from a general relativistic analysis of tidal phenomena, and stated that the effect should be extremely small. More recently, it has been shown that a lunar laser ranging experiment can set the impressive limit, $`h1.0\times 10^{21}cm^2g^1`$ Eckhardt . The most stringent laboratory limit on the gravitational shielding constant has its origin on a recent measurement of Newton’s constant carried out at the Physik-Institut, Universität Zürich, which yields $`h4.3\times 10^{14}cm^2g^1`$ Unnikrishnan .

In order to use an improved version of the Grenoble experiment for obtaining a limit on the gravitational screening one considers the effect of the Moon on the local gravity acceleration. This effect is about $`g_M3\times 10^6g`$. Of course, an experiment of this nature implies accurate estimates of various local gravitational perturbations related to human activity and geophysical nature. Comparing the contribution of the Moon when Moon, Earth and the experiment are in alignment and the former is “screened” by Earth, then $`g_{M2}=3.3\times 10^6g`$, from which follows that an effect on the neutron spectra will be felt only if:

$$h>1.7\times 10^{10}cm^2g^1.$$
(16)
Here we have assumed that Earth’s density is constant over its diameter: $`\rho _{}=5.51gcm^3`$. The result in Eq. (16) is much less stringent than the abovementioned bounds. Similar analysis using the Sun leads to a limit that is about $`4`$ times greater.

Naturally, any actual experiment to measure gravitational origin effects on ultracold neutrons requires besides the six orders of magnitude improvement in the energy resolution of the neutron spectrum, a series of precautions in order to eliminate local gravitational perturbations and asymmetries.

###### Acknowledgements.

The authors would like to thank Clovis de Matos for discussions and Prof. Edward Witten for valuable comments.


text:
# Dynamical state of superclusters of galaxies: do superclusters expand or have they started to collapse?

## 1 INTRODUCTION

Superclusters of galaxies are the largest coherent and massive structures known in the Universe. The existence of superclusters is known since the pioneering studies of Shapley (1930). The nearest example is the Local Supercluster with the Virgo cluster as the central cluster (de Vaucouleurs 1956). The shape, mass and dynamical state of different superclusters have been studied in several papers (see e.g. Einasto et al. 1997; Ettori, Fabian & White 1997; Small et al. 1998; Barmby & Huchra 1998; Batuski et al. 1999; Bardelli et al. 2000; Basilakos, Plionis & Rowan-Robinson 2001; and references therein). Ettori, Fabian & White (1997) studied the mass distribution in the Shapley Supercluster, using X-ray observations. To investigate the dynamical state of the observed overdense regions they used the spherical collapse model (Gunn & Gott 1972; Peebles 1980). Ettori, Fabian & White (1997) found that the core region of the Shapley Supercluster with the radius $`6.7h^1`$Mpc is close to the turnaround point, when the perturbed region ceases to expand and begins to collapse. Bardelli et al. (2000) also analyzed the dynamical state of the Shapley Supercluster, using the spherical collapse model. They estimated that the total overdensity of galaxies in the Shapley Supercluster is ($`N/\overline{N}`$) $`11.3`$ on a scale of $`10.1h^1`$Mpc and concluded that, if light traces mass and the density parameter $`\mathrm{\Omega }=1`$, the Shapley Supercluster has already reached its turnaround radius and has started to collapse. Small et al. (1998) studied the structure and dynamics of the Corona Borealis Supercluster. They found that this supercluster may have started to collapse.

In this paper we study the dynamical state of superclusters in Lambda cold dark matter ($`\mathrm{\Lambda }`$CDM) cosmological models, using N-body simulations. These models successfully explain many observations of the large- and small-scale structure including the mass function and the peculiar velocities of clusters of galaxies (see e.g., Ostriker & Steinhardt 1995; Gramann & Hütsi 2001). We study the radial velocity and the density contrast in superclusters. The spherically averaged radial velocity around a system, in the shell of radius $`R`$, can be written as

$$u=HRv,$$
$`(1)`$
where $`v_H=HR`$ is the Hubble expansion velocity and $`v`$ is the averaged radial peculiar velocity toward the centre of the system. At the turnaround point, the peculiar velocity $`v=HR`$ and $`u=0`$. If $`v<HR`$, the system expands and if $`v>HR`$, the system begins to collapse. In the spherical collapse model, the peculiar velocity $`v`$ is directly related to the density contrast $`\delta `$ (see below). We study the relation between $`v`$ and $`\delta `$ in N-body simulations and compare the results with the spherical collapse model. Can we estimate the dynamical state of the supercluster on the basis of its density contrast?

The velocity field around superclusters and clusters in different N-body simulations has been studied in several papers (Lee, Hoffman & Ftaclas 1986; Villumsen & Davis 1986; Lilje & Lahav 1991; van Haarlem & van de Weygaert 1993; Hanski et al. 2001). Villumsen & Davis (1986) studied the nature of the velocity field around large clusters in $`\mathrm{\Omega }=1`$ CDM models. They found that the flowfields when averaged over $`4\pi `$ sr fit spherical collapse model relatively well for density contrasts $`\delta <3`$, but the mean radial peculiar velocity is systematically low for increasing $`\delta `$. van Haarlem & van de Weygaert (1994) studied the evolution of the velocity profile around clusters in $`\mathrm{\Omega }=1`$ CDM models. They found that for different clusters the agreement with the spherical collapse model can be different. There are clusters where the agreement is very good, but often the predictions of the spherical collapse model compare badly with the actual velocity field in N-body simulations (see Fig. 16-17 in their study). Hanski et al. (2001) studied the velocity profile around four clusters in a $`\mathrm{\Omega }_0=0.3`$ $`\mathrm{\Lambda }`$CDM model. They found that the radially averaged velocity fields around simulated clusters are compatible with the spherical collapse model.

In this paper we study the velocity field for superclusters in flat $`\mathrm{\Lambda }`$CDM models with the density parameter $`\mathrm{\Omega }_0=0.20.4`$, the baryon density $`\mathrm{\Omega }_bh^2=0.02`$ and the normalized Hubble constant $`h=0.7`$. These values are in agreement with measurements of the density parameter (e.g. Bahcall et al. 1999), with measurements of the baryon density from abundances of light elements (O’Meara et al. 2001; Tytler et al. 2000) and with measurements of the Hubble constant using various distance indicators (Freedman et al. 2001; see also Parodi et al. 2000). To restore the spatial flatness in the low-density models, we assume a contribution from the cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }=1\mathrm{\Omega }_0`$. We also assume that the initial density fluctuation field is a Gaussian density field. In this case, the power spectrum provides a complete statistical description of the field.

We define superclusters as maxima of the density field smoothed on the scale $`R=10h^1`$Mpc. Smaller superclusters in the density field smoothed on the scale $`R=5h^1`$Mpc are also investigated. We study the number density of superclusters, where the radially averaged peculiar velocity is larger than a given value, $`n(>v)`$. In particular, we investigate the number density of superclusters which have radially averaged peculiar velocities $`v>HR`$.

This paper is organized as follows. In Section 2 we consider the spherical collapse model in a flat universe with cosmological constant. We find the relation between the peculiar velocity and the density contrast for the values of the density parameter $`\mathrm{\Omega }_0=0.20.4`$. In Section 3 we describe the N-body simulations and the algorithms that have been used to identify the superclusters. In Section 4 we investigate the density contrast and the radial peculiar velocity for the superclusters in the N-body simulations and compare the results with the spherical collapse model. The function $`n(>v)`$ is analyzed in Section 5 and a summary is presented in Section 6.

## 2 THE SPHERICAL COLLAPSE MODEL

The spherical collapse model describes the evolution of a spherically symmetric perturbation in an expanding universe. Under the assumption of sphericity, the nonlinear dynamics of a collapsing shell is determined by the mass interior to it. The spherical collapse model has been discussed in detail by Tolman (1934), Bondi (1947), Gunn & Gott (1972), Silk (1974, 1977), Peebles (1980) and Schechter (1980). This model has been used in the Press-Schechter (Press & Schechter 1974) formalism to evaluate the mass function of collapsed objects in the universe, as well as to determine the density parameter of the universe from the infalling flow (e.g. Davis & Huchra 1982). In recent years, the spherical collapse model has been used to study the dynamical state of superclusters of galaxies (Ettori, Fabian & White 1997; Small et al. 1998; Bardelli et al. 2000).

Although the Zel’dovich approximation (Zel’dovich 1970) generically yields pancakes and numerical simulations show that many collapsing regions are filamentary (e.g. Gramann 1988; Bertchinger & Gelb 1991), the spherical collapse model has been popular because of its simplicity. Despite its idealized nature, the spherical model appears to give a reasonable description of collapse of high peaks in a random Gaussian density field (Bernardeau 1994). Eisenstein & Loeb (1995) studied an analytical model for the triaxial collapse of cosmological perturbations. They found that the turn-around times of the short axes of a triaxial perturbation are fairly well predicted by the turn-around time of a spherical perturbation with the same initial density.

We consider the spherical collapse model in a flat universe with a cosmological constant, $`\mathrm{\Lambda }`$. This model was studied by Peebles (1984), Lahav et al. (1991), Eke, Cole & Frenk (1996) and Lokas & Hoffman (2001). In the spherical collapse model, the mass shell with a radius $`R(t)`$ contains a fixed mass $`m`$ and its dynamics satisfies the energy equation

$$\left(\frac{dR}{dt}\right)^2=\frac{2Gm}{R}+\frac{\mathrm{\Lambda }R^2}{3}K.$$
$`(2)`$
The constant $`K`$ is positive for a growing mass perturbation. In a flat universe, the evolution of a scale factor is described by

$$\left(\frac{1}{a}\frac{da}{dt}\right)^2=\frac{8}{3}\pi G\rho _b+\frac{\mathrm{\Lambda }}{3},$$
$`(3)`$
where $`\rho _b`$ is the mean background density. Using eq. (3) we can write eq. (2) as (Peebles 1984)

$$\left(\frac{dR}{da}\right)^2=\frac{a}{R}\frac{(\omega R^3\kappa R+1)}{(\omega a^3+1)},$$
$`(4)`$
where $`\kappa `$ is a new constant, the present value of $`a`$ is $`a_0=1`$ and

$$\omega =\frac{1}{\mathrm{\Omega }_0}1.$$
$`(5)`$
The density contrast can be expressed as

$$\delta =\frac{a^3}{R^3}1.$$
$`(6)`$
The radial velocity in units of the Hubble velocity is

$$\frac{u}{HR}=\frac{a}{R}\frac{dR}{da}=1\frac{v}{HR},$$
$`(7)`$
where $`v`$ is the peculiar velocity toward the centre of the system.

The perturbed region will turn around when

$$\omega R^3\kappa R+1=0.$$
$`(8)`$
Solving this cubic equation and requiring that a physically sensible root exists, we find that the smallest perturbation which will collapse has

$$\kappa _{\mathrm{min}}=\frac{3\omega ^{\frac{1}{3}}}{2^{\frac{2}{3}}}.$$
$`(9)`$
The corresponding turnaround radius is

$$R_{\mathrm{ta},\mathrm{max}}=(2\omega )^{\frac{1}{3}}.$$
$`(10)`$
Perturbations with larger value of $`\kappa `$ turn around and collapse earlier. Solving the cubic eq. (8) we can find the turnaround radius $`R_{ta}`$ as a function of the density parameter $`\omega `$ and constant $`\kappa `$. The solution can be expressed as (Eke, Cole & Frenk 1996)

$$R_{\mathrm{ta}}=2s^{\frac{1}{3}}\mathrm{cos}\left(\frac{\theta }{3}\frac{2\pi }{3}\right),$$
$`(11)`$
where

$$s=\left(\frac{3}{2}\right)^3\left(\frac{\kappa }{\kappa _{\mathrm{min}}^3}\right)^{\frac{3}{2}}$$
$`(12)`$
and $`\theta `$ satisfies

$$\mathrm{cos}\theta =\left(\frac{\kappa _{\mathrm{min}}}{\kappa }\right)^{\frac{3}{2}}(0<\theta <\frac{\pi }{2}).$$
$`(13)`$
In the limit where $`a1`$, we find that

$$\left(\frac{a}{R}\frac{dR}{da}\right)^2=\frac{a^3}{R^3}(1\kappa R)=1+\delta \kappa a$$
$`(14)`$
and therefore, the density contrast can be written as

$$\delta =2\left(\frac{v}{HR}\right)+\kappa a.$$
$`(15)`$
In linear perturbation theory the peculiar velocity for a spherical mass concentration is determined as

$$\frac{v}{HR}=\frac{1}{3}f(\mathrm{\Omega })\delta ,$$
$`(16)`$
where $`f(\mathrm{\Omega })`$ is the dimensionless growth factor. For $`a1`$, $`f(\mathrm{\Omega })=1`$. By using the relation (16), we find from eq.(15) that

$$\delta =\frac{3\kappa }{5}a.$$
$`(17)`$
In the limit where $`a1`$, the density contrast is directly related to the constant $`\kappa `$.

Fig.1 shows the relation between the peculiar velocity in units of the Hubble velocity, $`v/HR`$, and the density contrast $`\delta `$ for the flat models at the present moment. We integrated eq. (4) and found the present density contrast and peculiar velocity from eq. (6) and (7), respectively. For each value of the initial density contrast $`\delta _i`$ at $`a_i1`$, the initial radius $`R_i`$ was calculated using eq. (6) and the constant $`\kappa `$ was determined using eq. (17). For comparison, we show the peculiar velocity given by the approximation derived by Yahil (1985),

$$\frac{v}{HR}=\frac{1}{3}f(\mathrm{\Omega }_0)\delta (1+\delta )^{\frac{1}{4}}.$$
$`(18)`$
The growth factor $`f(\mathrm{\Omega }_0)`$ is $`0.41`$, $`0.51`$ and $`0.60`$ for the flat models with $`\mathrm{\Omega }_0=0.2,0.3`$ and $`0.4`$, respectively. Fig. 1 shows that the approximation (18) agrees well with the numerical solution for the density contrast $`\delta <10`$. For larger values of $`\delta `$, this approximation overestimates the peculiar velocity.

Consider the density contrast at the turnaround point, $`\delta _{\mathrm{ta}}`$, where $`v=HR`$. We found that $`\delta _{\mathrm{ta}}=15.48`$, $`11.19`$ and $`8.95`$ in the flat models with $`\mathrm{\Omega }_0=0.2`$, $`0.3`$ and $`0.4`$, respectively. For comparison, we calculated also the turnaround density contrast for the corresponding open models. The turnaround density contrast in the open models can be expressed as (Regös & Geller 1989)

$$\delta _{\mathrm{ta}}=\left(\frac{\pi }{x_1\mathrm{ln}x_2}\right)^2\omega ^31,$$
$`(19)`$
where $`x_1^2=(2\mathrm{\Omega }_0^11)^21`$ and $`x_2=x_1+2\mathrm{\Omega }_0^11`$. We found that $`\delta _{\mathrm{ta}}=16.22`$, $`11.57`$ and $`9.17`$ for the open models with $`\mathrm{\Omega }_0=0.2`$, $`0.3`$ an $`0.4`$, respectively. For the same value of $`\mathrm{\Omega }_0`$, the turnaround density in the open model is somewhat higher than in the flat model. However, the effect of cosmological constant on $`\delta _{\mathrm{ta}}`$ is small.

Next we study the relation between the radial peculiar velocity and density contrast for superclusters in N-body simulations.

## 3 NUMERICAL MODELS

### 3.1 Simulations

First, we studied superclusters in a N-body simulation carried out by the Virgo consortium for the flat $`\mathrm{\Lambda }`$CDM model with a cosmological constant. The Virgo simulations are described in detail by Jenkins et al. (1998). These simulations were created using an adaptive particle-particle/particle-mesh (AP3M) code as described by Couchman, Thomas & Pearce (1995) and Pearce & Couchman (1997). In the $`\mathrm{\Omega }_0=0.3`$ $`\mathrm{\Lambda }`$CDM model studied here, the power spectrum of the initial conditions was chosen to be in the form given by Bond & Efstathiou (1984),

$$P(k)=\frac{Ak}{[1+(aq+(bq)^{3/2}+(cq)^2)^\nu ]^{2/\nu }},$$
$`(20)`$
where $`q=k/\mathrm{\Gamma }`$, $`a=6.4h^1`$Mpc, $`b=3h^1`$Mpc, $`c=1.7h^1`$Mpc, $`\nu =1.13`$ and $`\mathrm{\Gamma }=\mathrm{\Omega }_0h=0.21`$. The normalization constant, $`A`$, was chosen by fixing the value of $`\sigma _8`$ (the linearly extrapolated mass fluctuation in spheres of radius $`8h^1`$Mpc) to be $`0.9`$. The evolution of particles was followed in the comoving box of size $`L=239.5h^1`$Mpc. The number of particles was $`N_p=256^3`$. Therefore, the mean particle separation $`\lambda _p=L/N_p^{1/3}=0.9355h^1`$Mpc. We denote this $`\mathrm{\Lambda }`$CDM model as the model M1.

We also investigated superclusters in a series of N-body simulations we carried out for flat $`\mathrm{\Lambda }`$CDM models, where the density parameter $`\mathrm{\Omega }_0=0.20.4`$ and $`\sigma _8=0.70.9`$. In these simulations, the evolution of particles was followed by using a particle-mesh (PM) code described by Gramann (1988) and Suhhonenko & Gramann (1999). The PM method is discussed in detail by Hockney & Eastwood (1981) and Efstathiou et al. (1985). We used the PM code, where the number of cells was the same as the number of particles. We denote these flat $`\mathrm{\Lambda }`$CDM models as the models M2-M3, S1-S4 and L1-L4.

In the models M2 and M3, the cosmological parameters were chosen similar to the model M1. We chose the flat $`\mathrm{\Omega }_0=0.3`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$ model and used the initial power spectrum $`P(k)`$ given in eq. (20) (with $`\sigma _8=0.9`$). In the model M2, we followed the evolution of $`256^3`$ particles in the box of size $`L=239.5h^1`$Mpc. In the model M3, the number of particles was $`N_p=128^3`$ and the box size $`L=192h^1`$Mpc. Table 1 lists the code, the number of particles and the box size used in models M1, M2 and M3, respectively. Suhhonenko & Gramann (2002) used the models M1 and M2 to investigate the rms peculiar velocity of galaxy clusters for different cluster masses and radii.

In the models S1-S4 and L1-L4, the number of particles was $`N_p=128^3`$. In these models, the initial power spectrum was calculated by using the fast Boltzmann code CMBFAST developed by Seljak & Zaldarriaga (1996). The normalized Hubble constant was chosen to be $`h=0.7`$ and the baryon density $`\mathrm{\Omega }_bh^2=0.02`$. Table 2 lists the density parameter $`\mathrm{\Omega }_0`$, $`\sigma _8`$ and the box size $`L`$ used in the models S1-S4 and L1-L4. We studied the flat models with a cosmological constant $`\mathrm{\Omega }_\mathrm{\Lambda }=1\mathrm{\Omega }_0`$.

For the models M3 and S3, all the parameters are the same, except of the form of the initial power spectrum. In the model M3, we used the approximation given in eq. (20) and in the model S3, the power spectrum was calculated numerically by using the code CMBFAST.

Fig. 2 shows the power spectra used in the N-body simulations for different values of the parameters $`\mathrm{\Omega }_0`$ and $`\sigma _8`$. The light curves show the numerical power spectra calculated by using the code CMBFAST. In the models studied, the power spectrum on large scales is quite different. For the same value of $`\sigma _8`$, the amplitude of the large-scale density fluctuations is highest for $`\mathrm{\Omega }_0=0.2`$ and lowest for $`\mathrm{\Omega }_0=0.4`$. The light solid line shows the power spectrum for $`\mathrm{\Omega }_0=0.3`$ and $`\sigma _8=0.9`$. For comparison, we show the power spectrum given in eq. (20) for $`\mathrm{\Omega }_0=0.3`$ and $`\sigma _8=0.9`$ (heavy solid line). We see that on large scales, the numerical power spectrum is considerably higher than the power spectrum given in eq. (20) ($`26`$ % at the wavelength $`\lambda =2\pi /k=192h^1`$Mpc).

### 3.2 Selection of superclusters

To select the superclusters in the N-body simulations, we used the following method:

(1) We calculated the density contrast on a grid. For each grid point, the density contrast was determined as

$$\delta =\frac{N}{\overline{N}}1,$$
$`(21)`$
where $`N`$ is the number of particles in the sphere of radius $`R`$ around the grid point, and

$$\overline{N}=\frac{4\pi }{3}\frac{N_pR^3}{L^3}$$
$`(22)`$
is the mean number of particles in the sphere of radius $`R`$. We denote the number of cells in the grid as $`N_g^3`$.

(2) We found the density maxima on the grid. The grid point was considered as a density maximum, if its density contrast was higher than the density contrast in all neighbouring grid points. In the following, we studied the density maxima, where the density contrast was $`\delta >1`$.

(3) We identified the candidate supercluster centres. For each grid point, where the density maximum was determined, we found all particles in a cube with the size $`l=L/N_g`$ around this grid point. For each particle we calculated the density contrast using eq. (21), where $`N`$ was the number of particles in the sphere of radius $`R`$ around this fixed particle. The location of the particle, where the density contrast had the maximum value, was identified as the candidate supercluster centre.

(4) The final supercluster list was obtained by deleting the candidate superclusters with lower density contrast in all pairs separated by less than the radius $`R`$.

In this way we define the superclusters as maxima of the density field smoothed with a top-hat window with radius $`R`$. In the models M1-M3 and S1-S4, we used the smoothing radius $`R=5h^1`$Mpc. In the models M1 and M2, the mean number of particles in spheres of $`R=5h^1`$Mpc is $`\overline{N}=639.4`$, in the models M3 and S1-S4 - $`\overline{N}=155.1`$. To select superclusters for the models M1 and M2, we used a $`80^3`$ grid. In this case, the cell size was $`l=L/N_g=2.99h^1`$Mpc. To identify superclusters for the models M3 and S1-S4, we used a $`64^3`$ grid ($`l=3h^1`$Mpc). In the models L1-L4, we studied the superclusters defined with $`R=10h^1`$Mpc. In these models, the mean number of particles in spheres of $`R=10h^1`$Mpc is $`\overline{N}=155.1`$. To select superclusters, we used a $`64^3`$ grid ($`l=6h^1`$Mpc).

## 4 RADIAL VELOCITY AND DENSITY CONTRAST IN THE SUPERCLUSTERS

Table 1 and Table 2 show the number of superclusters, $`N_s`$, identified in different models. For each supercluster, we investigated the density contrast, $`\delta `$, and the radially averaged peculiar velocity, $`v`$. The density contrast in a supercluster was determined as

$$\delta =\frac{N_f}{\overline{N}}1,$$
$`(23)`$
where $`N_f`$ is the number of particles in the sphere of radius $`R`$ around the centre of the supercluster. The radially averaged peculiar velocity for the supercluster, in a shell of radius $`R`$, was defined as

$$v=\frac{1}{N_l}\underset{i=1}{\overset{N_l}{}}(\stackrel{}{v}_i^{}\stackrel{}{v}_0^{})\widehat{\stackrel{}{r}_i},$$
$`(24)`$
where $`N_l`$ is the number of particles in the shell $`0.1R`$ wide, $`\stackrel{}{v}_i^{}`$ is the peculiar velocity of the particle $`i`$ at a radial distance $`\stackrel{}{r}_i`$ from the centre of the supercluster, $`\widehat{\stackrel{}{r}_i}`$ is the unit vector in the direction of the particle $`i`$ and $`\stackrel{}{v}_0^{}`$ is the mean peculiar velocity of all the particles within the sphere of radius $`R`$,

$$\stackrel{}{v}_0^{}=\frac{1}{N_f}\underset{i=1}{\overset{N_f}{}}\stackrel{}{v}_i^{}.$$
$`(25)`$
Fig. 3 shows the averaged radial peculiar velocity of a supercluster in units of the Hubble velocity, $`v/HR`$, vs. the density contrast of the supercluster in the models M1, M2 and M3. The superclusters were defined with $`R=5h^1`$Mpc. In the models M1, M2 and M3, the number of superclusters is $`1414`$, $`1374`$ and $`662`$, respectively. We see a close correlation between the peculiar velocity and the density contrast of a supercluster. We can estimate the dynamical state of the supercluster on the basis of its density contrast.

For comparison, we also show in Fig. 3 the peculiar velocity predicted in the flat spherical collapse model for $`\mathrm{\Omega }_0=0.3`$. The radially averaged peculiar velocity of superclusters is systematically smaller than that predicted by the spherical collapse model.

We investigated the mean peculiar velocity and the standard deviation for the superclusters in different density intervals. The superclusters were divided into subgroups according to their density contrast. We studied the velocities in fifteen subgroups, where the logarithm of the density contrast was in the range $`\mathrm{log}\delta =00.1`$, $`0.10.2`$, …, $`1.41.5`$. The mean peculiar velocity and standard deviation were determined for the intervals which contained at least eight superclusters. The lower right panel in Fig. 3 shows the mean peculiar velocity determined in the models M1, M2 and M3. The error bars show the standard deviation. For clarity, the standard deviation is only shown for the model M3. In the models M1, M2 and M3, the standard deviation is similar.

The models M1 and M2 differ from each other by the code used to follow the evolution of particles. The model M1 was created using a AP3M code and the model M2 by using a PM code. Fig. 3 shows that in the model M1, the mean peculiar velocities are somewhat higher than in the model M2. For example, in the range $`\delta =10.012.6`$, the mean peculiar velocities are $`\overline{v}=0.78HR`$ and $`\overline{v}=0.68HR`$ in the models M1 and M2, respectively. This difference is probably due to the smoothing inherent to the PM method. However, the effect of the code on the mean peculiar velocity is not large ($`15`$% in the range $`\delta =10.012.6`$).

In the model M3, the box size $`L`$ and the number of particles, $`N_p`$, are smaller than in the models M1 and M2. In this model for $`R=5h^1`$Mpc, the box size $`L=38.4R`$ and the mean particle separation $`\lambda _p=0.3R`$. The same values for $`L/R`$ and $`\lambda _p/R`$ are used in the models S1-S4 and L1-L4. Fig.3 shows that the mean peculiar velocities in the model M3 are similar to the mean velocities in the models M1 and M2. For example, in the range $`\delta =1012.6`$, the mean peculiar velocity is $`\overline{v}=0.74HR`$. Therefore, a box size $`L=38.4R`$ and a mean particle separation $`\lambda _p=0.3R`$ are sufficient to study the mean peculiar velocities in superclusters defined with radius $`R`$.

We can compare the mean peculiar velocities of superclusters with the predictions of the spherical collapse model. Lets consider the model M1. In this model, for the range $`\delta =1.001.26`$, the mean density contrast is $`1.14`$ and the mean peculiar velocity is $`\overline{v}=0.13HR`$. For comparison, in the spherical collapse model, the peculiar velocity for $`\delta =1.14`$ is $`v=0.17HR`$. For the range $`\delta =10.012.6`$, the mean density contrast for superclusters is $`11.2`$ and $`\overline{v}=0.78HR`$. In the spherical collapse model, the peculiar velocity for $`\delta =11.2`$ is $`v=1.0HR`$. Therefore, the mean peculiar velocities in the simulation are $`25\%`$ smaller than those predicted by the spherical collapse model. The difference in the peculiar velocities for $`\delta 1`$ and $`\delta 10`$ is similar. In the model M2 and M3, the deviations from the spherical collapse model are somewhat larger than in the model M1 ($`35`$% and $`30`$% for the model M2 and M3, respectively).

Fig. 4 shows the radially averaged peculiar velocities for the superclusters in the models S1-S4. The mean peculiar velocity with a standard deviation is determined similarly as in the models M1-M3. Panel (a) shows the peculiar velocities for the model S1 ($`\mathrm{\Omega }_0=0.2`$) and S4 ($`\mathrm{\Omega }_0=0.4`$). For the same value of the density contrast, the mean peculiar velocity is larger, if the density parameter $`\mathrm{\Omega }_0`$ is larger. For the range $`\delta =5.06.3`$, the peculiar velocity $`v/HR=0.37\pm 0.08`$ and $`v/HR=0.52\pm 0.11`$ in the model S1 and S4, respectively. For the range $`\delta =1012.6`$, the peculiar velocity $`v/HR=0.62\pm 0.08`$ and $`v/HR=0.87\pm 0.10`$, respectively.

Panel (b) in Fig. 4 shows the peculiar velocities for the superclusters in the model S2 and S3. In these models $`\mathrm{\Omega }_0=0.3`$. In the model S2, the parameter $`\sigma _8=0.7`$ and in the model S3, the $`\sigma _8=0.9`$. We see that the mean peculiar velocities in the model S3 are somewhat smaller than the peculiar velocities in the model S2. The difference between the models S2 and S3 is not large. For the range $`\delta =5.06.3`$, the peculiar velocity $`v/HR=0.48\pm 0.07`$ and $`v/HR=0.43\pm 0.10`$ in the model S2 and S3, respectively. For the range $`\delta =1012.6`$, the peculiar velocity $`v/HR=0.77\pm 0.14`$ and $`v/HR=0.73\pm 0.12`$, respectively. The peculiar velocities in the model S3 are similar to the peculiar velocities in the model M3.

We compared the mean peculiar velocities for the superclusters in the models S1-S4 with the predictions of the spherical collapse model. The mean peculiar velocities in the models S1, S2 and S4 are $`25`$% smaller than those predicted by the spherical collapse model. In the model S3, the difference is $`30`$%. In all models, this difference is similar for different values of the density contrast.

We also determined the turnaround density contrast $`\delta _{\mathrm{ta}}`$ for the $`R=5h^1`$Mpc superclusters. We found all the superclusters, where the peculiar velocity in units of the Hubble velocity, $`v/HR`$, is in the range \[0.9, 1.1\] and calculated the mean density contrast and the standard deviation for these superclusters. Table 3 shows the density contrast $`\delta _{\mathrm{ta}}`$ for the models studied. For the model M1, we found that $`\delta _{\mathrm{ta}}=16.2\pm 4.6`$. For comparison, in the spherical collapse model the turnaround density contrast is $`11.19`$ for $`\mathrm{\Omega }_0=0.3`$ (see Section 2). For the model S1 and S4, we found that $`\delta _{\mathrm{ta}}=23.7\pm 7.1`$ and $`\delta _{\mathrm{ta}}=12.3\pm 2.2`$, respectively. In the spherical model, $`\delta _{ta}=15.48`$ and $`8.95`$ for $`\mathrm{\Omega }_0=0.2`$ and $`\mathrm{\Omega }_0=0.4`$, respectively. The turnaround density contrast for the superclusters in N-body simulations is substantially larger than the turnaround density contrast in the spherical collapse model.

We also studied superclusters in the models L1-L4. In these models the box size was $`L=384h^1`$Mpc and we investigated the superclusters in the density field smoothed on the scale $`R=10h^1`$Mpc. Fig. 5 shows the radially averaged peculiar velocity of a supercluster in units of the Hubble velocity, $`v/HR`$, vs. the density contrast of the supercluster, for 687 superclusters of the model L3. For comparison, we show the peculiar velocities predicted in the spherical collapse model for $`\mathrm{\Omega }_0=0.3`$. We see that also for the superclusters studied with $`10h^1`$Mpc smoothing, the mean infall velocities are systematically smaller than those predicted by the spherical collapse model.

Fig. 6 shows the peculiar velocities for the superclusters in the models L1-L4. For comparison, the dotted line in panel (b) shows the peculiar velocities for the $`R=5h^1`$Mpc superclusters in the model M1. For the $`R=10h^1`$Mpc superclusters, the mean peculiar velocities in units of the Hubble velocity, $`\overline{v}/HR`$, are somewhat higher than for the $`R=5h^1`$Mpc superclusters. For the range $`\delta =5.06.3`$, the peculiar velocity $`v/HR=0.43\pm 0.04`$ and $`v/HR=0.63\pm 0.06`$ in the model L1 and L4, respectively. The peculiar velocity $`v/HR=0.54\pm 0.04`$ and $`v/HR=0.51\pm 0.05`$, in the model L2 and L3, respectively. The mean peculiar velocities for the superclusters in the models L1-L4 are $`15`$% smaller than those predicted by the spherical collapse model.

## 5 VELOCITY DISTRIBUTION FOR THE SUPERCLUSTERS

We investigated the number of superclusters, where the radially averaged peculiar velocity is larger than a given value, $`N(>v)`$. We denote the number of superclusters, which have peculiar velocities $`v>HR`$ and, therefore, have already started to collapse, as $`N_c`$. We also examined the number of superclusters, which have peculiar velocities $`v>0.5HR`$. We denote this number as $`N_{0.5}`$. Table 1 and Table 2 show the numbers $`N_c`$ and $`N_{0.5}`$ for the different models studied. In the model M1, we found $`63`$ superclusters which have already started to collapse.

Panel (a) in Fig. 7 shows the number density of superclusters, where the radially averaged peculiar velocity is larger than a given value, $`n(>v)`$. The number density is given for the superclusters in the models M1-M3. The function $`n(>v)`$ was determined as

$$n(>v)=\frac{N(>v)}{L^3}.$$
$`(26)`$
We also show the Poisson error bars for the number densities. Panel (b) in Fig. 7 demonstrates the density distribution for the superclusters in the models M1-M3. The density distribution in these models is similar, but the velocity distribution is different. The function $`n(>v)`$ for the model M1 is higher than for the model M2. In the model M1 and M2, the number density of superclusters that have peculiar velocities $`v>HR`$, is $`4.6\times 10^6h^3`$Mpc-3 and $`1.7\times 10^6h^3`$Mpc-3, respectively. The difference between the models M1 and M2 is smaller for smaller peculiar velocities. The velocity distribution in the model M3 is similar to the model M2. Therefore, the number density $`n(>v)`$ in the model M3 is probably underestimated.

Fig. 8 shows the number density, $`n(>v)`$, for the superclusters in the models S1-S4 and L1-L4. We can use models S1-S4 and L1-L4 to estimate the lower limit for the function $`n(>v)`$ in different models. Also, we can examine the relative differences between the different models. Panel (a) in Fig. 8 shows the number density for the $`R=5h^1`$Mpc superclusters in the models S1-S4 ($`L=192h^1`$Mpc). The function $`n(>v)`$ is lowest for the model S1 ($`\mathrm{\Omega }_0=0.2`$) and highest for the model S4 ($`\mathrm{\Omega }_0=0.4`$). In the model S1, we found $`4`$ superclusters, which have started to collapse. In the model S4, this number was $`25`$. The number density of superclusters, which have peculiar velocities $`v>0.5HR`$, is $`1.1\times 10^5h^3`$Mpc-3 and $`2.7\times 10^5h^3`$Mpc-3, in the model S1 and S4, respectively.

Panel (b) in Fig. 8 shows the number density for the $`10h^1`$Mpc superclusters in the models L1-L4 ($`L=384h^1`$Mpc). The velocity distributions in the models L1 and L2 are similar. In the model L1, the density parameter $`\mathrm{\Omega }_0`$ is smaller, but the amplitude of the large-scale density fluctuations is larger than in model L2 (see Fig. 2). The number density of superclusters, which have peculiar velocities $`v>0.5HR`$, is $`3.5\times 10^7h^3`$Mpc-3 and $`2.4\times 10^7h^3`$Mpc-3, in the models L1 and L2, respectively. The velocity distributions in the models L3 and L4 are also similar. In the model L4, the density parameter $`\mathrm{\Omega }_0`$ is higher, but the amplitude of the large-scale density fluctuations is smaller than in the model L3. The number density of superclusters, which have peculiar velocities $`v>0.5HR`$, is $`1.0\times 10^6h^3`$Mpc-3 and $`0.9\times 10^6h^3`$Mpc-3, in the models L3 and L4, respectively.

In the simulations studied, all the superclusters defined with $`10h^1`$Mpc smoothing are expanding by the present epoch. In the models L3 ($`\mathrm{\Omega }_0=0.3`$) and L4 ($`\mathrm{\Omega }_0=0.4`$), the maximum density contrast for the $`R=10h^1`$Mpc superclusters is $`11.2`$ and $`7.9`$, respectively. The turnaround density contrast for the $`R=10h^1`$Mpc superclusters is expected to be similar or somewhat smaller than for the $`R=5h^1`$Mpc superclusters (see Fig. 4 and Fig. 6). In the model M1 ($`\mathrm{\Omega }_0=0.3`$) and S4 ($`\mathrm{\Omega }_0=0.4`$), the turnaround density contrast for the $`R=5h^1`$Mpc superclusters is $`\delta _{ta}=16.2\pm 4.6`$ and $`\delta _{ta}=12.3\pm 2.2`$, respectively. Therefore, the $`R=10h^1`$Mpc superclusters in the models L3 and L4 are close to the turnaround point.

## 6 SUMMARY

In this paper we have investigated the dynamical state of superclusters in flat $`\mathrm{\Lambda }`$CDM models. We analyzed superclusters in the Virgo simulation for the $`\mathrm{\Lambda }`$CDM model with $`\mathrm{\Omega }_0=0.3`$ and $`\sigma _8=0.9`$. This simulation was carried out using the AP3M code. We also investigated superclusters in a series of N-body simulations we carried out for flat $`\mathrm{\Lambda }`$CDM models, where the density parameter $`\mathrm{\Omega }_0=0.20.4`$ and $`\sigma _8=0.70.9`$. In these simulations, a PM code was used. To identify superclusters in simulations we used a method where superclusters were defined as maxima of the density field smoothed on the scale $`R=10h^1`$Mpc. Smaller superclusters in the density field smoothed on the scale $`R=5h^1`$Mpc were also examined. For each supercluster identified, we determined the density contrast, $`\delta `$, and the radially averaged peculiar velocity, $`v`$ (see eq. (23-24)).

The results for the superclusters in the N-body simulations were compared with the spherical collapse model. We found that the radial peculiar velocities in the N-body simulations are systematically smaller than those predicted by the spherical collapse model ($`15`$% and $`25`$% for the $`R=10h^1`$Mpc and $`R=5h^1`$Mpc superclusters, respectively). The deviations from the spherical collapse model are similar for $`\delta 1`$ and $`\delta 10`$. We also studied the turnaround density contrast for superclusters. We found that the turnaround density contrast for superclusters is substantially larger than that predicted by the spherical collapse model. In the Virgo simulation with $`\mathrm{\Omega }_0=0.3`$, the turnaround density contrast is $`\delta _{\mathrm{ta}}=16.2\pm 4.6`$ for the $`R=5h^1`$Mpc superclusters. For comparison, in the spherical model the turnaround density contrast is $`11.19`$ for $`\mathrm{\Omega }_0=0.3`$.

We found the relations between $`\delta `$ and $`v`$ for different cosmological models. These relations can be used to estimate the dynamical state of a supercluster on the basis of its density contrast. The density contrast for superclusters on large scales can be determined from the new large galaxy redshift surveys. The results of our study show how the density contrast is related to the radial peculiar velocity of the supercluster. For example, in the model with $`\mathrm{\Omega }_0=0.3`$, $`\sigma _8=0.9`$ and $`R=10h^1`$Mpc, the peculiar velocity $`v/HR=0.54\pm 0.04`$ for the range $`\delta =5.06.3`$. The $`v/HR\delta `$ relation for different values of $`R`$ and $`\sigma _8`$ is similar.

On the other hand, recent improvements of the techniques for measuring distances of galaxies allow compilation of large samples of galaxies with measured peculiar velocities. We can use this information to study the radially averaged peculiar velocity of superclusters. On the basis of this peculiar velocity we can estimate the mass-density contrast in the superclusters and compare it with the galaxy-density contrast. In this way we can study the relation between the mass distribution and the galaxy distribution on large scales.

In the simulations studied, all the superclusters defined with $`10h^1`$Mpc smoothing have peculiar velocities $`v<HR`$ and, therefore, are expanding by the present epoch. A small fraction of the $`R=5h^1`$Mpc superclusters has already reached their turnaround radius and these superclusters have started to collapse. In the Virgo simulation, the number density of superclusters, which have started to collapse, is $`4.6\times 10^6h^3`$Mpc-3. The $`R=10h^1`$Mpc superclusters in the models with $`\mathrm{\Omega }_0=0.3`$, $`\sigma _8=0.9`$ and $`\mathrm{\Omega }_0=0.4`$, $`\sigma _8=0.8`$ are close to the turnaround point.

## ACKNOWLEDGEMENTS

We thank J. Einasto, M. Einasto, P. Heinamäki, G. Hütsi, J. Jaaniste and E. Saar for useful discussions. This work has been supported by the ESF grant 3601. The Virgo simulation used in this paper is available at http://www.mpa-garching.mpg.de/Virgo/virgoproject.html. This simulation was carried out at the Computer Center of the Max-Planck Society in Garching and at the EPCC in Edinburgh, as part of the Virgo Consortium project. The code CMBFAST used in this paper, is available at http://arcturus.mit.edu/ matiasz/CMBFAST/cmbfast.html.