text:
# Ab initio Monte Carlo simulations for finite-temperature properties: Application to lithium clusters and bulk liquid lithium

## I  INTRODUCTION

Molecular dynamics (MD) and Monte Carlo (MC) simulations have been the major techniques for calculating finite-temperature properties of condensed-matter systems. In particular, ab initio MD simulations, which combine classical molecular dynamics with first-principles quantum theory, determine the atomic and electronic structures of a system simultaneously and accurately, as well as the temperature-dependent properties and time evolution of the system. Car and Parrinello (1985); Payne et al. (1992) Ab initio MD methods have been widely used to provide reliable and accurate information for various systems since Car and Parrinello first introduced their work in 1985. Car and Parrinello (1985)

On the other hand, ab initio Monte Carlo simulations, which employ first-principles quantum theory to calculate the potential energy of a system at each classical MC step, have so far been used only to a very limited extent. The majority of MC simulations for condensed-matter systems have been performed with empirical or semiempirical atom-atom interaction potentials, which are fitted to available experimental data or ab initio calculations. While the atom-atom potentials may work well in some cases, they do not provide reliable results in other instances. Since the quantum many-body effects and the detailed electronic structures cannot be included in empirical or semi-empirical potentials, first-principles theory would obviously be desirable for providing reliable energetics. The major concern with ab initio Monte Carlo simulations is the high computational demand. In particular, Monte Carlo simulations with empirical potentials often generate millions of sampling configurations. It is a formidable task to determine the potential energies of such a huge number of configurations with ab initio total-energy calculations. In recent years, however, several groups have performed ab initio MC simulations for very small clusters, and they have tried to reduce the number of sampling configurations (to typically $`10^4`$). Keshari and Ishikawa (1994); Truong and Stefanovich (1996); Asada and Iwata (1996); Akhmatskaya et al. (1997); Estrin et al. (1997); Weht et al. (1998); Jellinek et al. (1998); Watanabe and Asada (1998); Bandyopadhyay et al. (1999); Bacelo et al. (1999); Srinivas and Jellinek (2000); Ishikawa et al. (2001) Ishikawa et al Keshari and Ishikawa (1994); Bacelo et al. (1999); Ishikawa et al. (2001) applied ab initio calculations to the Monte Carlo simulated annealing algorithm Kirkpatrick et al. (1983). They performed simulations for the Li6 cluster, small hydrogenated lithium clusters (Li5H, Li4H2, and Li7H), and HCl(H2O)n ($`n=3,4`$) clusters. Keshari and Ishikawa (1994); Bacelo et al. (1999); Ishikawa et al. (2001) Other groups Truong and Stefanovich (1996); Asada and Iwata (1996); Estrin et al. (1997); Watanabe and Asada (1998) used the ab initio molecular orbital (MO) method in classical MC simulations to study the interactions of ions and molecules with small water clusters. Structural and thermal properties of water dimers Bandyopadhyay et al. (1999) and lithium clusters (Li4, Li$`{}_{}{}^{+}{}_{5}{}^{}`$, and Li8) Weht et al. (1998); Jellinek et al. (1998); Srinivas and Jellinek (2000) were also investigated with ab initio MC simulations.

We emphasize that all the previous ab initio MC simulations were limited to very small clusters. Since experimental data for the structures of the clusters at nonzero temperatures were not available, direct comparisons of the above simulations with experiments were not possible. Therefore, the accuracy and efficiency of ab initio MC simulations remain to be determined. Furthermore, the feasibility of the extension of ab initio MC simulations to bulk systems, which require models containing many more atoms than the small clusters, has yet to be tested. As a first application of ab initio MC simulations to bulk systems, we have chosen the liquid phase of bulk lithium. We have also investigated the structural properties of a cluster containing 16 lithium atoms at different temperatures. The simulated results for the lithium cluster are physically reasonable, and the structure of liquid lithium obtained from our ab initio MC simulations is in good agreement with the available experimental measurements. We also demonstrate that the computational demands of simulating liquid lithium with our ab initio MC algorithms are comparable to those of ab initio MD simulations.

The remainder of this paper is organized as follows. Section II describes our ab initio MC algorithm. Section III describes the simulation results for a 16-atom lithium cluster (Sec. III A) and for liquid bulk lithium (Sec. III B), with a discussion of general aspects and possible improvements of our ab initio MC algorithm in Sec. III C. Sec. IV contains our conclusion.

## II  AB INITIO MONTE CARLO ALGORITHM

To simulate the cluster of 16 lithium atoms, the atoms were put in a periodically repeated cubic supercell. The size of the supercell ($`10\times 10\times 10`$ Å3) was much larger than the concentrated cluster (see Fig. 1.). Bulk lithium was represented by a cubic box (supercell) with periodic boundary conditions containing 54 Li atoms. The size of the supercell ($`10.698\times 10.698\times 10.698`$ Å3 and $`10.874\times 10.874\times 10.874`$ Å3 for simulations at $`T=523`$ K and $`T=725`$ K, respectively) was determined to match the experimental density of liquid lithium at the given temperature. Ohse (1985) (ed.)

The simulations were performed in canonical ensembles. The standard Metropolis algorithm, Metropolis et al. (1953); Landau and Binder (2000) which provides an efficient approach for simulating the equilibrium properties of an atomic system at a given temperature, was used in our ab initio MC simulations. Starting from a pre-selected initial configuration, the simulations repeated the following steps. (1) Given the current configuration $`I`$ of the system, a new configuration $`J`$ is generated by making a small random displacement of one randomly chosen particle (suitable for small systems) or random displacements of all the particles simultaneously (more efficient for systems containing a large number of particles). (2) Once the total energies ($`E_I`$ and $`E_J`$) of these two configurations are calculated, the acceptance probability of the new configuration $`J`$ is then determined as

$$P(J|I)=\mathrm{min}[1,\mathrm{exp}\left(\frac{(E_JE_I)}{k_BT}\right)],$$
(1)
where $`T`$ is the temperature and $`k_B`$ is Boltzmann’s constant. (3) If the configuration $`J`$ is accepted, it serves as the current configuration of the next MC step, and $`E_I`$ is set equal to $`E_J`$. If the configuration $`J`$ is not accepted, the configuration $`I`$ and its energy $`E_I`$ are retained and used to start the next step. In this way, the system will eventually reach equilibrium and evolve toward a Boltzmann distribution. Metropolis et al. (1953); Landau and Binder (2000) The resulting sequence of sampling configurations is then used to obtain the equilibrium properties of the system.

To calculate the total energy at each MC step, we employed an ab initio pseudopotential density-functional-theory (DFT) method with a plane-wave basis set. All the ab initio total-energy calculations were carried out using the Vienna ab-initio simulation package (VASP). Kresse and Hafner (1993); Kresse and Furthmüller (1996a, b) The electronic exchange-correlation effects were treated with the generalized gradient-corrected functionals given by Perdew and Wang. Perdew et al. (1992); Perdew and Wang (1992) We adopted the Vanderbilt ultrasoft pseudopotential to replace the core electrons Vanderbilt (1990); Kresse and Hafner (1994) and used the conjugate-gradient technique for performing electronic relaxations Payne et al. (1992). The total-energy calculations were conducted with 4 k-points in the three-dimensional Brillouin zone with a plane-wave energy cutoff of 120 eV. The convergence of the total energy differences between different configurations was checked for selected configurations with higher cutoff energies (up to 200 eV) and more k-points (up to 16 k-points), and the convergence was found to be within a few meV. Because we used a pseudopotential approach, the total-energy calculations for the configurations with too small distances between atoms, which result in core-core overlap, would give wrong total energies. We therefore defined a hard-core distance of 2.0 Å between any two Li atoms. Trial configurations with any interatomic distance smaller than the cutoff distance were always rejected in the MC simulations. This is physically reasonable because such configurations would have very high potential energies due to the strong core-core repulsion.

## III RESULTS AND DISCUSSION

### III.1 Small cluster

Starting from the configuration with the atoms at ordered bcc positions (shown in Fig. 1), we first performed our ab initio MC simulations for the cluster of 16 lithium atoms at temperatures of 1000 K and 5000 K. The MC moves in the simulations for the cluster consisted of random displacement of one randomly selected lithium atom at each MC step, with the maximum displacement adjusted to give an overall acceptance rate of about 50%. The chosen atom was displaced from its old position with equal probability to any position inside a sphere surrounding the old position, with the radius of the sphere being the maximum displacement. We obtained maximum displacements of 0.75 Å and 1.2 Å for the simulations at the temperatures of 1000 K and 5000 K, respectively. About 10000 MC steps were needed to reach the equilibrium state. Averages for calculating the structural properties were taken over the next 10000 sampling configurations after equilibration. It took approximately 1 minute of CPU time to obtain the total energy of each configuration with a single 375 MHz processor on an IBM SP3 supercomputer.

The pair correlation function $`g(r)`$ can be defined as Allen and Tildesley (1987)

$$g(r)=\frac{2V}{N(N1)}\frac{1}{4\pi r^2}\underset{i}{}\underset{j<i}{}\delta (rr_{ij}),$$
(2)
where $`V`$ is the system volume, $`N`$ is the number of particles, $`r_{ij}`$ is the distance between a pair of particles $`i`$ and $`j`$, and $``$ signifies averaging over the sequence of simulated configurations. The pair correlation function is proportional to the probability of finding a pair of particles a distance $`r`$ apart, normalized by the square of the particle density such that it approaches unity in the limit of a uniform random distribution at the same density. It hence provides a measure of the local spatial ordering of the system. As described in Ref. 26, we determined $`g(r)`$ by dividing the range of $`r`$ (\[2 Å , 10 Å\]) into 400 intervals of length 0.02 Å ($`\mathrm{\Delta }r`$), calculating every $`r_{ij}`$ for a given configuration, and recording in narrow bins $`[r,r+\mathrm{\Delta }r)`$ the frequency of occurrence of the different particle separations. In this way a histogram was built up and finally normalized at the end of the scan over the configurations.

Figure 2 shows the pair correlation functions for the cluster at temperatures of 1000 K and 5000 K, respectively. It can be observed that there are two peaks in the pair correlation function at $`T=1000`$ K. This indicates the existence of two atomic shells for the structure, and therefore the system locally has a liquid-like structure. At a temperature of 5000 K, however, there is no obvious second peak in the pair correlation function, and thus the system has a gas-like structure. Given that the boiling point of bulk Li is about 1600 K, the results from our simulations for the cluster of lithium atoms are physically reasonable. Note that the pair correlation function $`g`$($`r`$) at $`T=5000`$ K does not approach zero when $`r`$ is close to 2.0 Å, indicating that the configurations with very small separations (smaller than 2.0 Å) between atoms would contribute to the ensemble averages. This reflects the gas-like structure of the system at $`T=5000`$ K because the probability of occurrence of configurations with very small separations between particles is expected to be higher for a gas-like structure than for a liquid-like structure.

### III.2 Bulk system

Following the simulations for the lithium cluster, we carried out ab initio MC simulations for the liquid phase of bulk lithium. In contrast to the simulations for the lithium cluster, the Monte Carlo moves in the simulations for liquid lithium consisted of random displacements of all the atoms in the supercell simultaneously. The other aspects of the MC moves were the same as in the simulations for the lithium cluster. In order to achieve an overall acceptance rate of roughly 50%, we found that the maximum displacements at temperatures of 523 K and 725 K were 0.05 Å and 0.06 Å, respectively. We started our simulations at $`T=725`$ K from an initial configuration with all the Li atoms at their bulk bcc positions. Roughly 8000 MC steps were required to reach equilibrium. Since the initial configuration was a perfect crystalline structure, far away from the equilibrium structure of the liquid phase of bulk lithium, it required relatively long MC runs to reach equilibrium. In contrast, the initial configuration for the simulations at $`T=`$ 523 K was chosen as one of the equilibrium sampling configurations obtained in the simulations at $`T=725`$ K. A reduced number of MC steps (5000 MC steps) was then required for the equilibration period. After equilibration at either temperature, we generated 5000 sampling configurations for averages.

Our simulations were carried out in parallel on four 375-MHz processors of an IBM SP3 supercomputer. The total-energy calculation at each MC step required approximately 150 seconds of CPU time. On average, the electronic relaxation converged within seven self-consistency steps. The software (VASP) that we used for the total-energy calculations contains computations for the forces after the electronic relaxation converges. Such computations are not needed in ab initio MC simulations. Should VASP be modified to exclude the force calculations, the computational costs would therefore be decreased by approximatively 10% – 15%.

For comparison, we also performed ab initio MD simulations Kresse and Hafner (1993); Kresse and Furthmüller (1996a, b) for liquid lithium. Canonical ensembles were generated with the algorithm of Nosé, Nosé (1984) and the time step was set to 1 fs. The initial configuration for the ab initio MD simulations at both $`T=`$ 725 K and $`T=`$ 523 K was the same perfect bcc structure as that for the ab initio MC simulations at $`T=`$ 725 K. The equilibration periods extended over 5000 and 3000 time steps at $`T=`$ 523 K and $`T=`$ 725 K, respectively. Ensemble averages were calculated from the sampling configurations of the next 5000 steps. Each step required about 400 seconds of CPU time on an IBM SP3 supercomputer (in parallel on 4 processors), and the corresponding wall-clock time was roughly 100 seconds.

The pair correlation functions for liquid lithium at temperatures of 523 K and 725 K, obtained from our ab initio MC and ab initio MD simulations are shown in Fig. 3 and Fig. 4, respectively. For comparison, the experimental data for liquid lithium at $`T=523`$ K obtained by X-ray diffraction Waseda (1980) are also shown in Fig. 3. We observe that the agreement between ab initio MC, ab initio MD, and experiment is very good. In particular, the three peaks determined from our ab initio MC simulations are essentially located at the experimental positions. The intensities of the peaks obtained from the ab initio MC simulations are also basically in agreement with those from both the ab initio MD simulations and the X-ray data. The small differences of the intensities of the peaks between the theoretical results and the experimental data may be due to the small size of our supercells. The experimentally observed temperature effects, Olbrich et al. (1983) that is, broadening of the peaks and slight shifts to higher positions of the second and third peaks at the higher temperature are also well reproduced by our ab initio MC simulations (see Fig. 3 and Fig. 4).

### III.3 General discussion

Finally, we discuss some general aspects and possible improvements of our ab initio MC algorithm. The computational costs of ab initio MC simulations are dominated by the ab initio total-energy calculations. Since efficient iterative schemes for ab initio total-energy calculations with a plane-wave basis set have achieved order $`N^2`$ ($`N`$ is the number of the atoms in the supercell) scaling for systems containing up to 1000 electrons, Kresse and Furthmüller (1996b) the cost of simulations with our ab initio MC algorithm would increase essentially as the square of the number of atoms.

While the computational costs are comparable, our ab initio MC algorithm is not any faster than ab initio MD simulations. In particular, a larger number of steps are required to reach equilibrium in our ab initio MC simulations for bulk lithium than in the corresponding ab initio MD simulations. For the lithium systems that we have investigated, our ab initio MC algorithm produces the required sampling configurations at a speed similar to the ab initio MD simulations. In this respect, our ab initio MC simulations are as efficient as ab initio MD simulations.

In some cases, like simulations for strongly covalent systems and situations involving volume changes, however, our ab initio MC algorithm may not be efficient due to the low acceptance rate of the traditional MC moves. In such situations, therefore, our ab initio MC algorithm needs to be improved. A cluster MC move scheme proposed by Lee and Swendsen Lee and Swendsen (2001) may be used for this purpose. We have tested the cluster MC move scheme for the bulk lithium systems that we investigated with traditional MC moves as described in Sec. III. We found that the use of cluster MC moves in ab initio MC simulations provided almost the same performance as the traditional MC moves, indicating that both kinds of MC moves might work equally well for metallic systems. It is also possible that the relative performance of different kinds of MC moves somewhat depend on the number of particles included in the supercell. It would be interesting to see if the cluster MC move scheme would be more efficient than the traditional MC approach in simulations for systems containing large numbers (for example, hundreds) of particles. In the cases where the boundary conditions or/and the volume of the simulated cells change, which result in a low acceptance rate with the traditional MC moves, the cluster MC move scheme has been shown to increase the acceptance rate significantly. Lee and Swendsen (2001) We also expect that the cluster MC move scheme would be more efficient than the traditional MC move method in simulations for systems containing strong covalent bonds. Another approach, a hybrid of MD and MC simulation in which trial MC moves are generated with classical MD methods, may also enhance the efficiency of ab initio MC simulations. Duane et al. (1987) Other techniques aiming at sampling important configurations in the configuration space include the J-walking procedure proposed by Frantz et al. Jellinek et al. (1998); Srinivas and Jellinek (2000); Frantz et al. (1990) and a big-move method suggested by Akhmatskaya et al. Akhmatskaya et al. (1997).

## IV  CONCLUSION

In conclusion, we have shown that the use of ab initio calculations in classical Monte Carlo simulations provides structures of liquid lithium in agreement with experimental measurements and ab initio MD simulations. Our work also demonstrates that ab initio MC simulations for condensed-matter systems with up to several tens of atoms (and possibly with more than a hundred atoms) are feasible and may provide an alternative way to investigate finite-temperature properties of such systems.

###### Acknowledgements.

We thank G. Brown, Y.Z. Cao, Z. Liu, M.A. Novotny, and H.K. Lee for helpful discussions. This work was supported by the National Science Foundation under grants DMR-9981815, DMR-0111841, and DMR-0240078, and by Florida State University through the Center for Materials Research and Technology and the School of Computational Science and Technology. The calculations were performed on IBM SP3 computers at Florida State University and the National Energy Research Scientific Computing Center, which is supported by the U.S. Department of Energy.


format: text

text:
# Φ meson production at RHIC

## Abstract

The production of $`\varphi `$ mesons in Au+Au collisions at RHIC and their propagation in a hot and dense nuclear medium is studied within the microscopic quark-gluon string model. The inverse slope parameter of the transverse mass distribution agrees well with that extracted from the STAR data, while the absolute yield of $`\varphi `$ is underestimated by a factor 2. It appears that the fusion of strings alone cannot increase the $`\varphi `$ yield either. Less than 30% of detectable $`\varphi `$’s experience elastic scattering, this rate is insufficient for the full thermalization of $`\varphi `$. The directed flow of $`\varphi `$ at $`|y|2`$ demonstrates strong antiflow behavior, whereas its elliptic flow rises up to about 3.5% in the same rapidity interval. As a function of transverse momentum it rises linearly with increasing $`p_t`$, in agreement with the STAR data, and saturates at $`p_t>2`$ GeV/c.

The investigation of nuclear matter under extreme temperatures and densities, and the search for a predicted transition to a deconfined phase of quarks and gluons, the so-called Quark-Gluon Plasma (QGP), is one of the main goals of heavy ion experiments at ultrarelativistic energies. Both theorists and experimentalists are looking for genuine QGP fingerprints, that cannot be masked or washed out by processes on a hadronic level. The $`\varphi `$ meson was proposed about two decades ago  as one of the most promising QGP messengers because of the following reasons: Firstly, an enhancement of $`s\overline{s}`$ pairs in the QGP phase should lead to an enhancement of $`\varphi `$ mesons, as well as other strange hadrons, while the production of $`\varphi `$ in $`hh`$ interactions is suppressed due to the OZI rule . Secondly, since the $`\varphi `$ interaction cross section with non-strange hadrons is small, $`\varphi `$ will keep information about the early hot and dense phase. Thirdly, the $`\varphi `$ meson spectrum is not distorted by feeddown from resonance decays. Finally, among other channels $`\varphi `$ decays into a pair of kaons and, more rarely, into a lepton pair. Both channels have been explored for the reconstruction of the $`\varphi `$ yield in Pb+Pb collisions at SPS energy. It was found that the $`\varphi `$ multiplicity extracted from the di-kaon data  is significantly smaller than that obtained from the dimuons  (the so-called $`\varphi `$ meson puzzle). A similar effect has been observed recently by the PHENIX Collaboration  in Au+Au collisions at RHIC. A possible solution of this puzzle is the scattering of at least one of the daughter kaons in the nuclear medium  accompanied by in-medium modifications of kaon and $`\varphi `$ masses. Note also that the width and mass of the $`\varphi `$ meson are quite sensitive to the strangeness content of the surrounding medium.

In the present paper we study production, propagation, and freeze-out of $`\varphi `$ mesons in Au+Au collisions at $`\sqrt{s}=130`$ AGeV within the microscopic quark-gluon string model (QGSM) . The transverse mass spectra of $`\varphi `$ calculated for three different centrality bins are presented in Fig. 1 together with the experimental data reported by the STAR

data. The agreement between the calculations and the data can be improved by increasing the string tension or dropping the constituent quark masses , by including the OZI-violation mechanism in the model (at low $`\sqrt{s}`$ this is the only channel for $`\varphi `$ production), and by increasing the number of $`s\overline{s}`$-pairs inside nucleons in multipomeron diagrams .

The $`\varphi `$ mesons are produced either in primary nucleon-nucleon interactions via string fragmentation processes or in subsequent rescattering between produced strange hadrons. We found that in the first 4 fm/$`c`$ the $`\varphi `$ mesons are produced solely in string processes. The string channel dominates over the rescattering one up to 20 fm/$`c`$; in the later stages the rescattering channel prevails. This dynamical picture agrees well with that obtained for the $`\varphi `$ production at SPS energy within the UrQMD model . The $`dN/dt`$ distributions of $`\varphi ,N,\pi ,K,\mathrm{\Lambda }`$, which are decoupled from the fireball after the last elastic or inelastic collision, are shown in Fig. 3. Compared to AGS and SPS energies , a substantial part of hadrons leave the system immediately after their production. Mesonic distributions have sharp peaks at $`t=810`$ fm/$`c`$, while the distributions of baryons are wider due to the large number of rescatterings, which shift their $`dN/dt`$ maxima to later times. Figure 3 depicts the freeze-out distribution of the $`\varphi `$ mesons over rapidity and emission time. We see that the bulk production of $`\varphi `$ occurs in the central rapidity window $`|y|1.5`$ within the first 10 fm/$`c`$, while in $`|y|2`$ areas the $`\varphi `$ mesons are produced until the

late stages of the system evolution: Because of the rapid expansion the central zone becomes quickly dilute, its energy density drops, and early freeze-out of $`\varphi `$ mesons, which have small interaction cross section with non-strange matter, takes place. According to QGSM, about 51% of detectable $`\varphi `$ mesons decay after decoupling from the fireball. Less than 30% of $`\varphi `$ mesons experience at least one elastic rescattering before the freeze-out. This relatively low rate of elastic collisions may not be sufficient for the thermalization of $`\varphi `$ mesons, however, would it be enough for the development of noticeable anisotropic flow? - Other processes contributed to the flow evolution are production and absorption/decay of $`\varphi `$ in dense medium. While the elastic collisions and resonance production of $`\varphi `$ mesons increase the flow in the normal bounce-off direction, the $`\varphi `$ absorption channel effectively reduces this normal component thus elongating the flow in the antiflow direction.

Figure 5 shows the rapidity dependence of the directed flow of $`\varphi ,N`$, and $`K`$ in minimum bias Au+Au collisions at RHIC. The slopes of the all distributions are negative at $`|y|2`$, i.e. the antiflow component of the $`v_1`$ dominates over its normal counterpart (see ). Similar antiflow slopes of $`v_1(y)`$ are developed by $`\pi `$ and $`\mathrm{\Lambda }`$ ; its origin is traced to nuclear shadowing. At midrapidity $`|y|0.5`$ the directed flow of $`\varphi `$ is quite weak.

The elliptic flow of $`\varphi `$ mesons as a function of $`\eta `$ and $`p_t`$ is shown in Fig. 5 together with the model predictions for $`\pi `$ and $`N`$ , and recent experimental data . The magnitude of the $`\varphi `$ flow at midrapidity $`v_2^\varphi (\eta )=3.3\pm 1.1\%`$ is similar to that of pions and nucleons. At larger rapidity $`v_2^\varphi (y)`$ is weaker than the flow of other hadrons. As a function of $`p_t`$ the elliptic flow of $`\varphi `$ rises up to 2 GeV/$`c`$ almost linearly with $`p_t`$. The agreement with the experimental results is good, although large error bars do not allow more definite conclusions. The saturation and/or possible drop of the $`v_2^\varphi (p_t)`$ at higher $`p_t`$ is linked to the lack of rescattering in spatially anisotropic finite nuclear matter.

In summary, the QGSM successfully describes $`m_t`$-spectra of $`\varphi `$ in different centrality bins in Au+Au collisions at $`\sqrt{s}=130`$ AGeV, but underestimates the total $`\varphi `$ production. This problem cannot be cured alone by the fusion of strings. Mesons are frozen in the model earlier than baryons. The QGSM predicts unique antiflow slope at midrapidity for the directed flow of $`\varphi ,\pi ,N,K,\mathrm{\Lambda }`$, and reproduces correctly the elliptic flow $`v_2^\varphi (p_t)`$.

Acknowledgments. Fruitful discussions with A. Capella, A. Kaidalov, M. Krivoruchenko, D. Sousa and S. Voloshin are gratefully acknowledged. This work was supported by the EGK Basel-Tübingen under DFG contract GRK683, by the BMBF under contract 06TÜ986, and by the Bergen Computational Physics Laboratory (BCPL) in the framework of the European Community - Access to Research Infrastructure action of the Improving Human Potential Programme.


text:
# Weak Gravitational Lensing by a Sample of X-ray Luminous Clusters of Galaxies – III. Serendipitous weak lensing detections of dark and luminous mass concentrations.

## 1 Introduction

Weak gravitational lensing provides a powerful way to identify cluster-sized density peaks in the Universe, independent of their baryonic content. Given the currently modest sky coverage of optical imaging surveys with the depth and image quality required to detect new clusters by their weak lensing effect, the number of currently known mass-selected clusters is very small. It is still an open question whether the mass-selection will lead to the identification of a population of clusters which are physically different from optically selected clusters or X-ray detected clusters. If a population of “baryon-poor” clusters is found to exist, they may be very useful laboratories for the study of dark matter properties. For instance, some dark matter candidates such as sterile neutrinos may produce an observable signature from their decay (Abazajian, Fuller, & Tucker 2001). The best places to detect such a signature would be in baryon-poor clusters – if such objects exist – where the spectral line corresponding to dark matter decay would be relatively more conspicuous compared to the emission produced by bremsstrahlung in the hot intra-cluster gas (Hansen et al. 2002). In any case, the existence of baryon-poor clusters or sub-clusters would pose a challenge to current models for structure formation. Furthermore, any previously unrecognized population of clusters with high mass-to-light ratios would have to be taken into account when using the measured average mass-to-light ratios of clusters to estimate the density parameter $`\mathrm{\Omega }_m`$.

The present sample of weak lensing-detected clusters is small and contains both clusters with “normal” mass-to-light ratios and objects which appear to be optically dark. From weak lensing observations in the field of the cluster A1942 (catalog ), Erben et al. (2000) find a secondary mass peak $`7\mathrm{}`$ south of the cluster center which does not correspond to any strong concentration of bright galaxies. From deep $`H`$-band imaging of this region, Gray et al. (2001) constrain the bolometric mass-to-light ratio to be $`M/L_{\mathrm{Bol}}>1000h`$ in solar units for any reasonable lens redshift. Umetsu & Futamase (2000) find a dark mass concentration $`1\stackrel{}{\mathrm{.}}7`$ south of the high-redshift cluster CL1604+4304 (catalog ) ($`z=0.90`$) with an estimated mass of $`1.2\times 10^{14}h^1M_{}`$ and $`M/L_B1000h`$ in solar units, if it is located at the redshift of CL1604+4304 (catalog ). At present, it is not clear whether the dark clumps found by Erben et al. and Umetsu & Futamase are physically associated with the nearby clusters, or whether they represent chance alignments on the sky of objects at different redshifts. Most recently, Miralles et al. (2002) have reported evidence for another dark cluster from a conspicuous alignment of faint galaxies in a parallel STIS pointing adjacent to the local Seyfert galaxy NGC 625 (catalog ).

Wittman et al. (2001) report the discovery of a more “normal” cluster in a “blank sky” field through a combination of weak gravitational lensing and photometric data. Their $`BVRI`$ photometry shows a concentration of elliptical galaxies close to the lensing-derived mass peak corresponding to the cluster, and spectroscopic follow-up of candidate cluster members reveals a cluster with modest galaxy velocity dispersion ($`\sigma _v=615\pm 150\mathrm{km}\mathrm{s}^1`$) at $`z=0.28`$. The mass-to-light ratio of this cluster is $`M/L_R=560\pm 200h`$, which is somewhat high compared to average values of $`M/L_B300h`$ obtained from both virial and weak lensing analyses of X-ray selected clusters (e.g., Carlberg et al. 1997; Hoekstra et al. 2002; Dahle et al. 2003, in prep.), but there are some X-ray selected clusters with similar lensing-derived $`M/L`$, such as MS 1224.7+2007 (catalog ) at $`M/L_R=640\pm 150`$ (Fahlman et al. 1994; Fischer 1999), and A68 (catalog ) and A697 (catalog ) at $`M/L_R=680\pm 230h`$ and $`M/L_R=450\pm 115h`$, respectively (Dahle et al. 2003, in prep).

In Paper I in this series (Dahle et al. 2002) we presented weak lensing measurements of a sample of 39 X-ray selected clusters. The results were presented in the form of maps of the reconstructed projected matter density ($`\kappa =\mathrm{\Sigma }/\mathrm{\Sigma }_{\mathrm{crit}}`$, where $`\mathrm{\Sigma }_{\mathrm{crit}}`$ is the critical surface mass density for lensing) and in the form of radial mass profiles around each (lensing-determined) cluster center. Several mass maps show evidence for sub-peaks in the mass distribution which are not associated with obvious sub-clumping of optically luminous galaxies inside the cluster. In this paper, we investigate the properties of the most significant of these sub-peaks. We attempt to constrain the redshift and mass-to-light ratio of these systems and discuss whether they are likely to be physically associated with their apparent “host clusters”.

In §2 we describe the selection criteria for cluster candidates, which we describe individually in §3. In §4 we compare our new (sub-)cluster candidates to those found by other groups, and compare the observed abundances and physical properties of such objects to recent theoretical predictions.

The numbers in this paper are given for an Einstein-de Sitter ($`\mathrm{\Omega }_m=1`$, $`\mathrm{\Omega }_\mathrm{\Lambda }=0`$) cosmology. The Hubble parameter is given by $`H_0=100h\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$, and all celestial coordinates are given in J2000.0.

## 2 Identification of new cluster candidates

In Paper I we described a data set of weak lensing measurements of 39 clusters in the redshift range $`0.15<z<0.35`$, with selection based on very high X-ray luminosity ($`L_{\mathrm{x},0.12.4\mathrm{keV}}10^{45}`$ erg s-1). Most of the sky coverage of this survey comes from imaging of fourteen clusters with the UH8K mosaic CCD camera at the 2.24m University of Hawaii Telescope at Mauna Kea Observatory, which gives a $`19\mathrm{}`$ field when mounted at the f/10 Cassegrain focus (see Paper I for a complete list of the observations). The imaging data were reduced according to the “pipeline” reduction procedure for mosaic CCD data described by Kaiser et al. (1999), and the weak gravitational lensing was measured using the shear estimator introduced by Kaiser (2000). The reduction and analysis procedures used are described in detail in Paper I.

In some of the reconstructed mass maps of the fields observed with the UH8K camera, there are secondary mass peaks which appear significant when compared to peaks in randomized shear maps (see Figure 48 of paper I), or peaks seen in mass reconstructions of blank fields (see Wilson, Kaiser & Luppino 2001), but do not correspond to concentrations of early-type galaxies at the cluster redshift. These may be mass concentrations corresponding to high-redshift clusters that only contribute weakly to the light distribution in the field, compared to the more nearby target clusters. At redshifts $`>0.4`$, the early-type galaxies at a given redshift will be significantly bluer than predicted from the approximate linear $`VI`$ color-redshift relation used in Paper I. Thus, their contribution to the prediction for the dimensionless surface density given for each cluster in Paper I will be increasingly underestimated with increasing redshift.

An even more interesting possibility is that these are optically and/or X-ray dark clusters at intermediate redshifts (or dark sub-clumps physically associated with their apparent host clusters), i.e., mass concentrations that have much higher $`M/L`$ values than normal clusters.

### 2.1 Lensing Cluster Search

As noted in Paper I, the $`1\sigma `$ level of noise in the random mass maps is in the range $`0.035<\kappa <0.05`$, mostly depending on seeing. Given the large number of fields covered by our survey, we may expect to see maybe one or two $`3\sigma `$ peaks that are generated by random noise due to the intrinsic galaxy shapes. However, the reality of a mass peak at this level will become significantly more certain if it can be shown to be strongly correlated with clustering in the visible galaxy distribution, or with a peak in the X-ray luminous gas.

From the reconstructions of the dimensionless projected mass surface density $`\kappa `$ presented in Paper I, we select candidates for additional mass concentrations using the following criteria:

* The dimensionless projected surface density $`\kappa `$ should be $`>0.15`$ (corresponding to a $`>3\sigma `$ detection of mass even in the fields with the highest noise) at the position of the peak. This selection criterion was applied to mass maps derived separately from $`V`$\- and $`I`$-band imaging data, and to mass maps derived from a combination of the two, as detailed in Paper I. We found a total of three mass peaks which are above the $`\kappa `$ threshold in at least one of the different mass maps and which also satisfy the other three selection criteria below.
* The mass peak should be clearly separated (with offset $`>1\mathrm{}`$) from the mass peak that corresponds to the light peak/optical cluster center (or any other significant peak in the light distribution of galaxies belonging to a previously known cluster, as determined from their $`VI`$ colors).
* Both $`V`$\- and $`I`$-band photometry should be available at the position of the clump. Given the data set of Paper I, this implies that the total sky coverage of our search is 1.0 deg2. The dominant contribution to the sky coverage comes from $`19\mathrm{}`$ wide UH8K fields, and all the mass peaks discussed here are located in such fields.
* Candidates within $`1\mathrm{}`$ from the field edges are discarded, as the noise in weak lensing reconstructions are known to increase close to the edges.

The limit on the value of the mass surface density $`\kappa `$ is partly motivated by the study of Wilson, Kaiser & Luppino (2001), who present weak lensing data of similar quality to the data set used in this paper for 1.5 deg2 of “blank sky”. See Kaiser, Wilson & Luppino (2000) and Paper I for more details about the characteristics of each data set. After smoothing their reconstructed mass maps with a $`45\mathrm{}`$ Gaussian filter (we use a $`47\mathrm{}`$ Gaussian filter) they do not find any pixel values above $`\kappa =0.15`$. If the sky density of optically dark high-density mass peaks was found to be significantly higher in the fields of known clusters (e.g., the fields studied in Paper I) than in blank fields (e.g, the fields of Wilson et al. 2001), it would be a clear indication of a physical association between the dark clumps and “normal” clusters.

We note that the mass maps presented in this paper differ slightly from the mass maps presented for the same clusters in Paper I. The difference is in the relative weighting of shape measurements of faint galaxies in the two passbands when estimating the shear. In Paper I, the mass reconstruction was based on shear estimates from a combined “$`V+I`$” catalog which was made by selecting the catalog entry from the passband in which the object was detected with the highest significance. The maps presented in this paper are generated by taking a weighted average (where the weights are given by the “figure of merit” values $`Q^2/d\mathrm{\Omega }`$; see Paper I and Kaiser 2000) of this “$`V+I`$” mass map and the mass maps generated separately from the $`I`$-band and $`V`$-band catalogs. The $`1\sigma `$ level of noise in these maps at the positions of the newly discovered mass peaks is $`\kappa =0.03`$.

### 2.2 Photometric cluster search

For the UH8K fields which have both $`I`$-band and $`V`$-band data covering the whole field, the nature of the lensing-derived mass peaks can be investigated by looking for associated concentrations of red galaxies, corresponding to galaxy clusters up to $`z1`$.

Our photometric procedures are described in Paper I and are only briefly recounted here. We used SExtractor (Bertin & Arnouts 1996) to derive total $`V`$\- and $`I`$-band magnitudes for galaxies in the fields, and we also measured a $`VI`$ color within a fixed aperture of $`2\stackrel{}{\mathrm{.}}7`$. The reason for calculating the color within this smaller aperture rather than based on the total magnitudes is that the aperture magnitudes are less sensitive to contamination by other objects and to color gradients within the galaxies themselves and thus yield a tighter red sequence in a color-magnitude diagram. The magnitudes were corrected for Galactic extinction using the dust maps of Schlegel, Finkbeiner & Davis (1998), assuming an $`R_V=3.1`$ extinction curve.

To identify optical cluster candidates based on our two-filter photometry, we used a simplified version of the cluster-red-sequence (CRS) method used by Gladders & Yee (2000) for their Red-Sequence Cluster Survey. This cluster finding algorithm takes advantage of the fact that early-type galaxies in dense cluster environments form a remarkably tight, almost horizontal sequence in their color-magnitude diagrams, and these sequences are strongly homogeneous at a given redshift (e.g., Smail et al. 1998; Gladders & Yee 2000 and references therein). Also, at redshifts $`z1`$ the early-type galaxies that form the red cluster sequence will have a redder $`VI`$ color than any other normal galaxies at redshifts equal to, or lower than, the redshift of the cluster.

The basic idea behind the CRS method is to create a series of smoothed maps of the galaxy density in the field by selecting galaxies in a series of partially overlapping slices in the color-magnitude diagram. Gladders & Yee (2000) used color slices defined by theoretical predictions of the red sequence, but we choose here to neglect the slope of the red sequence and use a series of horizontal slices evenly spaced in $`VI`$. Both theoretical predictions from population synthesis models and empirical results based on our data set show the absolute value of the red-sequence slope to be less than 0.1 per $`I`$-magnitude. Given the relatively narrow magnitude range (particularly at high redshifts) of cluster galaxies contributing to a given slice, and the fact that our color slices are approximately twice as wide as those used by Gladders & Yee (2000), the effect of the slope will not significantly influence our ability to detect cluster candidates with this algorithm. However, it would probably bias redshift estimates of cluster candidates, but only at levels comparable to the systematic uncertainty introduced by the photometric calibration of the UH8K mosaic camera (we note that the photometric stability of the CCD chips in the UH8K is questionable at the $`0.1`$ magnitude level).

For each slice, we include all galaxies at $`I<23`$ that have a probability higher than 10% of really belonging to that slice when taking the photometric uncertainties into account. The width of the slices was $`\mathrm{\Delta }(VI)=0.3`$ at bright magnitudes and slightly larger at faint magnitudes due to the increasing photometric uncertainties. This is about twice the measured intrinsic width of the red sequence (Stanford, Eisenhardt, & Dickinson 1998).

Rather than using a more elaborate galaxy weighting scheme such as the one outlined by Gladders & Yee (2000), we chose a simple faint galaxy cutoff at $`I=23`$. Also, to generate smoothed galaxy number density plots comparable to the weak lensing-derived mass maps, we employed a Gaussian smoothing kernel of fixed scale $`47\mathrm{}`$ rather than a kernel tailored to fit the expected cluster density profile at a given redshift. The primary purpose of generating these plots is to search for galaxy density enhancements at the position of the lensing-detected cluster candidates.

### 2.3 Lensing masses and mass-to-light ratios

For each (sub-)cluster candidate, the (minimum) mass and the mass-to-light ratio were estimated within an aperture centered on the corresponding peak found in the $`\kappa `$ map (see Table 1). A lower limit on the aperture mass was found using the aperture densitometry statistic $`\zeta `$ of Kaiser et al. (1994):

$$\zeta (R_1,R_2)=\overline{\kappa }(R_1)\overline{\kappa }(R_1,R_2)=\frac{2}{1R_1^2/R_2^2}_{R_1}^{R_2}\gamma _Td\mathrm{ln}r,$$
(1)
which measures the average surface mass density within an aperture radius $`R_{\mathrm{ap}}=R_1`$ minus the average surface mass density in a surrounding annulus $`R_1<R<R_2`$ (here, $`\gamma _T`$ is the tangential component of the shear). Hence, a lower limit on the aperture mass can be simply calculated as $`M_{\mathrm{ap}}=\pi r^2\zeta \mathrm{\Sigma }_{\mathrm{crit}}`$. We chose $`R_2=2R_1`$, and our estimate of $`M_{\mathrm{ap}}`$ will in this case e.g., underestimate the true aperture mass by 33% for an isothermal sphere ($`\mathrm{\Sigma }R^1`$) mass distribution.

In order to minimize the bias caused by the presence of a more massive cluster in each field we chose to use an outer radius $`R_2`$ much smaller than the $`R_2`$ values used for the aperture mass measurements presented in Paper I. For this reason, $`R_2`$ was chosen to be less than half the angular separation between the new cluster candidate and the Abell cluster. For this geometry, we find from simple simulations that the aperture mass measurement is affected by $`<1`$% , provided that the mass of the Abell cluster is not more than twice that of the new cluster candidate. In the case of the A959 (catalog ) and A1705 (catalog ) fields, we used a control annulus with $`R_2=2\mathrm{}`$ and inner radius $`R_1=1\mathrm{}`$. The cluster candidate in the A1722 (catalog ) field had a somewhat larger separation from its host cluster, and it was therefore possible to use a larger aperture with $`R_2=4\mathrm{}`$ and $`R_1=2\mathrm{}`$. By defining an aperture measurement of the surface light density $`\zeta _L`$ in a manner similar to equation (1) (i.e., generating an aperture luminosity estimate by subtracting the average galaxy luminosity density in the annulus $`R_1<R<R_2`$ from the average galaxy luminosity density at $`R<R_1`$), the mass-to-light ratio of the mass peak can be estimated. Thus, in contrast to $`M_{\mathrm{ap}}`$, this ratio will not be systematically underestimated.

### 2.4 Hot gas content

In X-ray luminous clusters of galaxies, such as those in the cluster sample of Paper I, most of the baryonic mass is in the form of hot, X-ray emitting gas (e.g., Ettori & Fabian 1999). In order to constrain the hot gas content of the (sub-)cluster candidates we analyzed the best quality archive X-ray data (standard screened EINSTEIN IPC, ROSAT PSPC and ASCA SIS data, respectively) for these candidates. No obvious X-ray emission is associated with any of the (sub-)cluster candidates so we determine upper limits for their X-ray luminosity as described below.

An upper limit for the (sub-)cluster X-ray count rate is extracted within two apertures: The aperture as determined from the weak lensing maps ($`R_{\mathrm{ap}}`$) and the radius within which the average (sub-)cluster mass density is 500 times the critical density of the Universe ($`r_{500}`$). Special care is taken in estimating the background.

For EINSTEIN IPC data (A1705) and ROSAT PSPC data (A959) the background is estimated from the data set including the new mass concentration itself. The background is calculated as the average of three source-free regions of the same size as the (sub-)cluster region, and at the same optical axis distance.

For ASCA SIS data (A1722) the background count rate is evaluated from the publically available “blank sky” data in a two stage process: A first background estimate is obtained using exactly the same region on the SIS chip as the (sub-)cluster region. Then, the general background level in the “blank sky” data is compared to the general background level in the (sub-)cluster data by calculating the average count rate within three regions of the same size as the region of interest, at the same optical axis distance and devoid of obvious X-ray sources. The inital background estimate is then corrected for the average difference in background level between the “blank sky” data and the (sub-)cluster data. Finally, data from the SIS0 detector and the SIS1 detector are merged into one common SIS data set.

For each (sub-)cluster, we derive net count rates (corrected for point source contributions) within the two apertures and convert these to 0.1 – 2.4 keV fluxes using PIMMS and an absorbed Raymond-Smith spectral model with the following parameters: Abundance is $`Z=0.25Z_{\mathrm{}}`$, absorption by neutral hydrogen is fixed at the Galactic value (from Dickey and Lockman 1990), and for each (sub-)cluster the temperature is fixed at the value derived from the weak lensing mass within $`r_{500}`$ using the mass-temperature relation of Finoguenov et al. (2001).

Assuming the redshift of the (sub-)cluster given in Table 1 based on optical data, the flux is converted to a luminosity for the desired cosmology. In order to determine the X-ray luminosity of the (sub-)clusters relative to a sample of nearby “standard” clusters, the upper limit on the 0.1 – 2.4 keV luminosity within $`r_{500}`$ was compared to the expected value based on the lensing mass within $`r_{500}`$ and the mass-luminosity relation of Reiprich & Böhringer (2002) (best bisector fit).

## 3 Candidate Mass Concentrations

The selection criteria outlined in § 2.1 produced three new cluster candidates, in the fields of A959 (catalog ), A1705 (catalog ) and A1722 (catalog ). To make the following discussion clearer and more concise, we introduce the naming convention “WL HHMM.M(+/-)DDMM” for these objects, where “WL” denotes a weak lensing-detected mass concentration and is followed by numbers that are based on the celestial coordinates of the objects (for epoch J2000.0), where the right ascension is given in units of hours and minutes (with one decimal) and the declination is given in units of degrees and arcminutes. For each candidate, we analyze its optical and lensing properties based on the data presented in Paper I along with X-ray data to search for emission from any associated hot gas. Table 1 summarizes the most important parameters for each cluster candidate.

### 3.1 The A959 Field

From spectroscopic observations at the WHT, we have determined a redshift $`z=0.286`$ for this cluster (Irgens et al. 2002), significantly lower than the value usually quoted in the literature ($`z=0.353`$). After properly revising its X-ray luminosity and temperature, it is found to have $`L_{\mathrm{x},0.12.4\mathrm{keV}}=1.43\times 10^{45}`$ erg s-1 (Böhringer et al. 2000) and $`T=6.6`$ keV (Mushotzsky & Scharf 1997). A singular isothermal sphere (SIS)-model fit to the tangential shear measured out to a radius of $`1.5h^1`$Mpc shows that the cluster has a mass corresponding to a velocity dispersion $`\sigma =990_{110}^{+100}`$km s-1. The fit to a SIS-model is however not very good, and weak lensing mass measurements at large radii indicate that this cluster is very massive, with a projected 2D aperture mass $`M_{\mathrm{ap}}(<1.5h^1\mathrm{Mpc})2\times 10^{15}M_{}`$ (see Paper I).

Optically, the cluster center is not dominated by any single galaxy, but it has a core region consisting of many early-type galaxies of similar brightnesses. A highly significant mass peak (with some evidence for substructure) is seen in the weak lensing mass map (see Figure 1), and the dark matter distribution appears to resemble the cluster light distribution, as apparent in the first four panels of Figure 2. However, there is a dark matter concentration, hereafter denoted as WL 1017.3+5931 (catalog ), $`6\mathrm{}`$ southwest of the cluster center which does not correspond to any peak in the luminosity distribution. Its peak value111As noted in § 2.1, the mass maps used for the selection of mass peaks were different from those displayed in this paper. The reason for not using the new maps presented here (and lowering the selection threshold from $`\kappa =0.15`$ to $`0.14`$, which would have selected exactly the same three objects!) is that the noise is more clearly defined in the old maps. While the new maps have lower noise, this is not achieved by increasing the available sample of background galaxies, but instead by in effect taking a weighted average of the shape measurements of these galaxies in the two passbands. Hence, the photon counting noise for the shape measurements of faint galaxies is somewhat reduced, while the noise due to the finite number density of background galaxies does not change. is $`\kappa =0.143`$.

Due to the very low background, the ROSAT PSPC 15.4 ksec exposure yields the tightest current constraints on diffuse X-ray emission from WL 1017.3+5931 (catalog ). In Figure 1, a couple of X-ray sources are seen close to WL 1017.3+5931 (catalog ). Archival ROSAT HRI images show that sources “C” and “D” are extended sources while source “A” is a point source and source “B” consists of two point sources on top of some extended emission. This extended X-ray emission is possibly linking WL 1017.3+5931 (catalog ) to the main cluster of A959 (catalog ). The mass morphology of WL 1017.3+5931 (catalog ) provides clearer evidence for this link, since there appears to be a bridge of matter extending from the A959 (catalog ) cluster center toward WL 1017.3+5931 (catalog ), and the peak itself appears highly elongated in the direction toward the center of A959 (catalog ). Even at the lowest “saddle point” in the apparent mass bridge between A959 (catalog ) and WL 1017.3+5931 (catalog ), the projected mass density is still as high as $`\kappa =0.122`$, which is $`85\%`$ of the peak value at the location of WL 1017.3+5931 (catalog ). Simulations of the region between two cluster mass peaks with a Gaussian smoothing scale equal to that of the real mass map and with random realizations of realistic noise indicate that this feature is (marginally) significant, at almost the $`2\sigma `$ level. Although the reality of this feature is currently debatable, the apparent mass bridge and the ellipticity of WL 1017.3+5931 (catalog ) probably constitute the only significant evidence so far for a physical link between a dark sub-clump found by weak lensing in a cluster field and its apparent host cluster.

No strong galaxy overdensities are seen at the position of WL 1017.3+5931 (catalog ) for any of the color slices plotted in Figure 2. The lower limit value for $`M/L_B`$ in Table 1 is based on an aperture light estimate using galaxies in the $`2.0<VI<4.0`$ color range. These are all redder than the A959 (catalog ) cluster galaxies (the red cluster sequence has $`VI1.7`$). If we instead consider galaxies in the color range of the red cluster sequence of A959 (catalog ), we find that the luminosity density in the control annulus at $`1\mathrm{}<R<2\mathrm{}`$ is higher than the luminosity density within $`1\mathrm{}`$, implying that there is no overdensity of light associated with A959 (catalog ) cluster galaxies at the location of WL 1017.3+5931 (catalog ). Two of the brightest galaxies in the field ($`I16.5`$) are visible $`2\mathrm{}`$ from the center of WL 1017.3+5931 (catalog ), but these are likely to be foreground ellipticals ($`VI1.3`$), and no significant galaxy clustering is seen around these two galaxies. The estimated velocity dispersion of WL 1017.3+5931 (catalog ) listed in Table 1 comes from a SIS-model fit to the tangential weak shear averaged in radial bins centered on the dark clump. If WL 1017.3+5931 (catalog ) is at the A959 (catalog ) main cluster redshift, within a radius of $`r_{500}`$ it is a factor three X-ray under-luminous relative to the expectation from the mass-X-ray luminosity relation of Reiprich and Böhringer (2002). However, given the scatter in the mass-luminosity relation of an order of magnitude, even though WL 1017.3+5931 (catalog ) is relatively X-ray dark it is still consistent with the relation. Chandra or XMM-Newton observations are required to pin down the X-ray emission of WL 1017.3+5931 (catalog ) and accurately constrain the hot gas content, ie. if it is really X-ray dark.

In general, there is very good morphological agreement between the mass map and hot gas as mapped by ROSAT PSPC: The mass extension towards the northeast has a clear X-ray counterpart, and at the position of the extended X-ray source “D” a mass concentration is seen. This feature also corresponds to an overdensity of galaxies (most prominent in the $`1.6<VI<1.9`$ color slice in Figure 2) with colors similar to, or slightly redder than, the early-type galaxies in the core of A959 (catalog ). The general southeast-northwest elongation of the A959 (catalog ) cluster galaxy distribution is matched in X-rays.

We also note that there are several concentrations of red ($`2.0<VI<3.0`$) galaxies in the field, most clearly visible in the lower three panels of Figure 2. One of these is only $`1\mathrm{}`$ north of the A959 (catalog ) cluster center (see the $`2.65<VI<2.95`$ color slice in Figure 2) and may contribute to the substructure seen in the mass peak of A959 (catalog ).

### 3.2 The A1705 Field

Optical images show that the cluster center of A1705 (catalog ) is dominated by a single cD galaxy. The $`\kappa `$ reconstruction of this field (see Figure 3.1) shows a number of secondary peaks in addition to the most significant mass peak, which is associated with the center of A1705 (catalog ). The position given for A1705 (catalog ) by Abell (1958) is $`5\mathrm{}`$ off from what we take to be the cluster center, based on the peaks in the mass, light and galaxy number density distributions, while the Zwicky cluster Zw 1312.1+7311 (catalog ) is more consistent with the position we determine for the A1705 (catalog ) cluster center (at $`\alpha =13^\mathrm{h}13^\mathrm{m}.0,\delta =72\mathrm{°}55\mathrm{}`$). The color slicing in Figure 4 shows several interesting effects: Firstly, the morphology-density relation manifests itself in the effect that the density peak associated with A1705 (catalog ) appears more concentrated as one moves to progressively redder color intervals (which contain the most evolved early-type cluster galaxies). This effect is also apparent in the similar plots for A959 (catalog ) and A1722 (catalog ). Secondly, the density peak of A1705 (catalog ) in Figure 4 shows a significant shift in the centroid position as a function of color. The galaxy number density peak is seen to move closer to the position of the mass peak for redder color intervals, and this again demonstrates that the galaxies with the most evolved stellar populations are concentrated in the densest environments. The centroid shift may indicate that A1705 (catalog ) is dynamically young. A secondary mass peak situated about $`4\mathrm{}`$ toward the northeast from the main peak appears to be associated with galaxies at the cluster redshift (the first five panels of Figure 4 appear to show an extension toward northeast from the main galaxy overdensity peak associated with A1705 (catalog )). Close to this secondary mass peak, and possibly contributing to it, is the $`z=0.112`$ IRAS galaxy F13121+7315 (catalog ) (at $`\alpha =13^\mathrm{h}13^\mathrm{m}32.0^\mathrm{s},\delta =72\mathrm{°}59\mathrm{}11\mathrm{}`$), associated with the AGN X-ray source RX J1313.5+7259 (catalog ).

The two lowest rows of the galaxy number density plots in Figure 4 shows a strong peak southwest of the location of the main cluster center. This is a galaxy concentration at a redshift well beyond A1705 (catalog ). Figure 3.2 also shows this strong overdensity in the distribution of the reddest ($`VI>2.2`$) galaxies in the field. A strong concentration of galaxies with similar colors is centered on what appears to be a distant giant elliptical galaxy located at $`\alpha =13^\mathrm{h}12^\mathrm{m}31.1^\mathrm{s},\delta =72\mathrm{°}50\mathrm{}54\mathrm{}`$. A peak in the mass distribution, which we have named WL 1312.5+7252 (catalog ), is situated less than $`1\mathrm{}`$ away from this galaxy. Its peak value is $`\kappa =0.143`$. A rough cluster redshift estimate can be derived for this peak from the median color ($`VI=2.45`$) of the early-type cluster galaxy sequence seen in a color-magnitude diagram. Using our empirically calibrated color-redshift relation for early-type cluster galaxies (see Figure 47 of Paper I), we estimate the redshift of the cluster associated with WL 1312.5+7252 (catalog ) to be $`z=0.55\pm 0.05`$. This may be a slight underestimate of the true redshift, since cluster ellipticals at $`z0.5`$ are less evolved and intrinsically slightly bluer than their counterparts at lower redshifts (Lubin 1996).

A curved blue object, apparently a gravitationally lensed arc, is situated at $`\alpha =13^\mathrm{h}12^\mathrm{m}32.6^\mathrm{s},\delta =72\mathrm{°}50\mathrm{}55\mathrm{}`$, $`7\mathrm{}`$ east of the distant giant elliptical galaxy (see Figure 3.2). This arc is quite bright (at $`I=20.90`$, $`V=21.39`$), and it is worth noting that the comparatively bright arc magnitude places WL 1312.5+7252 (catalog ) among the very small number of clusters at $`z>0.5`$ that are known to produce strongly lensed arcs that are sufficiently bright for spectroscopic studies. The only other such examples known to us are CL2236-04 (catalog ) (Melnick et al. 1993) and RCS 0224-0002 (catalog ) (Gladders, Yee, & Ellingson 2002).

If this arc is indeed a strongly lensed galaxy, the separation between the arc and the giant elliptical galaxy can be used to make a rough estimate of the mass of WL 1312.5+7252 (catalog ).

Assuming a SIS-type mass distribution centered on a giant elliptical galaxy at redshift $`z_l=0.6`$, and assuming that the lensed arc is situated at the Einstein radius $`\theta _E`$, we estimate a velocity dispersion $`\sigma _{\mathrm{arc}}=c\sqrt{(\theta _Ew_s)/(4\pi w_{ls})}900`$km s-1 for $`z_s=`$ 1 and $`\sigma _{\mathrm{arc}}700`$km s-1 for $`z_s=`$ 2. The numbers above are calculated for the case of an Einstein - de Sitter Universe, where $`w=1(1+z)^{1/2}`$. The estimated velocity dispersion of $`1150`$km s-1 for WL 1312.5+7252 (catalog ) (see Table 1) comes from a SIS-model fit to the tangential weak shear averaged in radial bins centered on the mass peak. This high value and the high aperture mass value $`M_{\mathrm{ap}}`$ provides evidence for a relatively low arc redshift. The mass-to-light ratio listed in Table 1 is measured within $`R_{\mathrm{ap}}`$ and is based on galaxies in the color range $`2.0<VI<2.8`$. The lower limit of this color interval was chosen to avoid foreground galaxies associated with A1705 (catalog ), and the upper limit was chosen to avoid background galaxies that are significantly redder than the red sequence of WL 1312.5+7252 (catalog ). The $`M/L_B`$ value is consistent with a normal luminous cluster, and assuming a redshift-dependence $`L_B(z+1)`$ for the luminosity evolution, the present-day value would be $`M/L_B=(268\pm 84)h`$.

The only pointed X-ray observation of WL 1312.5+7252 (catalog ) is a 3.8 ksec EINSTEIN IPC exposure. In the IPC image there is diffuse emission from the main A1705 cluster extending in the direction of WL 1312.5+7252 (catalog ). Also, there is a hint of diffuse emission extending towards the mass structure southeast of the main cluster. The bright X-ray point source, RX J1313.5+7259 (catalog ), located northeast of the main cluster masks any signs of diffuse emission that might be associated with the feature in the mass map in this direction.

If WL 1312.5+7252 (catalog ) is at a redshift of $`z=0.55`$, within a radius of $`r_{500}`$ it is at least a factor three X-ray under-luminous relative to the expectation from the mass-X-ray luminosity relation of Reiprich and Böhringer (2002). Again, given the scatter in the mass-luminosity relation of an order of magnitude, WL 1312.5+7252 (catalog ) could be a rather normal cluster by X-ray standards.

In conclusion, information from galaxy photometry, (probable) strong lensing and weak lensing all indicate the presence of a new massive, high-redshift cluster in the field of A1705 (catalog ). Although the cluster is not detected in existing X-ray data, the X-ray constraints on the cluster are consistent with the expected cluster X-ray luminosity, but places the cluster at the faint end of the observed distribution of X-ray luminosities for the given mass. It would of course be of great interest to accurately measure the redshift of this cluster, and of the arc, to obtain definite confirmation of the lensed nature of the arc and to increase the accuracy of mass estimates of WL 1312.5+7252 (catalog ) from strong and weak lensing. Furthermore, Chandra or XMM-Newton observations are crucial to determine the hot gas content of the cluster and to independently map the mass distribution.

### 3.3 The A1722 Field

The cluster center of A1722 (catalog ) is dominated by a single cD galaxy. A $`17\mathrm{}`$ long, thin blue gravitationally lensed arc is visible $`18\mathrm{}`$ southwest of the cD. As shown in the mass map (Figure 3.2), there is a second peak in the mass distribution, WL 1320.4+6959 (catalog ), which is almost as significant as the mass peak corresponding to A1722 (catalog ) (its peak value is $`\kappa =0.166`$). Figure 3.3 shows the distribution of red ($`VI>2.2`$) galaxies in the field. A concentration of galaxies is seen close to the position of WL 1320.4+6959 (catalog ), but with a larger spread in galaxy colors than the WL 1312.5+7252 (catalog ) cluster in the field of A1705 (catalog ). There appears to be a general overdensity of red galaxies in the southern part of the A1722 (catalog ) field, and there are several galaxy concentrations along a filamentary structure running along the southern part of the field. The same general appearance is evident in the panels corresponding to $`VI>2.0`$ in Figure 9.

The apparent early-type galaxy sequence associated with the galaxy concentration around WL 1320.4+6959 (catalog ) is less well-defined than for WL 1312.5+7252 (catalog ), but there appears to be a concentration of galaxies with colors around $`VI=2.25`$, corresponding to an estimated redshift of $`z=0.45\pm 0.05`$. A corresponding strong density peak is seen for galaxies in the range $`1.9<VI<2.5`$ in Figure 9. As shown in Figure 10, there is an offset of $`1\mathrm{}`$ between this galaxy density peak and the

WL 1312.5+7252 (catalog ) mass peak. A red galaxy with an extended envelope is located within the innermost contour of the galaxy density peak, and appears to be the central galaxy of the cluster. The two last panels of Figure 9 show a moderate number density peak of very red ($`2.65<VI<3.1`$) galaxies located even closer to the WL 1320.4+6959 (catalog ) mass peak. A color of $`VI=2.9`$ places the galaxies in this peak at $`z1`$ or higher, but any cluster at this redshift would have to be unrealistically massive to produce the observed weak lensing signal, since most of the source galaxies used in the lensing analysis are situated at lower redshifts. We therefore believe that the clustering of galaxies at $`z0.45`$ is chiefly responsible for producing the peak in projected mass density distribution. The first five panels in Figure 9 ($`1.15<VI<2.05`$) also show moderate galaxy density peaks close to WL 1320.4+6959 (catalog ). Thus, it is possible that a projection of multiple structures along the line of sight, some of them possibly associated with A1722 (catalog ), are contributing to the WL 1320.4+6959 (catalog ) mass peak. Alternatively, we may just be seeing galaxies of different colors and morphological types at a single redshift. In this case, the significant scatter seen in the red cluster sequence suggests that the cluster is less dense and has a larger fraction of late-type galaxies than WL 1312.5+7252 (catalog ), which would be expected for a cluster which is still forming. Spectroscopic determinations of galaxy redshifts in the area of WL 1320.4+6959 (catalog ) would be required to conclusively distinguish between these two scenarios.

From an ASCA SIS 66 ksec exposure, the upper limit on the X-ray luminosity from WL 1320.4+6959 (catalog ) is a factor two below the expectation from the mass-luminosity relation of Reiprich & Böhringer (2002). However, the upper limit on the X-ray luminosity is uncertain by a factor two due to the uncertain subtraction of X-ray emission from the main cluster and a bright point source (seen in an archival 27 ksec ROSAT HRI exposure) southwest of the main cluster. Although WL 1320.4+6959 (catalog ) is X-ray underluminous, existing X-ray data are consistent with a cluster at redshift of $`z0.45`$.

The mass-to-light ratio listed in Table 1 is measured within $`R_{\mathrm{ap}}`$ and is based on galaxies in the color range $`2.0<VI<2.6`$. The lower limit of this color interval was chosen to avoid foreground galaxies associated with A1722 (catalog ) (see Figure 9), and the upper limit was chosen to avoid background galaxies that are significantly redder than the apparent red sequence associated with WL 1320.4+6959 (catalog ). The $`M/L_B`$ value is consistent with a normal luminous cluster, and assuming $`L_B(z+1)`$ luminosity evolution, the present-day value would be $`M/L_B=(429\pm 171)h`$.

## 4 Discussion

Using weak gravitational lensing data, we have in the previous paragraphs identified three prominent (projected) mass concentrations representing new galaxy cluster (or sub-cluster) candidates. We use two-color, $`V`$\- and $`I`$-band photometry to investigate the nature of these structures and make rough redshift estimates.

The first mass concentration, WL 1017.3+5931 (catalog ), is the most enigmatic of these objects. The morphology of the mass peak and the associated X-ray emission suggest a possible association with the nearby $`z=0.29`$ cluster A959 (catalog ), but we find no overdensity of early-type cluster galaxies at the location of WL 1017.3+5931 (catalog ). There is also no strong evidence for clustering of galaxies at any other redshift at this position. This object remains a good candidate for an optically “dark” (sub-)cluster. The upper limit on the X-ray luminosity leaves the possibility that it is an X-ray (underluminous) cluster at or beyond the redshift of A959 (catalog ). Deep X-ray data of this system woulde be particularly interesting in order to accurately measure its hot gas content.

The second mass concentration, WL 1312.5+7252 (catalog ), appears to be associated with a rich cluster at $`z0.55`$ which also acts as a strong lens. This structure is associated with a prominent overdensity of red galaxies and has a prominent red galaxy sequence. Thus, it is likely to constitute a single strong physical overdensity of dark matter and galaxies, rather than being caused by a line-of-sight projection of lesser structures. Its $`M/L_B`$ value is similar to the typical values of clusters selected from baryonic tracers. The conservative upper limit on the X-ray luminosity of WL 1312.5+7252 (catalog ) is consistent with a cluster at $`z0.55`$ with the mass determined from weak lensing. However, WL 1312.5+7252 (catalog ) could well be an X-ray dark mass concentration, containing only a small amount of hot gas.

The third mass concentration, WL 1320.4+6959 (catalog ), is associated with an overdensity of galaxies at an estimated redshift of $`z0.45`$, but these galaxies show a larger spread in $`VI`$ color than the galaxies in the WL 1312.5+7252 (catalog ) cluster. However, the upper limit on its X-ray emission does not allow us to rule out that WL 1320.4+6959 (catalog ) is a single cluster at $`z0.45`$ which is X-ray underluminous and possibly X-ray dark. It is thus unclear whether this density peak represents a chance superposition of objects at different redshifts or whether it represents a single cluster at $`z0.45`$ which is still forming and is not yet virialized. Weinberg & Kamionkowski (2002) estimate the abundances of non-virialized, X-ray underluminous protoclusters that will be detectable in weak lensing surveys, and they find that such systems are likely to be detected in surveys comparable to ours. The $`M/L_B`$ value we find is typical for a normal optically luminous cluster.

At this point we may speculate whether it is significant that we find three mass concentrations with $`\kappa >0.15`$ based on 1.0 deg2 of imaging of fields containing massive clusters, while Wilson et al. (2001) find no such objects from similar data (1.5 deg2 of imaging) of blank fields. We also recall that the two optically dark mass concentrations previously discovered by Erben et al. (2000) and Umetsu & Futamase (2000) were both found in the fields of massive clusters. In the case of WL 1312.5+7252 (catalog ) we are seeing a mass concentration which is clearly physically unrelated to the nearby Abell cluster. For WL 1320.4+6959 (catalog ) it is hard to draw an equally firm conclusion, since, as noted in § 3.3, there may be a significant contribution to the projected mass density from structures at the redshift of A1722 (catalog ).

The photometric data for WL 1017.3+5931 (catalog ) do not provide any strong evidence for an association with nearby A959 (catalog ), but the morphology of the mass peak does suggest a possible physical link between the two objects. In a recent study, White, van Waerbeke & Mackey (2002) show that significant peaks in the projected density distribution, resembling clusters with virial masses of $`10^{14}3\times 10^{14}h^1M_{}`$ may be generated by line-of-sight projections of multiple correlated structures with $`M<10^{14}h^1M_{}`$. It is possible that WL 1017.3+5931 (catalog ) is such an object, generated by a superposition of $`10^{13}h^1M_{}`$ halos within the overdensity associated with A959 (catalog ). This scenario would naturally explain the low X-ray luminosity, but does not predict an unusually high $`M/L_B`$ value or the lack of an associated peak in the galaxy density distribution.

Clearer answers to the nature of weak lensing-detected mass concentrations may come soon from systematic cluster searches in the deep wide-field imaging data sets currently used for measurements of “cosmic shear”. The sky area collectively probed by such surveys to similar depth is now at least an order of magnitude larger than the sky area we study here (see e.g., van Waerbeke et al. 2001).

Our results also demonstrate the usefulness of multi-color photometry and color slicing techniques when interpreting results from weak lensing cluster searches. At redshifts $`z0.5`$ and higher, even rich clusters do not represent strong galaxy density enhancements in single-passband imaging data. In the case of WL 1320.4+6959 (catalog ), an optical counterpart to the mass concentration could not be identified before data in a second passband had been obtained. It is also clear that X-ray data from either Chandra or XMM-Newton are required in order to detect or tightly constrain the hot gas content of the detected (sub-)clusters.

The methodology for cluster searches may be further refined in the future by developing a more objective and quantitative search algorithm that combines weak gravitational lensing information with e.g, the cluster-red-sequence method of Gladders & Yee and/or X-ray data. The X-ray data would be particularly useful for separating line-of-sight superpositions of less massive objects (which would be a significant source of noise and bias for both optical and weak lensing data; see e.g., Hoekstra 2001; White et al. 2002) from genuine deep potential wells. This would greatly improve the power of such surveys to constrain cosmological models.

We thank the anonymous referee for suggestions that have improved our work and its presentation. We also thank Harald Ebeling, Jens Hjorth, Henk Hoekstra, Gerry Luppino, Somak Raychaudhury, and Gillian Wilson for useful comments and discussions. We thank the staff of the University of Hawaii 2.24m telescope and the Nordic Optical Telescope for support during our observing runs. HD gratefully acknowledges support from a doctoral research fellowship awarded by the Research Council of Norway. HD and PBL thank the Research Council of Norway for travel support. KP acknowledges support from the Danish National Research Council. This research has made use of the NASA/IPAC Extragalactic Database (NED) which is operated by the Jet Propulsion Laboratory, California Institute of Technology, under contract with the National Aeronautics and Space Administration. This research has made use of data obtained from the High Energy Astrophysics Science Archive Research Center (HEASARC), provided by NASA’s Goddard Space Flight Center.


text:
# Hamiltonian fractals and chaotic scattering of passive particles by a topographical vortex and an alternating current

## I Introduction

Large-scale structures such as vortices, background and tidal currents are the main constituents in geophysical flows that strongly influence all the transport and mixing phenomena. Transport of heat and mass plays a crucial role in dynamics of the ocean and atmosphere. To study the transport and mixing processes the Lagrangian representation is convenient to adopt. The motion of a test particle (a tracer) satisfies the differential equation

$$\frac{d𝐫}{dt}=𝐯(𝐫,t),$$
(1)
where $`𝐫=(x,y)`$ is the position of the particle and $`𝐯`$ represents the incompressible Eulerian velocity field. It has been recognized that particle trajectories even in simple flows can be chaotic in the sense that nearby trajectories separate exponentially in time. While in a three-dimensional flow it is possible even if the flow is time independent A65 , Lagrangian chaos or chaotic advection in a two-dimensional flow may occur only in time-dependent flows A84 ; O89 ; Z91 ; C91 .

In this paper we study the dynamics of passive particles (which take on the velocity of the flow very rapidly and do not influence the flow) in a simple two-dimensional open flow composed of a fixed point vortex and a background current with a periodic component. This model is inspired by a very interesting natural phenomena, topographical vortices over mountains, to be found in the ocean and atmosphere K83 ; Z95 . In the deep sea, they have been found in different regions of the World Ocean with the help of buoys of neutral buoyancy Col . A laboratory prototype of a topographical vortex is a cylindric anticyclonic vortex to be found by Taylor T23  in homogeneous fluid over an underwater obstacle. In the homogeneous ocean the topographical vortices are cylindric whereas they are of cone-like form in the stratified ocean.

We adopt, of course, an oversimplified model with a point vortex in an attempt to catch and quantify the main features of the tracer dynamics in the flow from the dynamical standpoint. The fixed point vortex is embedded in a background planar flow composed of steady and time-dependent components. If the flow satisfies the incompressibility condition, $`\mathrm{div}𝐯=0`$, the evolution equation (1) can be written in the Hamiltonian form

$$\begin{array}{c}\frac{dx}{dt}=v_x(x,y,t)=\frac{\mathrm{\Psi }}{y},\hfill \\ \frac{dy}{dt}=v_y(x,y,t)=\frac{\mathrm{\Psi }}{x},\hfill \end{array}$$
(2)
where the streamfunction $`\mathrm{\Psi }(x,y,t)`$ plays the role of a Hamiltonian. Thus, the configuration space of an advected particle is a Hamiltonian phase space of an associated dynamical system. It means that coordinate-space trajectories of tracers coincide with phase-space trajectories, and more importantly, the configuration space contains all the typical structures of Hamiltonian phase space including islands of regular motion, islands around islands, cantori, stochastic layers, etc. (for a recent review of Hamiltonian chaos see, for example Z98 ). This prominent non-uniformity of the phase space implies anomalous transport of tracers, including Lévy flights, their dynamical trapping and fractal properties, stickiness to islands boundaries, etc. These phenomena have been studied theoretically in the field of chaotic advection in a large number of systems (see, for example, papers A84 ; C91 ; R90 ; Z88 ; P99 ; EA ; BB01 ; JZ  which represent a small part of relevant publications) and observed in laboratory experiments SW93 ; SK96 ; ST98 ; CC87 .

In our model flow, advected particles from the inflow region enter the region where the fixed point vortex is located and then wash out to the outflow region. While the tracer trajectories in the inflow and outflow regions, where the influence of the vortex is negligibly small, are simple as they just follow the streamlines of the background current, the trajectories in an influence zone of the vortex, called the mixing region in the following, are typically chaotic BP01 . So, the problem of tracer dynamics in the open flow we address in this paper is, in fact, a problem of chaotic scattering Ott93 ; JZ ; SK96 .

The paper is organized as follows. In Section 2, we introduce the model streamfunction and the equations of motion of advected particles, consider the integrable time-independent version of the model flow and prove analytically transversal intersections of stable and unstable manifolds of the saddle point producing homoclinic chaos in the time-dependent system. A numerically simulated image of the homoclinic structure and a typical Poincaré section of the particles trajectories in the mixing region provide a general picture of chaotic advection. In Section 3, we report our main results. We observe and discuss the phenomenon of stickiness of tracers trajectories to the boundaries of the vortex core and islands of regular motion which act as dynamical traps for tracers. We show that distribution of the dynamical traps over the space of initial tracer’s positions is fractal-like, and the trapping time distribution demonstrates initially exponential decay followed by a power-law decay at the probability distribution function (PDF) tail with a characteristic exponent $`\gamma 2`$. We compute also dependence of tracer’s trapping time on their initial positions and find it to be a typical fractal scattering function with an uncountable number of singularities. The mechanism of generating the fractal is shown to resemble the famous Cantor-set generation with the Hausdorff fractal dimension of the scattering function to be equal to $`d1.84`$. Finally, in Section 4, we discuss some possible applications of the results.

## II Model flow

We consider a point fixed vortex embedded in a planar flow of an ideal incompressible fluid with a stationary and periodic components. The respective dimensionless streamfunction

$$\mathrm{\Psi }=\mathrm{ln}\sqrt{x^2+y^2}+x(\epsilon +\xi \mathrm{sin}\tau )=\mathrm{\Psi }_0+\xi \mathrm{\Psi }_1,$$
(3)
generates the Lagrange equations of motion of a passive particle in the flow BP01

$$\begin{array}{c}\dot{x}=\frac{y}{x^2+y^2},\hfill \\ \dot{y}=\frac{x}{x^2+y^2}+\epsilon +\xi \mathrm{sin}\tau ,\hfill \end{array}$$
(4)
where dot denotes differentiation with respect to dimensionless time $`\tau `$, $`\xi `$ and $`\epsilon `$ are the normalized velocities of a particle in the stationary and periodic components of the current flowing in the $`y`$ direction from the south to the north. Without perturbation, $`\xi =0`$, the phase portrait of the dynamical system (4) consists of a collection of finite and infinite trajectories (streamlines) separated by a loop passing through a saddle point with the coordinates $`(1/\epsilon ;\mathrm{\hspace{0.17em}0})`$. In the polar coordinates $`x=\rho \mathrm{cos}\phi `$ and $`y=\rho \mathrm{sin}\phi `$, the unperturbed equations can be solved in quadratures

$$\epsilon d\tau =\left[1\left(\frac{E\mathrm{ln}\rho }{\epsilon \rho }\right)^2\right]^{1/2}d\rho ,$$
(5)
where $`E=\epsilon \rho \mathrm{cos}\phi +\mathrm{ln}\rho `$ is the conserved energy. Depending on initial positions and the values of the control parameters, particles either are trapped inside the separatrix loop and move along closed streamlines or move along infinite streamlines outside the loop. In the integrable system, the stable manifold of the saddle point, along which particles move towards the saddle, and its unstable manifold, along which they move outward the saddle, coincide. Taking into account the separatrix value of the energy $`E_s=1\mathrm{ln}\epsilon `$, one can easily derive the travel time for a particle between two points with the coordinates $`(\rho _0,\rho )`$ nearby the separatrix

$$T(\rho _0,\rho )=\underset{\rho _0}{\overset{\rho }{}}\frac{d\rho ^{}}{\sqrt{\epsilon ^2\left({\displaystyle \frac{E_s+\delta \mathrm{ln}\rho ^{}}{\rho ^{}}}\right)^2}},$$
(6)
where $`\delta `$ is a small deviation of the energy from its separatrix value.

Under perturbation, the stationary saddle point becomes a saddle periodic motion which is represented on a Poincaré section with the dimensionless period $`2\pi `$ by a stationary point. In the extended phase space $`(x,y,\tau )`$, there exists a stable (unstable) manifold of the saddle point with trajectories approaching the saddle periodic motion when $`\tau \pm \mathrm{}`$. Under a typical perturbation, stable and unstable branches of the unperturbed separatrix loop split and transverse each other on the Poincaré section infinitely many times producing a complicated homoclinic structure. It can be analytically proved with the help of the Poincaré-Melnikov integral

$$I(\alpha )=\underset{\mathrm{}}{\overset{\mathrm{}}{}}\{\mathrm{\Psi }_0,\mathrm{\Psi }_1\}[x_s(\tau \alpha ),y_s(\tau \alpha )]𝑑\tau ,$$
(7)
where $`\{\mathrm{\Psi }_0,\mathrm{\Psi }_1\}`$ is the Poisson bracket, $`x_s`$ and $`y_s`$ are the separatrix solutions of the unperturbed problem parameterized by a real number $`\alpha `$. With the streamfunction (3) we get the function

$$I(\alpha )=\mathrm{sin}\alpha (\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑\tau \dot{x}_s\mathrm{cos}\tau )\mathrm{cos}\alpha (\underset{\mathrm{}}{\overset{\mathrm{}}{}}𝑑\tau \dot{x}_s\mathrm{sin}\tau ),$$
(8)
that obviously has an infinite number of simple zeroes if $`\dot{x}_s(\tau )0`$. Transversal intersections of stable and unstable manifolds of a hyperbolic point has been analytically proven with two-dimensional flows under almost arbitrary nontrivial perturbation K95 .

A simplified image of such intersections is given in FIG. 1 where a fragment of the Poincaré section with a small perturbation $`\xi =0.01`$ in the time moments $`2\pi n(n=0,\mathrm{\hspace{0.17em}1},\mathrm{\hspace{0.17em}2},\mathrm{})`$ for a bunch of trajectories, starting from a neighborhood of the saddle point, is shown. It is a “seed” of Hamiltonian chaos that arises under an arbitrary small perturbation around the unperturbed separatrix loop.

A phase portrait of chaotic advection in our simple flow is shown in FIG. 2. The Poincaré section of passive particle trajectories in FIG. 2a reveals a non-uniformity of the phase space that is typical for Hamiltonian systems with $`3/2`$ degree of freedom. Fluid around the point vortex placed at ($`x=0,y=0`$) forms a coherent structure, a vortex core, which is filled with regular, almost elliptic orbits. Another large coherent structure is seen to the west from the vortex core, it is filled with regular orbits as well. A magnification of the northern part of this long island is shown in FIG. 2b which demonstrates a chain of small islands along its border line. Further magnification of one of these small islands shown in the inset of FIG. 2b reveals chains of smaller islands and so on. A zone with invariant KAM tori embedded in a stochastic sea, is called the mixing region. While passive particle trajectories just follow streamlines outside the mixing region (in the far downstream and upstream regions) they can be chaotic inside of it.

## III Dynamical traps and fractals

Chaotic advection of passive particles in a time-dependent open flow is a typical problem of chaotic scattering. The southern background current with a comparatively small alternating component transports particles along the streamlines into the mixing region, where the particle trajectories can be very complicated and typically chaotic, and then washes out most of them to the north where they again follow the streamlines. We show in this section that in spite of the simplicity of our flow there exists a surprisingly rich variety of different types of trajectories inside the mixing region. Particles with small variations of initial positions are trapped there for trapping times with a broad spectrum extended from a particle transit time with the average velocity of the current to infinity. The dynamical trapping cannot be characterized only by the trapping time $`T`$, the ways by which advected particles move in the mixing region are very different.

Passive particles are placed outside the mixing region on the line segment with $`y_0=6`$ and different values of $`x_{0i}`$. We compute the time when the particles reach the line $`y=6`$ and positions of the particles, $`x_i(\tau )`$ and $`y_i(\tau )`$. The control parameters are chosen to be $`\epsilon =0.5`$ and $`\xi =0.1`$. In FIG. 3 we show the motion of a particle with the initial position $`(x_0=4.358533,y_0=6)`$ trapped in the mixing region for the comparatively long time $`T115`$. In the upper panel of the figure, the dependence of the $`x`$ coordinate of the particle on time $`\tau `$ (a) and the particle’s trajectory (b) are shown. The Poincaré plot of the particle positions at the moments of time $`\tau =2\pi n`$ (black circles) is shown in FIG. 3c. The unperturbed separatrix and the Poincaré section for a collection of particles initially placed inside the mixing region are shown for reference. The particle moves quite uniformly without any preference to visit special zones in the mixing region.

The second particle to be placed initially close to the first one at the point $`(x_0=4.358034,y_0=6)`$ demonstrates cardinally different behavior shown in FIG. 4. It spends in the mixing region much more time, $`T1581`$, than the first particle does. The figure clearly demonstrates the phenomenon of stickiness well-known in chaotic Hamiltonian dynamics D90 ; ZSW ; Z98 . In contrary to the first particle, the second one spends most of time in the mixing region nearby the boundaries of the vortex core and the long island. Motion of the particle in these areas is almost regular with zero maximal Lyapunov exponent. In Hamiltonian systems, there may exist partially broken invariant curves (cantori) in a vicinity of island boundaries which prevent particles to quit the area and act as a kind of permeable barriers. Both these reasons result in stickiness of the islands boundaries that is seen in FIG. 4.

To study the trapping time distribution we have carried out a number of numerical experiments. A trapping map referred to initial positions of $`N=22801`$ passive tracers, distributed homogeneously on a rectangular grid $`(x_0[5,\mathrm{\hspace{0.17em}0}];y_0[6,\mathrm{\hspace{0.17em}1}])`$, is shown in FIG. 5.

Color modulates values of the trapping time $`T`$. The inset shows the respective zone of this map magnified by the factor of $`30`$ on the $`x`$ axis and factor of $`17`$ on the $`y`$ axis. A patchiness exists for all the scales revealing a fractal-like structure of the underlying phase space. Tracers started initially in the areas of the upstream region, corresponding to black patches on the map, will be trapped in the mixing region for longer times. The inset demonstrates that there exist tracers with very close initial positions whose trapping times differ largely. To quantify this set of dynamical traps with a wide range of trapping times we compute the respective distribution function $`N(T)`$ for $`N=110011`$ tracers (positioned initially on a rectangular grid $`(x_0[4.65,4.35];y_0[6.05,5.95])`$. The double logarithmic plot in FIG. 6 demonstrates initially exponential decay followed by a power-law decay at the PDF tail, $`N(T)T^\gamma `$, with the characteristic exponent $`\gamma 2`$. While a Poissonian distribution is expected for strong chaotic mixing a power-law dependence indicates the presence of dynamical traps in the flow. Particles trapped for longer times give an enhanced contribution to the PDF tail. Different choices of initial positions of the tracers visiting the mixing region do not change significantly the trapping-time distribution.

The well-known manifestation of chaotic scattering is fractal-like scattering functions with an uncountable number of singularities Ott93 ; JZ ; SK96 . We take the trapping time to be the scattering function, $`T(x_0)`$, and compute it carefully for $`N=10^4`$ particles distributed initially on the line segment with $`y_0=6`$ and $`x_0[4.65,4.35]`$. The results in FIG. 7 demonstrate a typical scattering function $`T(x_0)`$ with an uncountable number of singularities which are unresolved in principle. The inset in FIG. 7a shows a zoom by the factor of 20 of one of these singularity zones. Successive magnifications confirm a self-similarity of the function with increasing values of the trapping times $`T`$. FIG. 4 gives an example of a trajectory with very long trapping time $`T1581`$. Remind, please, that we measure time in units of the period of the alternating current. The trapping map in FIG. 5 provides a two-dimensional image of the fractal structure of the scattering function $`T(x_0,y_0)`$. As we have shown above, sticky boundaries of islands of regular motion, embedded in a stochastic sea, act as a kind of dynamical traps providing all the spectrum of values of the trapping times up to infinity. Namely these dynamical traps are the ultimate reason for appearing fractals in our simple model flow. To give an insight into a mechanism of generating the fractal we plot in FIG. 7b segments of the initial string with $`10^6`$ tracers which are trapped in the mixing region after $`n`$ rotations around the fixed point vortex. After each rotation, a portion of tracers is washed out in the downstream region with $`y6`$. This process resembles the mechanism of generating the famous Cantor set but is more complicated.

In conclusion, we compute the length of the curve $`T(x_0)`$ shown in FIG. 7a

$$L(1/N)=\underset{i=1}{\overset{N}{}}|T_{i+1}T_i|,$$
(9)
as a function of the inverse number of tracers, $`N^1`$, distributed initially on the same line segment with $`y_0=6`$ and $`x_0[4.75,4.25]`$ in the upstream region of the flow. It is a decreasing (in average) function with a negative slope whose double lagarithmic plot is presented in FIG. 8. The Hausdorff dimension of the fractal-like scattering function is computed to be $`d1.84`$.

## IV Conclusion and discussion

We have shown that chaotic advection of passive particles by the open flow, composed of a fixed point vortex and a background current with a time-dependent component, has a fractal nature. Appearance of fractals and dynamical traps in geophysical flows with topographical vortices should fluence strongly transport and mixing of heat and mass. Instead of homogeneous mixing, one expects a highly-structural, hierarchical fractal distribution of advected particles that is of great interest in oceanography and atmospheric sciences. It should be emphasized that the fractals and dynamical traps are not specific for our particularly simple model. They are typical in low-dimensional Hamiltonian systems and should exist in more realistic models with elliptic topographical vortices and different boundary conditions considered, for example, in KK .

If tracers are not passive but active (chemical reagents or biological species), their activity is superimposed on chaotic mixing in open background flows. Typically, there exists a local coupling between a spatial character of the activity and the spatial character of the fluid with coherent structures, dynamical traps and fractals. Say in chemistry, molecules will react if they are in a close contact with each other. So, the imperfect mixing should largely increases the chemical or biological activity in some regions of the flow and suppress it in other regions. In marine biology, the paradoxical observational fact of coexisting thousands of phytoplankton species in spite of a limited number of resources for which they compete has been known for a long time H61 . Chaotic mixing of active particles in open nonstationary flows in the ocean with small-scale spatial inhomogeneities and fractal patterns can provide a natural explanation of the plankton paradox K00 . Homogeneous and complete mixing leads to an equilibrium state in the system of competing species resulting in the survival of the most fitted species, whereas chaotic mixing generates, typically, a persistent nonequilibrium state in the system providing the coexistence of much more species than the number of ecological niches they occupy.

## V Acknowledgments

This work is supported by the Russian Federal Program “World Ocean” in the framework of the Project “Modelling variability of hydrophysical fields” and by the Russian Foundation for Basic Research under project numbers 02-02-17796 and 02-02-06841. We gratefully acknowledge useful discussions with V.F. Kozlov.