text:
# Multiresolution Analysis of Substructure in Dark Matter Halos

## 1 Introduction

The hierarchical clustering hypothesis provides an attractive paradigm for the formation of structure in a universe dominated by cold dark matter. Small-scale objects form first and merge to yield systems of increasing size. This highly non-linear process has been studied extensively using N-body simulations with particular attention paid to the survival of subhalos once a merger event has occurred. Early results suggested that substructures (i.e., subhalos within halos) are erased efficiently (White 1976; Frenk et al. 1988). This so-called overmerging problem plagued investigations of galaxy cluster formation since it lead to the conclusion that the constituent galaxies do not survive. However, recent high resolution simulations (Ghigna et al. 1998, Klypin et al. 1999) together with analytic work (Moore, Katz, & Lake 1996) demonstrated that the overmerging problem was due entirely to the poor mass and spatial resolution of early simulations. Indeed, Moore et al. (1999) found that on galactic scales, simulated halos may have too much substructure. Their high resolution simulation of a $`10^{12}M_{}`$ (i.e., Milky-Way sized) halo revealed over 500 $`M10^8M_{}`$ subhalos, a factor of 50 greater than the number of visible satellites in the Milky Way.

An essential element in the analysis of cosmological N-body simulations is the identification of physical structures, namely halos and subhalos. There are now a number of algorithms available to do this such as friends-of-friends (FOF; Davis et al. 1985), DENMAX (Bertschinger & Gelb 1991; Gelb & Bertschinger 1994), and SKID (Governato et al. 1997). The FOF algorithm identifies structures by linking all particle pairs separated by less than a user-supplied distance known as the linking parameter. For a particular choice of linking parameter, the algorithm produces a unique list of structures. DENMAX, and the closely related algorthm SKID, are based on contour surfaces of the three-dimensional density field. Each local maximum of the density field is assumed to correspond to the center of a halo or subhalo: All particles interior to the last closed contour surrounding a given maximum are assigned to the corresponding halo.

Two challenges, brought to the fore by the dramatic improvements in the mass and force resolution of present-day simulations, now confront these methods: (1) identification of subhalos within the dense regions of a parent halo; and (2) analysis of multiple levels of substructure. A given simulation particle may be a member of a small clump that is, in turn, gravitationally bound to a larger (sub)halo. Likewise, the subhalo may be gravitationally bound to a galactic or cluster-sized halo. Such a particle is most accurately described as being a member of three distinct structures. Thus, any scheme which assigns a given particle to at most one structure cannot hope to capture the hierarchical nature of dark matter halos.

Early incarnations of the FOF algorithm relied on a single linking length. If the linking length is set to be too large, subhalos in the inner parts of large halos are missed. With a small linking length, on the other hand, the algorithm picks out substructure but loses information on large scales. These problems can be avoided if one uses a “hierarchical” version of FOF (Klypin et al. 1999) wherein the algorithm is run several times with different linking lengths. DENMAX and SKID are best suited to finding substructure since they locate all maxima in the density field. The hierarchical nature of clustering can be studied by applying smoothing filters to the density field and rerunning the algorithm.

In this paper, we describe a multiresolution analysis (MRA) that handles, in a natural way, the multiple levels structure found in cosmological N-body simulations. MRA refers to a general class of tools that provide a simple hierarchical representation of a signal, the signal, in our case, being the density field in a simulation. At each resolution, the analysis picks out the details of the signal at a characteristic scale. Thus, MRA can be thought of as a “mathematical microscope”: Coarse resolution (low magnification) probes large-scale structures in the signal while fine resolution (high magnification) probes small-scale structures.

We employ a specific MRA that is based on the wavelet transform known as the $`\mathrm{𝑎𝑙𝑔𝑜𝑟𝑖𝑡ℎ𝑚𝑒}\stackrel{`}{a}\mathrm{𝑡𝑟𝑜𝑢𝑠}`$. A wavelet transform allows one to analyze a signal simultaneously in scale and position. The transform accomplishes this task by convolving the signal with a special type of window function known as a wavelet. Wavelets must have compact support and integrate to zero. The wavelet transform therefore probes local properties of the density field and is insensitive to a mean background. The wavelets used in a particular implementation of the transform are chosen to be translations and dilations of a single prototype known as the mother wavelet. A low resolution analysis of a signal is achieved by using large-scale versions of the mother wavelet: high resolution is achieved by using small-scale versions. Thus, wavelet analysis fits naturally within the framework of MRA (Mallat 1989a, 1989b).

The wavelet transform can provide a complete description of a density field found in a cosmological N-body simulation. Segmentation analysis produces, from the results of the wavelet transform, a catalog of structures and substructures. Dynamical information can then be used to cull, from this catalog, unbound associations of particles.

The algorithm described in this paper was first used in the context of cosmological N-body simulations by Lega et al. (1995) and Lega et al. (1996) who were interested in the morphology of large scale structure in various cosmological models (e.g., hot vs. cold dark matter). It has also been used by Gambera et al. (1997) to analyze observational data for the Coma cluster. We adapt this algorithm to the specific task of sorting out multiple levels of substructure that arise within individual dark matter halos. Our implementation of the algorithm is outlined in Section 2. Results from the analysis of a cluster-sized dark matter halo are presented in Section 3. In Section 4 we describe our somewhat unorthodox procedure for identifying bound systems of particles. Statistics of the subhalos are discussed in Section 5 where we also present an example of a region of the simulation that exhibits multiple levels of substructure.

## 2 Multiresolution Analysis of an N-body Simulation

In this section, we describe the two steps which comprise the heart of the algorithm: (1) calculation of the wavelet transform of a density field; and (2) application of segmentation analysis to identify individual physical structures.

### 2.1 From the Continuous Wavelet Transform to the $`\mathrm{𝐚𝐥𝐠𝐨𝐫𝐢𝐭𝐡𝐦𝐞}\stackrel{\mathbf{`}}{𝐚}\mathrm{𝐭𝐫𝐨𝐮𝐬}`$

The wavelet transform of a given signal represents the detail in the signal as a function of position and scale. Operationally, the wavelet transform is computed by convolving the signal with a window function known as a wavelet. For a particular wavelet transform, the set of analyzing wavelets is constructed by applying translations and dilations to the mother wavelet. Thus, the mother wavelet characterizes the transform.

In a continuous wavelet transform (CWT), the set of analyzing wavelets is chosen to include all possible translations and dilations of the mother wavelet. If the signal depends on one variable (say position) the CWT will be a continuous function of two variables, position and scale. It follows that the information contained in a CWT is highly redundant. By contrast, the discrete wavelet transform (DWT) utilizes a countable set of analyzing wavelets constructed by taking discrete translations and dilations of the mother wavelet. Provided certain conditions are met, the redundancy in the wavelet transform can be eliminated. In other words, the set of analyzing wavelets in a DWT can be chosen to form a complete, orthonormal basis for a large class of continuous functions.

For large data sets, the number of computations required in the standard DWT becomes exceedingly large. There are now a number of efficient algorithms which make it possible to perform DWTs on large data sets. These algorithms are related to the standard DWT much in the way the fast Fourier transform is related to the standard Fourier transform. The $`\mathrm{𝑎𝑙𝑔𝑜𝑟𝑖𝑡ℎ𝑚𝑒}\stackrel{`}{a}\mathrm{𝑡𝑟𝑜𝑢𝑠}`$(Holschneider et al. 1989, Dutilleux 1989) is an example of such an algorithm. Our description of this algorithm is preceded by a general discussion of continuous and discrete wavelet transforms.

The CWT of a one-dimensional function $`f(x)`$ is given by

$$w(\lambda ,a)=_{\mathrm{}}^{\mathrm{}}f(x)\psi ^{}(\lambda ,a;x)𝑑x$$
(1)
where the parameters $`\lambda `$ and $`a`$ characterize the scale and position of the analyzing wavelet $`\psi (\lambda ,a;x)`$ and $`\psi ^{}`$ is the complex conjugate of $`\psi `$. $`\psi (\lambda ,a;x)`$ is constructed by applying a rescaling (by $`\lambda `$) and a translation (by $`a`$) to the user-supplied mother wavelet $`\psi _0`$:

$$\psi (\lambda ,a;x)=\frac{1}{\sqrt{\lambda }}\psi _0\left(\frac{xa}{\lambda }\right).$$
(2)
The original function $`f`$ can be reconstructed via the inverse transform

$$f(x)=w(\lambda ,a)\psi (\lambda ,a;x)𝑑\lambda 𝑑a$$
(3)
provided the mother wavelet satisfies the admissibility condition

$$\left|\widehat{\psi }_0(k)\right|^2\frac{dk}{|k|}<\mathrm{}$$
(4)
where $`\widehat{\psi }_0`$ is the Fourier transform of $`\psi _0`$. From this condition it follows that $`\widehat{\psi }(0)𝑑x\psi (x)=0`$ and therefore the wavelet must be oscillatory, i.e., wavelike in nature.

A CWT resembles a windowed Fourier transform (WFT) in the sense that both operations utilize a set of window functions or filters to analyze a signal simultaneously in scale and position (see, for example, Kaiser 1994; Burrows, Gopinath, & Guo 1998). However, the width of the window function in the CWT scales with $`\lambda `$ whereas the width of the window function in a WFT is fixed. Thus, the CWT is scale-independent: Rescaling lengths in the data set leaves the CWT unchanged modulo a shift along the scale axis.

The CWT defined in Eq. (1) maps a function of one variable into a function of two variables. Thus, the resultant transform, $`w(\lambda ,a)`$, contains redundant information (i.e., $`\psi (\lambda ,a;x)`$ forms an overcomplete basis). Consider instead the DWT and associated reconstruction formula:

$$w_{n,i}=_{\mathrm{}}^{\mathrm{}}\psi _0^{}\left(\frac{xa_i}{\lambda _n}\right)f(x)$$
(5)
$$f(x)=\underset{n,i}{}\psi _0\left(\frac{xa_i}{\lambda _n}\right)w_{n,i}$$
(6)
where $`i`$ and $`n`$ are integers, $`a_i=i\mathrm{\Delta }a`$, and $`\lambda _n=\sigma ^nDeltaa`$. $`\sigma `$ and $`\mathrm{\Delta }a`$ are constants representing respectively the ratio between different levels $`n`$ of the DWT and the sampling interval. The $`w_{n,i}^{}{}_{}{}^{}s`$ are referred to as the wavelet amplitude coefficients (WACs). In what follows, we take $`n=1`$ to correspond to the smallest scale (highest resolution) level in the transform. Increasing $`n`$ corresponds to dilating the wavelet and thus degrading the resolution. In practice, the signal is sampled at regular intervals — for the case at hand, the sampling interval corresponds to the grid spacing of the discretized density field — and it is natural to choose $`\mathrm{\Delta }a`$ to be equal to the grid spacing. The integral in Eq. (5) is then replaced by a sum over $`j`$ where $`xx_j=j\mathrm{\Delta }a`$.

The goal in wavelet analysis is to find a set of analyzing wavelets such that any square-integrable function can be expressed by Eq. (6). Many examples of suitable wavelet bases have been discovered. However, for large data sets, the computation requirements become prohibitively large. Consider a signal that consists of $`N`$ data points. If the wavelet has compact support so that it takes $`m`$ computations to perform a single convolution at level $`n=1`$, it will take $`m\times \sigma ^{(n1)}`$ computations to perform a convolution at level $`n`$ and therefore $`N\times m\times \sigma ^{(n1)}`$ computations in total for level $`n`$ of the transform. In principle, we may consider wavelets whose scale is comparable to the size of the data set, i.e., $`m\times \sigma ^{n_{\mathrm{max}}1}N`$ where $`n_{\mathrm{max}}`$ refers to the highest level of the transform. The number of computations for this level will be $`O(N^2)`$. Since $`n_{\mathrm{max}}\mathrm{log}_\sigma \left(N/m\right)`$ the complexity of the transform will be $`N^2\mathrm{log}_\sigma \left(N/m\right)`$.

The $`\mathrm{𝑎𝑙𝑔𝑜𝑟𝑖𝑡ℎ𝑚𝑒}\stackrel{`}{a}\mathrm{𝑡𝑟𝑜𝑢𝑠}`$is an efficient scheme for performing discrete wavelet transforms. The algorithm assumes a scaling parameter $`\sigma =2`$. Wavelet coefficients are obtained recursively in that the WACs at level $`n`$ are calculated by taking an appropriate sum over WACs at level $`n1`$. Thus, the number of computations necessary for the highest level of the transform is at most $`O(N)`$ (or more accurately, $`O(N\mathrm{ln}N)`$). Since the algorithm is an example of the application of wavelets to MRA (Mallat 1989a, 1989b) we first interject a few general remarks about MRA.

Consider a function $`f(x)`$ sampled at a fixed regular interval $`\mathrm{\Delta }a`$. This highest resolution sample of $`f`$, which we refer to as $`f_0`$, is constructed by convolving $`f`$ with a window function $`\varphi `$:

$$f_{0,i}=\underset{j}{}\varphi \left(x_ix_j\right)f\left(x_j\right).$$
(7)
(The second subscript on $`f_0`$ refers to position, i.e., $`x_i=i\mathrm{\Delta }a`$, etc.). Smoother representations of $`f`$ can be constructed with filters $`\varphi _n`$ that are dilated versions of $`\varphi `$:

$$f_{n,i}=\underset{j}{}\varphi _n\left(x_ix_j\right)f\left(x_j\right)$$
(8)
where

$$\varphi _n\left(x\right)=\frac{1}{2^n}\varphi \left(\frac{x}{2^n}\right).$$
(9)
For this reason, $`\varphi `$ is referred to as the scaling function.

The “details” $`w_{n,i}`$ of the signal correspond to the part of $`f`$ that is removed as one degrades the resolution from one level to the next. Schematically, we have

$$f_0=w_1+w_2+\mathrm{}w_{N1}+f_N$$
(10)
where $`f_N`$ is the lowest resolution version of $`f`$ and

$`w_{n,i}`$ $``$ $`f_{n1,i}f_{n,i}`$ (11)
$`=`$ $`{\displaystyle \underset{j}{}}\left(\varphi _{n1}\left(x_ix_j\right)\varphi _n\left(x_ix_j\right)\right)f_{x_j}.`$
Comparing Eqs. (6) and (11) we see that the mother wavelet can be calculated by a simple subtraction

$$\psi (x)=\varphi (x)\frac{1}{2}\varphi \left(\frac{x}{2}\right).$$
(12)
For MRA, we require that the scaling function satisfies the recursion relation

$$\varphi \left(\frac{x}{2}\right)=\underset{j}{}h_j\varphi \left(xx_j\right)$$
(13)
where $`h_j`$ are constants (Mallat 1989a, 1989b; see, also Kaiser 1994). It is this equation that enables us to calculate the WACs recursively.

The $`\mathrm{𝑎𝑙𝑔𝑜𝑟𝑖𝑡ℎ𝑚𝑒}\stackrel{`}{a}\mathrm{𝑡𝑟𝑜𝑢𝑠}`$(Holschneider et al. 1989, Dutilleux 1989) is an MRA that has a number of attractive features. The algorithm is computationally efficient and easy to program. For example, reconstruction of the original function involves a simple sum over scale at each position. In three dimensions, the algorithm is approximately isotropic. Finally, the wavelet coefficients at each level are calculated for all points of the highest resolution (interval $`\mathrm{\Delta }a`$) grid and therefore the shapes and sizes of structures are well-determined. By contrast, the original MRA developed by Mallat (1989a, 1989b) is computationally intensive, the reconstruction formula is complicated, and the transform in 2 or more dimensions is not isotropic. Moreover, wavelet coefficients are calculated on a decimated grid (in going from level $`n`$ to level $`n+1`$ the number of wavelet coefficients is reduced by a factor of 2 in one dimension or 8 in three dimensions) and therefore the representation of larger structures becomes rather crude. A detailed comparison of Mallat’s MRA and the $`\mathrm{𝑎𝑙𝑔𝑜𝑟𝑖𝑡ℎ𝑚𝑒}\stackrel{`}{a}\mathrm{𝑡𝑟𝑜𝑢𝑠}`$can be found in Shensa (1992).

Following Lega et al. (1995), we choose the scaling function to be the cubic B-spline:

$$\varphi (x)=\frac{1}{12}\left(|x2|^34|x1|^3+6|x|^34|x+1|^3+|x+2|^3\right)$$
(14)
where distances are measured in units of $`\mathrm{\Delta }a`$. Note that $`\varphi (x)=0`$ for $`|x|2`$. The recursion formula, Eq. (13), takes the form

$$\varphi \left(\frac{x}{2}\right)=\frac{1}{8}\underset{j=2}{\overset{2}{}}C_{2j}^4\varphi \left(xx_j\right)$$
(15)
where $`C_n^mn!/m!(nm)!`$ are the usual binomial coefficients. Eq. (15) can be verified by a straightforward but tedious calculation. The mother wavelet, calculated from Eq. (12) and shown in Figure 1, is similar to the popular Mexican Hat wavelet (the Laplacian of a Gaussian). Explicit formulae used in the calculation of the wavelet coefficients are given in Appendix A.

As discussed above, $`𝑑x\psi _0(x)=0`$. In addition, the first moment of $`\psi _0`$ also vanishes: $`𝑑xx\psi _0(x)=0`$. It follows that the leading contribution to the wavelet transform of $`f(x)`$ is from terms quadratic in $`x`$: The wavelet transform is evidently unaffected by a constant or linearly varying background. It is this property of the wavelet transform that makes it ideal for identifying structures in cosmological simulations.

In three dimensions, the wavelet is constructed from the scaling function $`\mathrm{\Phi }`$, assumed to be separable in Cartesian coordinates:

$$\mathrm{\Phi }\left(𝐫\right)\varphi (x)\varphi (y)\varphi (z).$$
(16)
Once again, explicit formulae for the WACs are given in Appendix A. The wavelets constructed in this manner are quasi-isotropic, i.e., the profile of the mother wavelet depends on the direction in space, though it is the same along any three mutually orthogonal directions. Moreover, the differences along non-orthogonal directions are relatively minor as can be seen in Figure 1 where we compare the profile of the wavelet along one of the three original axes with that along one of the diagonals $`x=\pm y=\pm z`$.

### 2.2 Application of the $`\mathrm{𝐀𝐥𝐠𝐨𝐫𝐢𝐭𝐡𝐦𝐞}\stackrel{\mathbf{`}}{𝐀}\mathrm{𝐓𝐫𝐨𝐮𝐬}`$to N-body Simulations

The first step in applying the $`\mathrm{𝑎𝑙𝑔𝑜𝑟𝑖𝑡ℎ𝑚𝑒}\stackrel{`}{a}\mathrm{𝑡𝑟𝑜𝑢𝑠}`$to an N-body simulation is to discretize the density field by approximating the density at each point of a three-dimensional Cartesian grid. We use the nearest-grid-point scheme frequently employed in particle-mesh simulations (Hockney & Eastwood 1988). Simply put, the density at a particular grid point is given by the sum of the masses of all particles within a cell centered on that point divided by the volume of the cell.

For the analysis presented in this paper was use a grid $`128^3`$ in size, the maximum one allowed given memory constraints of the available computers. The number of grid points is roughly a factor of four greater than the number of particles in the simulation and therefore we do not expect that a finer mesh will lead to an improvement in the results. Indeed, the highest resolution wavelets are able to pick out associations of very small numbers of particles.

The $`n=1`$ wavelet has a width of approximately five cells. Thus, only wavelets up to $`n=5`$ (80 cells across) fit into our grid. The $`n=5`$ level essentially provides a map of the gross features of the halo and therefore useful information on subhalos is contained in levels $`n=14`$.

With the discretized representation of the density field in hand, the WACs can be calculated according to the formulae in Appendix A.

### 2.3 Segmentation Analysis

The multiresolution analysis described above provides a complete description of the (discretized) density field obtained in an N-body simulation. The next step in the analysis is to develop a prescription for identifying physical structures. For this step, we take, as a working definition of structure, an association of particles in space that is unlikely to have occurred by chance. This definition includes both gravitationally bound clumps and the debris of systems that have been tidally disrupted by the parent halo. In the next section, we describe how one can distinguish between these two possibilities.

In the wavelet description, structures are identified as connected regions where the WACs rise above a predetermined threshold. The introduction of a user-supplied threshold is a feature common to all clump-detection algorithms. For example, in the FOF algorithm, the threshold corresponds to the linking parameter. In an MRA, a different threshold is chosen for each level in the analysis. Following Lega et al. (1995) we calculate the wavelet transform for a random distribution of particles where the number of particles and “simulation volume” are chosen to be the same as in the actual simulation to be analyzed. At each level, the threshold is chosen to be five times the rms value of the wavelet coefficients.

At first glance, the $`5\sigma `$ threshold would seem to eliminate virtually all chance associations of particles. However, the statistical fluctuations in the particle number are different in a highly nonuniform halo simulation than in a random distribution. Indeed, the wavelet transform maps not only substructure but also the large scale structure of the main halo. In any case, application of a $`5\sigma `$ threshold is just the first step in producing a catalog of substructures: Dynamical information is also used to remove unbound (supposedly chance) associations of particles. In the next section we present results for different choices of the threshold which demonstrate that the algorithm, taken as a whole, is relatively robust.

The thresholding procedure leads to a binary representation of the wavelet transform with a value $`1`$ ($`0`$) assigned to gridpoints where the WAC is above (below) the threshold. The next step it to identify individual structures. Initially, lists of structures (connected region of above-threshold WACs) are generated for each level in the transform. These structures are identified by a method known as segmentation analysis (Rosenfeld 1969; Lega et al. 1995), the details of which are presented in Appendix B.

Typically, the larger structures in a simulation lead to significant WACs at multiple levels in the transform. The spherically averaged radial density profiles for halos and subhalos are well-approximated by piecewise power-law functions of radius. Therefore no single scale can characterize a halo. The upshot is that the combined list of structures from all levels of a transform will contain numerous redundancies. For example, a single smooth subhalo (i.e., one devoid of substructure) that is identified at the $`n=4`$ level will almost certainly also be represented by entries in the $`n=1,\mathrm{\hspace{0.17em}2},`$ and $`3`$ levels. These unwanted redundances are eliminated as follows: For a given structure identified at the top level of the transform ($`n=4`$ for the case at hand), we remove a single entry, namely the one closest to the center of the $`n=4`$ structure, from each of the lists of structures at the lower levels. The procedure is repeated for the $`n=3`$ and $`n=2`$ lists. The end result is an irreducible list of subhalos and sub-subhalos where each distinct object is represented by one and only one entry.

The final list of objects will contain not only gravitationally bound systems but also unbound associations of particles (e.g., tidal streams) and density enhancements that occur because of fluctuations in the particle distribution. Dynamical information enables us to distinguish among these possibilities. Our prescription for identifying bound systems is described in Section 4.

## 3 Substructure in a Cluster-Sized Halo

Our substructure-finding algorithm is applied to the N-body simulations of cluster-sized halos that appeared in Dubinski (1996) and Ghigna et al. (1998). These simulations were performed using parallelized versions of the Barnes-Hut treecode (Barnes & Hut 1986). In each case, the main halo exhibits a high degree of substructure with hundreds of subhalos.

In what follows, we describe the results from our analysis of the Gigna et al. (1998) simulation (Figure 2). A more extensive discussion of these results as well as the results for the Dubinski (1996) simulation can be found in Seymour (2000). In the Gigna et al. (1998) simulation, the virial radius of the main halo is $`2\mathrm{Mpc}`$. There are approximately $`5\times 10^5`$ particles inside this radius corresponding to a mass of $`4.6\times 10^{14}M_{}`$.

We begin by constructing a $`128^3`$ pixelized map of the density field. The DWT for levels $`n=14`$ are then calculated using the formulae in Appendix A. Two-dimensional projections of the wavelet coefficients, constructed by selecting the maximum wavelet coefficient along each ‘line-of-sight’, are shown in Figure 3. This figure illustrates the manner in which substructure on different scales is captured in the different levels of the DWT. Note that the central object and many of the larger satellites register strong signals in all four levels illustrating the way in which the DWT captures the internal structure of these systems.

The results of the segmentation analysis are shown in Figure 4 where circles, superimposed on the projected particle distribution, are drawn for the structures detected at each level. The radii of the circles corresponds to the size of the region found to be above threshold (i.e., proportional to the cube root of the number of pixels above threshold). The fact that most of the large subhalos are detected in multiple levels in the transform is reflected in the appearance of concentric circles. When this redundancy is removed (see Section 2.3) 773, 223, 72, and 5 subhalos remain in levels $`14`$ respectively.

## 4 Identifying Bound Systems

The MRA described above identifies structures in the density field but does not distinguish between particles that are members of a subhalo and interlopers, i.e., particles that are passing through the subhalo but are not associated with it in any true dynamical sense. In addition, some of the structures that we have identified are extremely small and may in fact be chance associations that arise from Poisson fluctuations in the density field. Dynamical information makes it possible to eliminate interlopers and also cull the list of structures of unbound particle groups.

The standard working definition for a bound structure in an N-body simulation is “the largest set of particles that are mutually gravitationally bound” (see, for example, Bertschinger & Gelb 1991; Gigna et al 1998). With this definition, the procedure for identifying bound structures given the results of our wavelet and segmentation analysis is as follows: (1) For each entry in our list of structures, calculate its radius and center-of-mass. For the radius, we take $`r=\left(3m/4\pi \right)^{1/3}\mathrm{\Delta }a`$ where $`m`$ is the number of cells associated with the structure, i.e., connected cells with WACs above the threshold. Likewise, the center-of-mass is calculated from the WACs that rise above the threshold. (2) Next calculate the total energy (kinetic plus potential) for each particle within a distance $`d<r`$ from the structure’s center of mass. Note that the kinetic energy is calculated in the rest frame of the particles. (3) The highest energy particle is identified and if its energy is greater than zero, it is removed. One then returns to steps (2) and (3) with the proviso that in calculating the potential energy and center of mass frame, only the remaining particles are used. The procedure is repeated until the highest energy particle that remains has negative energy.

In excluding unbound particles when calculating the gravitational potential, one underestimates the (negative) potential energy of the remaining particles and may therefore inadvertently remove bound particles from a subhalo. An alternative procedure is to recalculate the center-of-mass velocity but not the potential energies each time an unbound particle is removed. In this way, all of the particles in the neighborhood of the clump contribute to the gravitational potential, as they should!

Figure 5 shows an example of a subhalo that is identified as a bound system by the second method (interlopers included in the potential) but rejected by the first method (i.e., no mutually gravitationally bound (sub)system of particles was found). The spatial distribution of particles for the three Cartesian projections are shown in Figure 5a. Note that only a small fraction of the particles (26 out of 532) are identified as being members of the gravitationally bound system thus explaining why, in this case, the background particles are so important. To check that we have indeed found a bona fide structure, we show, in Figure 5b, the velocity space distribution of particles. The bound clump forms a tight group in velocity space: Their velocity dispersion is much smaller than that of all of the particles in this region ($`\sigma _{\mathrm{bound}}/\sigma _{\mathrm{all}}=0.19`$). We can estimate, from Poisson statistics, the probability that such a velocity space clump will occur by chance. The probability of finding a single particle in the velocity-space region of the clump is $`p=\left(\sigma _{\mathrm{bound}}/\sigma _{\mathrm{all}}\right)^30.0063`$ and thus the probability of finding $`n_b=26`$ particles out of $`n_t=532`$ in this region is $`e^xx^{n_b}/n_b!`$ where $`x=pn_t`$. Thus, the probability of finding $`n_b`$ out of $`n_t`$ particles with velocity dispersion $`<\sigma _b`$ in any region of velocity space is

$$P=\frac{1}{p}\frac{e^xx^{n_b}}{n_b!}.$$
(17)
In the example of Figure 5, $`P3.7\times 10^{13}`$.

For many of the systems found by the second method, the case that they are bona fide bound structures is not nearly so strong. In particular, for dense regions of the halo, it is not so unlikely to find, by chance, small groups of particles that appear to be bound. Thus, we need some additional set of criteria to further cull, from the list of structures, chance associations of particles. The criteria we choose are as follows: For each bound system we calculate the probability $`P`$ as above and also the fraction of local particles $`F=n_b/n_t`$ that are found to be part of the bound system. A bound clump is deemed to be genuine if either $`P<10^4`$ or $`F>0.75`$. The first condition applies to bound systems moving through dense regions of the halo (as in Figure 5) while the second condition applies to relatively isolated systems in which $`n_bn_t`$ and $`\sigma _b\sigma _t`$. Our criteria are somewhat arbitrary and other choices are possible. For example, Gigna et al. (1998) only consider groups of 16 or more particles. This particle number cutoff may be overly conservative as our algorithm identifies numerous systems of fewer than 16 particles that are clearly clustered in both velocity and configuration space. Nevertheless, our results at the low-mass end of the subhalo distribution may be suspect. Clearly, the relevant numerical experiment is to run a set of simulations with idential initial conditions but with different numbers of particles.

The results for both methods of finding bound clumps are summarized in Table 1. With method 1, nearly three quarters of the level 1 structures and over one half of the level 2 structures are removed as being “unbound” while with method two, less than half of the level 1 structures and less than one third of the level two structures are removed. The additional constraint eliminates (as potentially being chance associations) an additional 100 level 1 structures.

| $`\sigma =5`$ | | | | |
| --- | --- | --- | --- | --- |
| $`n`$ | $`SA`$ | $`B1`$ | $`B2`$ | $`\overline{B2}`$ |
| 1 | 773 | 177 | 433 | 328 |
| 2 | 223 | 119 | 158 | 158 |
| 3 | 72 | 72 | 72 | 72 |
| 4 | 5 | 5 | 5 | 5 |

We have repeated the procedure outlined above using both $`4\sigma `$ and $`6\sigma `$ thresholds and the results are shown in Tables 2 and 3. As expected, the number of structures identified by our wavelet and segmentation analysis using a $`4\sigma `$ threshold is higher than the number identified with a $`5\sigma `$ threshold and likewise the number identified using a $`6\sigma `$ threshold is lower. However, the unbinding procedure seems to identify the greatest number of structures when a $`5\sigma `$ threshold is used. The unbinding procedure begins with a region specified by the segmentation analysis where it can search for a bound object. If the region is much larger than the clump, and the clump center-of-mass velocity is much different from the center-of-mass velocity of the background, the unbinding procedure may miss the clump altogether. These results illustrate a truism common to all substructure-finding algorithms: On an object by object basis near the mass resolution of the simulation subhalo identification is often ambiguous. That is, one is likely to miss some bound clumps while including, in the catalog of structures, some chance associations. However, as we demonstrate below, the statistical properties of the subhalo population (namely fraction of main halo in bound clumps) appears to be insensitive to the choice of threshold.

| $`\sigma =4`$ | | | |
| --- | --- | --- | --- |
| $`n`$ | $`SA`$ | $`B2`$ | $`\overline{B2}`$ |
| 1 | 820 | 335 | 244 |
| 2 | 212 | 200 | 187 |
| 3 | 77 | 77 | 77 |
| 4 | 10 | 10 | 10 |

| $`\sigma =6`$ | | | |
| --- | --- | --- | --- |
| $`n`$ | $`SA`$ | $`B2`$ | $`\overline{B2}`$ |
| 1 | 439 | 313 | 240 |
| 2 | 134 | 134 | 134 |
| 3 | 63 | 63 | 63 |
| 4 | 3 | 3 | 3 |

We conclude this section by considering a region of the simulation $`800\mathrm{kpc}`$ on a side ($`1.5\%`$ of the simulation volume) centered on a large subhalo that has been detected at level 4 of the DWT. The upper left panel of Figure 6a displays all of the particles in this region while the upper right panel displays only those particles gravitationally bound to the main subhalo. We have also identified twelve smaller subhalos within the region defined by the main subhalo. These are shown in the lower left and lower right panels (levels 2 and 1 of the analysis respectively). (Many of the clumps that appear in the upper left panel are in the foreground or background of the main subhalo and therefore are not shown.) One of these small subhalos, indicated by the arrow, is gravitationally bound to the main subhalo. In addition, a second subhalo appears to be associated dynamically with the main subhalo. The remaining ten subhalos are interlopers as is evident from Figure 6b where the 12 small subhalos are plotted in velocity space. The circle represents the velocity dispersion of the large subhalo. Note that the relative velocities of the subclumps are typically much higher than their internal velocities.

The gravitationally bound subhalo is an example of a third level of substructure in that its constituents can be regarded as members of three distinct systems. Examples of this type are rare, a result that is not surprising given that the systems which constitute the second level of substructure (the large subhalos orbiting the main halo) are consist of only a few thousand particles and therefore subject to the purely numerical overmerging problem.

## 5 Substructure Statistics

Our subhalo detection algorithm has identified over 500 subhalos and it is therefore possible to study, in a statistical sense, their characteristics. In Figure 7, for example, the rms internal velocity dispersion of the clumps are plotted as a function of their position within the main halo. As noted in Ghigna et al. 1998, a wide range of clump velocities are found throughout the main halo though there does appear to be a trend toward larger internal velocities for subclumps closer to the center of the main halo. This result may be an indication that subhalos are heated by the tidal field of the main halo. In addition, there appears to be a slight enhancement of subclumps with large rms velocities at $`r1.5\mathrm{Mpc}`$. These subhalos are probably associated with the large subhalo evident in Figure 2 at $`x1.2\mathrm{Mpc},y0.1\mathrm{Mpc}`$. In Figure 8, the masses of the subhalos (shown in terms of particle number) are plotted as a function of rms velocity. One finds a reasonably tight correlation between $`M`$ and $`v_{\mathrm{rms}}`$ with $`Mv_{\mathrm{rms}}^\alpha `$ where $`\alpha 3`$. The correlation appears to hold over nearly 3 orders of magnitude in mass. This result is to be compared with Ghigna et al. 1998 (see their Figure 20) who find $`\alpha 33.4`$ over approximately two orders of magnitude in mass.

In Figure 8, there appears to be a distinct population of objects that have, for fixed mass, velocities a factor of 3-5 greater than those in the main distribution. These objects are found preferentially in the inner regions of the main halo once again suggesting that they are heated by the tidal fields of the main halo. This result brings us to Figure 9 where we plot the virial ratio, $`2T/|W|`$ as a function of position in the main halo. $`T`$ and $`W`$ are respectively the total kinetic and potential energies of the bound system. $`T`$ is simply the sum of the kinetic energies of the particles in the system. Calculation of $`W`$ is somewhat more involved since a sum of the potential energies of the individual bound particles double counts for pairs of bound particles. The formula for $`W`$ is:

$$W=\underset{i\{B\}}{}\left(\frac{1}{2}\underset{j\{B\};ji}{}\frac{m_im_j}{|𝐫_𝐢𝐫_𝐣|}+\underset{j\{U\}}{}\frac{m_im_j}{|𝐫_𝐢𝐫_𝐣|}\right)$$
(18)
where $`\{B\}`$ and $`\{U\}`$ refer respectively to the sets of bound and unbound particles or interlopers. The population of subhalos exhibit a large scatter in the virial ratio about the equilibrium value of $`1`$. This result implies that the halos are, for the most part, not fully virialized, perhaps because that are constantly being disturbed (if not disrupted) by the tidal field of the main halo. Indeed, there appears to be a trend toward a larger dispersion in the virial ratio toward the center of the main halo.

Finally, in Figure 10, we plot the fraction of the virial mass contained in subhalos of a given size. We find that in total, $`12\%`$ of the halo is bound in subhalos in good agreement with results from Gigna et al. (1998). The figure includes results for different choices of the threshold. The fact that the results do not change much is an indication of the robust nature of the algorithm.

## 6 Conclusions

Recent advances in numerical cosmology lead to the conclusion that the more particles one uses in an N-body simulation, the more substructure one finds in dark matter halos. This result is consistent with the picture that halos form through hierarchical clustering since an increase in the number of simulation particles translates into an increase in the dynamic range of the hierarchy accessible to the simulation. What is perhaps surprising is the extent to which substructure survives.

The identification of substructure inside dark matter halos is a nontrivial task. While large subhalos are easily located by eye, small subsystems, and especially ones inside dense regions of the halo, can be difficult to spot. Moreover, there can be many levels of substructure.

The substructure-finding algorithm presented in this paper is intrinsically hierarchical and therefore perfectly suited to the task at hand. The DWT analyzes the density field at different resolutions while segmentation analysis produces a catalog of structures at each level of the transform. Further analysis is required to eliminate redundances and cull the catalog of unbound systems.

One advantage of a wavelet-based analysis is that it is insensitive to a constant or linearly varying background density. This feature makes it possible to pick out structures such as the one shown in Figure 5. Note that the density enhancement for this structure is only $`5\%`$. To pick such a structure out using FOF would require an extremely finely tuned linking parameter.

Our results indicate that the MRA algorithm is competitive with FOF and DENMAX/SKID, the two most widely used methods. In particular, our results for the number of subhalos found, the fraction of mass in subhalos, and subhalo statistics are consistent with those of Gigna et al. (1998). A detailed comparison of the three methods for a sequence of simulations performed at increasing particle resolution is required to determine if there is any true advantage of one method over the others.

We are grateful to B. Moore and J. Dubinski for providing us with their simulations. We thank A. Babul, J. Navarro, R. Henriksen, K. Perrett, J. Stadel, and D. Stiff for useful conversations. LMW would like to thank the University of Chicago and the Canadian Institute for Theoretical Astrophysics for their hospitality during recent visits. This work was supported, in part, by the Natural Science and Engineering Research Council of Canada.

## Appendix A WACs for the $`\mathrm{𝐀𝐥𝐠𝐨𝐫𝐢𝐭𝐡𝐦𝐞}\stackrel{\mathbf{`}}{𝐀}\mathrm{𝐓𝐫𝐨𝐮𝐬}`$

In this appendix, we provide the details required to implement the $`\mathrm{𝑎𝑙𝑔𝑜𝑟𝑖𝑡ℎ𝑚𝑒}\stackrel{`}{a}\mathrm{𝑡𝑟𝑜𝑢𝑠}`$. Consider a one-dimensional function $`f(x)`$. The highest resolution sample of $`f`$ is constructed from Eq. (7) or equivalently

$$f_{0,i}=\frac{1}{6}\left(f\left(x_{i1}\right)+4f\left(x_i\right)+f\left(x_{i+1}\right)\right).$$
(A1)
$`f_1`$ can be calculated from $`f_0`$ using Eqs. (8) and (9) as follows:

$`f_{1,i}`$ $`=`$ $`{\displaystyle \frac{1}{2}}{\displaystyle \underset{j}{}}\varphi \left({\displaystyle \frac{x_ix_j}{2}}\right)f\left(x_j\right)`$ (A2)
$`=`$ $`{\displaystyle \frac{1}{16}}{\displaystyle \underset{j}{}}{\displaystyle \underset{k}{}}C_{2k}^4\varphi \left(x_ix_jx_k\right)f\left(x_j\right)`$ (A3)
$`=`$ $`{\displaystyle \frac{1}{16}}{\displaystyle \underset{k}{}}C_{2k}^4f_{0,ik}.`$ (A4)
Likewise, the coefficients corresponding to smoother samples of $`f`$ can be calculated recursively:

$$f_{n,i}=\frac{1}{16}\underset{k}{}C_{2k}^4f_{n1,m}$$
(A5)
where $`mik2^{n1}`$. The WACs are calculated by a simple subtraction (Eq. (11).

Next, we consider a three-dimensional density field $`\rho (x,y,z)`$. As in the one-dimensional case, the first step is to compute the auxiliary coefficients, here written as $`c(n,i,j,k)`$. The $`n=0`$ coefficients are computed as follows:

$$c(1,i,j,k)=\underset{i^{},j^{},k^{}=2}{\overset{2}{}}\rho (i^{},j^{},k^{})\varphi \left(i^{}i\right)\varphi \left(j^{}j\right)\varphi \left(k^{}k\right)$$
(A6)
where we abbreviate $`x_i=i\mathrm{\Delta }a`$ as $`i`$. The $`n>1`$ coefficients are given by

$$c(n,i,j,k)=\underset{i^{},j^{},k^{}=2}{\overset{2}{}}h\left(i^{}\right)h\left(j^{}\right)h\left(k^{}\right)c(n1,ii^{}2^{n1},jj^{}2^{n1},kk^{}2^{n1})$$
(A7)
where $`h(i)\frac{1}{16}C_{2i}^4`$. Finally, the wavelet coefficients are found by performing the subtraction

$$w(n,i,j,k)=c(n,i,j,k)c(n+1,i,j,k)$$
(A8)
## Appendix B Segmentation Analysis

This appendix provides an outline of the segmentation analysis routine used in the paper. Further details can be found in Rosenfeld (1969) and Lega et al. (1995). The algorithm is introduced by way of a worked example in two spatial dimensions. The extension to three dimensions is straightforward.

Figure 11a presents the results for a single level of a mock wavelet transform of a two-dimensional density field. The number in each cell represents the WAC (for convenience, rounded off to the nearest whole number). Application of a threshold between 4 and 5 leads to Figure 11b.

The core of the segmentation algorithm is the scan. The scan identifies connected segments of non-zero elements along each of the two (or three, for the actual data) Cartesian directions. Starting in upper left corner of the data array, the first scan proceeds along rows successively from top to bottom labelling each segment of $`1^{}s`$ by a unique integer. The result of the left-to-right scan is shown in Figure 11c.

Next, one scans from top to bottom along successive columns. In this scan one assigns the lowest integer contained in each segment. The result is shown in Figure 11d. For a three-dimensional data set, there is an additional to be performed.

Iterating the above procedure (i.e., alternating between horizontal and vertical scans) will eventually yield the desired result: a unique label for each connected structure. However, for large data sets and odd-shaped structures, the computation time can become prohibitively long. For a data set containing $`N^2`$ elements, each scan takes $`O(N^2)`$ operations ($`N`$ operations along each row times $`N`$ columns or vice versa). One can imagine particularly perverse examples (for example, staircase-like structures) in which one requires $`O(N^2)`$ scans, in other words, $`O(N^4)`$ computations. In three dimensions, one could in principle require $`O(N^9)`$ computations which is unreasonable for even modest values of $`N`$.

One can obviate the need for repeated scans by means of the following scheme. The scan is performed once in each of the Cartesian directions. By virture of the scanning process, each cluster of connected segments is described by a unique set of integers. That is, no cluster of connected segments contains values of another cluster. In our example, the cluster on the left contains values $`1,\mathrm{\hspace{0.17em}4},\mathrm{\hspace{0.17em}6},`$ and $`8`$ while the cluster on the right contains values $`2,\mathrm{\hspace{0.17em}7},`$ and $`11`$. Our goal is to group together the integers represented in each cluster. This step is accomplished by first constructing a dictionary as follows: For each non-zero element in the array, check all non-zero neighboring elements and list the element and its neighbor if they are different. Note that the last direction scanned (in our example, the vertical direction) need not be check since all non-zero neighbors are necessarily the same. The dictionary is then pruned for duplicates (Figure 12a).

The remaining step is to sort the dictionary into connected pairs. First, one picks out pairs that have the same number in either the first or second position. The sort for the cluster on the left of our data is shown in Figure 12b. We then relabel the second digit of those pairs selected by the sort (Figure 12c). The group is then pruned to eliminate redundancy (Figure 12d) to yield a table that identifies all of the integers for cluster 1 with the single integer 4. Figure 12e shows the result when the remaining dictionary is sorted.


format: text

text:
# The Normal Cluster Weak Lensing Survey: Mass Profiles and M/L Ratios of Eight Clusters at z=0.2

## 1 Motivation and sample selection

Everyone at this conference agrees that cluster masses are important as cosmological probes and as clues to the formation process of clusters and larger structures. Yet most of the work on cluster masses presented here has been based on the assumption that clusters are in hydrostatic or virial equilibrium, while at the same time we have seen many clusters that are clearly not in equilibrium. Gravitational lensing can play an important role here because it is one of the few methods that does not make dynamical assumptions. However, most of the clusters examined with lensing to date have been rare, spectacular clusters, and thus are not representative of clusters in general.

Therefore we are conducting a weak lensing survey of “typical” clusters, which we define here in terms of X-ray luminosity: $`5\times 10^{43}<L_X`$ (0.3—3.0 keV) $`<7\times 10^{44}`$ $`h^2`$ ergs s-1, or $`0.23.0L_X^{}`$. We imposed three selection criteria based on observing considerations: $`60<\delta <0`$ (i.e., visible from Cerro Tololo Inter-American Observatory \[CTIO\]), $`|b|>20`$, and no star brighter than $`R=10`$ lying within a 45 diameter field centered on the cluster (to minimize the area lost to bright star haloes). Finally, we imposed a redshift cut motivated by the desire for as uniform a sample as possible, at $`z0.2`$, while retaining $`10`$ clusters. This resulted in the criterion $`0.15<z<0.30`$. We applied these criteria to three catalogs: the Extended Medium Sensitivity Survey (EMSS), the Einstein pointed archive, and an early version of the Serendipitous High-redshift Archival ROSAT Cluster (SHARC) survey (the two SHARC clusters were not included in their final statistical cluster sample, but were found in the course of the SHARC survey using the methods cited). The resulting eight clusters are listed in Table 1.

## 2 Observations and data reduction

We took deep multicolor images of each cluster at the CTIO’s 4-m Blanco telescope, with the Big Throughput Camera mounted at prime focus. The BTC was designed to eliminate concerns about mass sheet degeneracy with its large field. The camera has a 2 $`\times `$ 2 array of back-illuminated 2048$`\times `$2048-pixel CCDs. With pixels subtending 0.43′′ each and large gaps between the CCDs, the total field of view comes to 35. We covered the gaps and increased the total field to 45 with large dithers between exposures. We observed in the $`B_jRI`$ photometric system, with at least $`B_J`$ and $`R`$ coverage on each cluster.

This paper is based on preliminary results from the eight clusters listed. We imaged two additional clusters, Abell 2744 and MS0508.8-4523 (Abell 3322), which fit our criteria, but no results are available for them yet. These two clusters were also re-observed with the Mosaic camera at the same telescope, and will be used to ensure that no instrumental artifacts survive the data processing. We also observed most clusters in $`I`$ and some in Johnson-Cousins $`V`$ for the sake of studying other cluster properties, but the current work considers the $`B_J`$ and $`R`$ information only. A forthcoming paper presents the observations, lensing results, and comparison with X-ray data and with simulations in great detail; future papers will study the cluster luminosity functions and alignment of cluster members.

The data reduction steps, especially those aspects related to optical distortions and object shape measurement, are essentially the same as that described elsewhere in some detail for “blank” fields observed with the BTC. The main difference is that here, we do not attempt to combine shape measurements from different filters. Since the $`R`$ images are deeper than the other filters, we simply use the shapes from the $`R`$ images. Because these fields are centered on clusters, the analysis diverges from Ref. 9 beginning with the background galaxy selection.

## 3 Background galaxy selection and lensing analysis

With clusters in the foreground, we make some effort to select background galaxies based on color. The color locus of cluster ellipticals is easily identifiable in the color-magnitude diagram, so we used only galaxies which were significantly bluer in an attempt to eliminate cluster members from the lensing analysis. A uniform criterion $`B_jR<1.5`$ was used for all cluster fields. We also imposed a magnitude cut $`22<R<26.5`$ to eliminate foreground objects and objects too faint to measure accurately. We are currently working on estimates of the residual contamination by cluster members.

Given this background galaxy sample, we use the model-independent aperture densitometry method to estimate the mass. In discrete form, the crucial equation relating the galaxy shapes to the mean surface mass density $`\overline{\mathrm{\Sigma }}`$ inside a given radius, minus the mean surface density in a control annulus outside that radius is:

$$\overline{\mathrm{\Sigma }}(r<\mathrm{r}_0)\overline{\mathrm{\Sigma }}(\mathrm{r}_0<\mathrm{r}<\mathrm{r}_{\mathrm{max}})=\mathrm{\Sigma }_\mathrm{c}\underset{\mathrm{i}}{\overset{\mathrm{N}}{}}\gamma _\mathrm{i}\mathrm{g}(\mathrm{r}_\mathrm{i})\frac{\mathrm{r}_{\mathrm{max}}}{\mathrm{r}_\mathrm{i}}^2,(\mathrm{r}_0<\mathrm{r}_\mathrm{i}<\mathrm{r}_{\mathrm{max}})$$
(1)
Here $`\mathrm{\Sigma }_c`$ is the lensing critical density, $`\gamma _i`$ is the tangential ellipticity of the $`i`$th background galaxy corrected for seeing effects, $`r_i`$ is the projected distance from that galaxy to the center of the mass distribution, g($`r`$) is a weighting function, and N is the number of background galaxies in the annulus $`(\mathrm{r}_0<\mathrm{r}<\mathrm{r}_{\mathrm{max}})`$. We found the center of the mass distribution from two-dimensional maps built up by repeatedly applying a similar equation over a grid in x and y. In most cases this coincided with the brightest cluster galaxy and the center of the X-ray emission, but in a few cases it was displaced by up to 1 from the X-ray center. $`\mathrm{\Sigma }_c`$ depends on the source redshift distribution and was calibrated for Abell 3364 by simulations with full 3-dimensional ray tracing from model galaxies with a redshift and morphology distribution consistent with the Hubble Deep Field and recent Keck spectroscopy. These simulations were in turn checked against mass densitometry of clusters with known mass.

To obtain a comparable estimate of luminosity surface density, we also measured the total light overdensity in $`R`$ band (with no color selection of objects) in the same successive apertures, with respect to the same control annulus. $`R`$ band in the observed frame is very nearly $`V`$ band in the rest frame, so we will quote M/L in rest-frame $`V`$.

## 4 Mass and M/L ratios

Cluster masses and M/L ratios within 1$`h^1`$ Mpc are summarized in Table 2. The 1$`h^1`$ Mpc radius was chosen to give reasonable signal-to-noise even for the least massive clusters in the sample. The separation between low M/L and high M/L clusters is roughly three times the measurement error. While a cluster-to-cluster scatter in M/L has been seen in the literature, it is not clear whether the differences in the literature are entirely real or may be due in part to the widely varied approaches to selection and analysis. Here an ensemble of clusters of similar X-ray luminosity, observed and analyzed in a uniform way, shows similar scatter, strengthening the case for real differences in M/L from cluster to cluster.

What might cause real differences in M/L? It has been argued that because the highest overdensities collapsed first, they now have the oldest populations and therefore the highest M/L. The highest initial overdensities should now be the most massive clusters, so that the clusters with higher M/L should also tend to be more massive. We face two problems in testing for this effect. First, our clusters were selected to lie in a narrow range of X-ray luminosity and therefore presumably in a narrow mass range. Second, the effect of passive luminosity evolution is fairly small, while our errors in the masses are about 20%. There is a hint that the more massive clusters in Table 2 have higher M/L, but no conclusion can be reached given the large error bars. A sample including both massive and low-mass clusters should be able to answer this question definitively. We note that our clusters do fall nicely on a curve of enclosed M/L versus radius for all types of systems.

Figure 1 shows the mean mass and light profiles of Abell 3364. The other clusters have similar profiles but generally have lower signal-to-noise. Mass drops more steeply than isothermal at $`r>1`$ Mpc (the apparent rise in the outermost point may be due to systematic errors from the point-spread function at this very low level of shear). In these normal clusters, mass follows light surprisingly well on all scales to which this technique is sensitive, out to at least 2$`h^1`$ Mpc; there is no evidence for an increase in M/L in cluster outskirts. This lends some legitimacy to extrapolating to lower density environments since the mass and light for these normal clusters are followed out to unprecedented small values: $`\mathrm{\Sigma }30`$ M pc-2 mass density, and $`0.01`$ L pc-2, factors of 1000 and 10,000 smaller than customarily encountered in the rich clusters studied in the past.

## 5 Estimate of $`\mathrm{\Omega }_{\mathrm{matter}}`$

There are two ways to extrapolate from cluster mass information to $`\mathrm{\Omega }_{\mathrm{matter}}`$, the mean mass density of the universe as a fraction of the critical density. Both methods assume that clusters are a fair sample of the universe, so we must use some caution in interpreting the results. However, such an assumption is more plausible with this sample of typical clusters than with previous samples which have been based on a few very rich clusters.

One method can be applied immediately: if the cluster M/L holds for the universe in general, we can simply scale the cluster M/L by the local luminosity density, which is well known, and obtain the local mass density. In practice we scale the local luminosity density back to z=0.2 using the recent luminosity evolution of galaxies. Our result is $`\mathrm{\Omega }_{\mathrm{matter}}=0.19\pm 0.03`$, but we caution that the true uncertainty lies in the M/L fair sample hypothesis rather than in the statistical error given. The previous section has shown that mass is indeed well-correlated with light in the typical cluster environment, but extremely low-density environments presumably emit less than their share of light, so this might be considered a lower limit on $`\mathrm{\Omega }_{\mathrm{matter}}`$.

The second method will be applied in the near future: if the cluster baryon fraction holds for the universe in general, we can simply divide $`\mathrm{\Omega }_b`$ by this fraction to obtain $`\mathrm{\Omega }_{\mathrm{matter}}`$. $`\mathrm{\Omega }_b`$, the mean density of baryons in the universe, is well constrained by Big Bang nucleosynthesis and so does not introduce great uncertainty into the result. However, this is generally considered an upper limit to $`\mathrm{\Omega }_{\mathrm{matter}}`$ because any census of baryons in a cluster is likely to miss baryons (in currently undetectable forms such as black holes, planets, etc.) rather than overestimate them. We don’t have any cluster baryon fraction measurements as yet, but soon the combination of Sunyaev-Zel’dovich effect measurements and lensing will provide these. For now we simply note that values in the literature yield $`\mathrm{\Omega }_{\mathrm{matter}}0.3h^2`$, which as an upper limit is a nice complement to our value derived from luminosity scaling.

## 6 Shear-selected clusters

In the two-dimensional mass maps of the cluster fields, we found significant additional mass clumps apart from the target clusters. The density of these objects is roughly 4-7 deg-2. In many cases, moderately bright galaxies can be seen clumped in these areas, indicating the possible presence of a cluster or group. In other cases, there is no obvious optical counterpart. We followed up several of the cluster candidates with CTIO 4-m spectroscopy, and confirmed that each is a cluster with a well-defined redshift, not a chance projection. In no case are they associated with the target cluster, but they have the same general range in redshift, which is not surprising considering that the background galaxy selection did not vary. These discoveries are an important step toward a shear-selected sample of clusters and will be described in a forthcoming paper.

## 7 Summary

Our weak gravitational lensing survey of ten clusters is still in progress, but several conclusions are emerging. We find a significant cluster-to-cluster scatter in M/L, despite uniform observing, reduction, and analysis procedures. Mass follows light very closely on scales from $`200h^1`$ kpc (our resolution limit) to at least 2$`h^1`$ Mpc (the limit of our field size), and to unprecedentedly low projected mass and light surface densities. This gives us some confidence in extrapolating to even larger scales. If the M/L of these clusters is representative of the universe, we can scale by the local luminosity density to get $`\mathrm{\Omega }_{\mathrm{matter}}`$ $`\stackrel{>}{}`$0.2. In the near future, cluster baryon fractions will be available from lensing plus Sunyaev-Zel’dovich effect measurements, providing an estimate of $`\mathrm{\Omega }_{\mathrm{matter}}`$ from baryon scaling arguments. The two scaling arguments really provide lower and upper limits respectively; together they will bracket $`\mathrm{\Omega }_{\mathrm{matter}}`$ if these “normal” clusters are a fair sample with respect to M/L and baryon fraction.

The 45 fields are large enough to include several several serendipitous clusters, but more importantly we discovered them on the basis of shear alone. This demonstrates the feasibility of a shear-selected sample of clusters, which would be immensely valuable in determining the cluster mass function and avoiding the biases toward emitted light which have always accompanied the study of clusters. Several surveys with this goal are already underway.

## Acknowledgments

We thank Alistair Walker, John Filhaber, David Rojas, Steve Heathcote, and the staff at CTIO for their help in making the $`BTC`$ and 4-m Blanco telescope system work in an optimal way and for their continued support. Thanks also to Bob Lee, Morley Blouke, and Pat Waddell for their engineering help on the BTC project. Initial development of the BTC was partially funded by NSF SBIR grant AST86-17058 to Lassen Research.


text:
# Strongly Non-Equilibrium Bose-Einstein Condensation in a Trapped Gas

## Abstract

We present a qualitative (and quantitative, at the level of estimates) analysis of the ordering kinetics in a strongly non-equilibrium state of a weakly interacting Bose gas, trapped with an external potential. At certain conditions, the ordering process is predicted to be even more rich than in the homogeneous case. Like in the homogeneous case, the most characteristic feature of the full-scale non-equilibrium process is the formation of superfluid turbulence.
PACS numbers: 03.75.Fi, 32.80.Pj, 67.40.-w

Kinetics of Bose-Einstein condensation (BEC) in a weakly interacting Bose gas is one of the most fundamental problems of non-equilibrium statistical mechanics. The exciting progress in the experiment on BEC in ultracold gases initiated by the pioneer works  opens up an opportunity of laboratory study of BEC kinetics. Though an important step in this direction has been already made , most non-trivial regimes of ordering kinetics have not been achieved yet.

In Refs.  (see also Ref. ) a scenario of BEC in a strongly non-equilibrium self-evolving homogeneous gas was developed. The scenario involves a number of qualitatively different stages. The evolution starts with an explosion-like wave in the wavenumber space, propagating from higher energies towards the lower ones, and resulting in a formation of a specific power-law distribution of particles over the wavenumbers . Immediately after its formation, this distribution starts to relax, the relaxation having the form of a back wave in the wavenumber space . Simultaneously, in the low-energy region the so-called coherent regime sets in that leads to the formation of the quasicondensate state . Basically, the quasicondensate state corresponds to what is known in the theory of superfluidity as the state of superfluid turbulence. It can be viewed as a condensate containing a tangle of vortex lines (plus a sharply non-equilibrium distribution of long-wave phonons) . The formation of the quasicondensate occurs very rapidly (characteristic time is much smaller than the time of the wave formation) . In contrast to it, the final stage of long-range ordering, associated with relaxing superfluid turbulence and long-wave phonons, takes a macroscopically large (the size of the system enters the answer) time .

From a common wisdom one may expect that at certain conditions strongly non-equilibrium BEC in a trapped gas should reproduce main qualitative features of the homogeneous scenario. It is not surprising thus that the explosion-like initial kinetic stage  is quite relevant to the case of the trapped gas  and is consistent with the experimental observations of Ref. . Nevertheless, full-scale experimental realization of the scenario  in a trapped gas is not guaranteed yet by just considerable deviation from the equilibrium in the initial state. The point is that the situation with the BEC kinetics in ultracold gases is determined by a competition between large number of particles and very small interaction. The latter circumstance leads to a rather large coherent length associated with the ordering process, which for a trapped gas can turn out to be just larger than the size of the condensate, in which case the character of the condensate formation will be drastically changed: Instead of quasicondensate state there will take place a build-up of a genuine condensate with the initial form pre-determined by the external potential. The theoretical analysis of the BEC kinetics in a trapped gas has been dealing mostly with this physical picture , ignoring thus the nontrivial regime associated with the superfluid turbulence formation.

The main goal of the present letter is to consider the superfluid-turbulence regime of BEC kinetics in an external potential. We will demonstrate that the necessary and sufficient condition for this regime to occur corresponds to violation of the Knudsen condition at a certain stage of self-evolution within a certain region of coordinate space. While the external potential together with the initial number and energy of particles determines the initial parameters of the quasicondensate: the size $`R_{}`$, the number of quasicondensate particles $`\stackrel{~}{N}_0`$, and typical initial separation between vortex lines in the vortex tangle, $`r_0`$, the very process of the quasicondensate formation goes in a quasi-homogeneous anti-Knudsen regime when the external potential becomes irrelevant and, in principle, may be even turned off without any considerable effect on the process of the quasicondensate formation.

In the experiments with trapped ultracold gases normally there takes place the Knudsen regime, when the free path length of a particle with the energy $`\epsilon `$, $`l_{\mathrm{free}}(\epsilon ,n_\epsilon )`$ \[which depends also on occupation numbers $`n_\epsilon `$ because of the stimulated scattering\], is much larger than the typical radius of the particle’s trajectory, $`R_\epsilon `$. For definiteness we consider a parabolic trap with all the three frequencies of the same order $`\omega _0`$. So that $`R_\epsilon \omega _0v(\epsilon )`$ \[$`v(\epsilon )`$ is the typical velocity corresponding to the energy $`\epsilon `$\], and the condition $`l_{\mathrm{free}}(\epsilon ,n_\epsilon )R_\epsilon `$ is equivalent to $`\tau _{\mathrm{coll}}(\epsilon ,n_\epsilon )\omega _01`$, where ($`\mathrm{}=1`$)

$$\tau _{\mathrm{coll}}(\epsilon ,n_\epsilon )1/m(a\epsilon n_\epsilon )^2$$
(1)
is the collisional time. Here $`a`$ is the scattering length and $`m`$ is the particle mass.

Knudsen regime is a convenient starting point for analyzing kinetics in a potential. Almost in all qualitative aspects it corresponds to an isotropic homogeneous case, since the distribution of particles depends only on the two variables, $`\epsilon `$ and $`t`$ (ergodic approximation). The scaling of the collision time remains the same, so that the main quantitative difference comes from the difference in the density of states.

During some initial period of evolution the correlations between different single-particle modes are vanishingly small. Such a regime admits a description in terms of kinetic equation. This stage is thus referred to as kinetic stage. If the deviation from equilibrium is strong enough, the evolution at the kinetic stage results in the formation of a self-similar wave in the wavenumber space propagating in an explosion-like fashion from the high-energy region (where the particles are initially distributed) towards lower energy scales. Corresponding solution of the kinetic equation has the form

$`n_\epsilon (t)=A\epsilon _0^\alpha (t)f(\epsilon /\epsilon _0(t)),tt_{},`$ (2)
$`\epsilon _0(t)=Bt_{}t^{1/2(\alpha 1)}.`$ (3)
Here $`A`$ and $`B`$ are dimensional constants depending on the initial condition and related to each other by the formula $`B=\mathrm{const}(ma^2A^2)^{1/2(\alpha 1)}`$. The dimensionless function $`f`$ is defined up to an obvious scaling freedom. The explosion character of the evolution guarantees that the wave reaches the point $`\epsilon =0`$ at some finite time moment $`t=t_{}`$ \[$`t=0`$ corresponds to the beginning of evolution\], the value of $`t_{}`$ being on the order of the collisional time at $`t=0`$. Physically, the explosion-like evolution is supported by the stimulation of the collision rate at the head of the wave, $`\epsilon \epsilon _0`$, by ever growing occupation numbers. A remarkable feature about the solution (2)-(3) is its universality: The coefficient $`A`$ is the only quantity related to a particular cooling regime or a particular form of the initial condition. Having in mind the most characteristic statement of the problem when the self-evolution starts with an essentially non-equilibrium state of $`N`$ particles in a parabolic trap of typical frequency $`\omega _0`$, with a typical occupation numbers $`\gamma >1`$, we readily estimate $`A`$ by extrapolating the solution (2)-(3) to $`t=0`$:

$$A\gamma ^{1\alpha /3}\omega _0^\alpha N^{\alpha /3}.$$
(4)
Generally speaking, the index $`\alpha `$ in (2)-(3) cannot be established from the scaling properties of the collision term of the kinetic equation, being related thus to the particular form of the latter . It is possible, however, to specify lower and upper limits for $`\alpha `$ following from the consistency of (2)-(3) with the requirement that these formulae describe an explosion-like singularization of distribution (rather than infinite-time shrinking). To this end we note that from the scale invariance of the collisional term it follows that $`f(x)`$ should behave at $`x1`$ like some power of $`x`$. Then, from the fact that at $`t=t_{}`$ the occupation numbers are finite at $`\epsilon >0`$ one concludes that

$$f(x)x^\alpha \mathrm{at}x\mathrm{}.$$
(5)
The requirement that the particle distribution does not shrink as a whole implies that the number-of-particles integral for the distribution (2), (5) is divergent at $`\epsilon \mathrm{}`$. The number of particles $`N=w(\epsilon )n_\epsilon 𝑑\epsilon `$, where the density of states $`w(\epsilon )`$ is given by

$$w(\epsilon )=\{\begin{array}{cc}Vm^{3/2}\epsilon ^{1/2}/\sqrt{2}\pi ^2\text{(homogeneous)},\hfill & \\ 2\epsilon ^2/\omega _x\omega _y\omega _z\text{(parabolic potential)}.\hfill & \end{array}$$
(6)
Here $`V`$ is the volume of a homogeneous system and $`\omega _x`$, $`\omega _y`$, and $`\omega _z`$ are the frequencies of the 3D harmonic oscillator. This immediately yields $`\alpha <3/2`$ for the homogeneous case and $`\alpha <3`$ for the case of parabolic potential. The condition $`\alpha >1`$ is necessary for $`\epsilon _0(t)`$ to approach zero at $`t=t_{}`$. So we have $`1<\alpha <3/2`$ in the homogeneous case and $`1<\alpha <3`$ in the parabolic potential. The most accurate up-to-date numeric analysis of $`\alpha `$ in the homogeneous case was performed in Ref.  with the result $`\alpha 1.24`$. In the case of parabolic potential $`\alpha 1.6`$ .

Formally, at $`t=t_{}`$ the front of the wave $`\epsilon _0(t)`$ reaches the energy $`\epsilon =0`$. However, somewhat earlier the solution (2)-(3) becomes inapplicable in the region of very small energies. In the homogeneous case this occurs due to the change of the evolution regime (in corresponding low-energy region) from weak turbulence to strong turbulence (the so-called coherent regime) . In the external potential, however, the scenario is modified. As we demonstrate below, what changes the solution (2)-(3) at a certain stage of evolution of a trapped gas is either (i) discreteness of the levels, or (ii) violation of the Knudsen regime in a certain region of the coordinate space. Apart from a cross-over region, the two circumstances are mutually exclusive: If discreteness of the levels becomes relevant, the Knudsen regime will not be violated, and if the Knudsen regime is violated, the discreteness of levels will not become relevant. Hence, the change of the evolution in a trapped gas from the self-similar wave (2)-(3) is governed by the following two parameters:

$$\xi =\tau _{\mathrm{coll}}\omega _0\text{(Knudsen parameter)},$$
(7)
$$d=\tau _{\mathrm{coll}}/w(\epsilon )\text{(discreteness parameter)}.$$
(8)
Knudsen regime corresponds to $`\xi 1`$; the discreteness of the levels is irrelevant when $`d1`$. Note that $`\tau _{\mathrm{coll}}`$ essentially depends on $`\epsilon `$ and $`n_\epsilon `$, which ultimately means that both $`\xi `$ and $`d`$ are the functions of energy and time.

First, let us consider the case (i), where the condition $`d1`$ occurs when the front of the wave $`\epsilon _0(t)`$ reaches some energy $`\epsilon _d`$. As it follows from (1)-(3), (6), the parameter $`d`$ at the front of the wave can be estimated as

$$d_{\epsilon _0}\omega _0^3\epsilon _0^{2(\alpha 2)}/ma^2A^2.$$
(9)
From Eq.(9) we see that $`d_{\epsilon _0}`$ decreases with decreasing $`\epsilon _0`$ only if $`\alpha <2`$, which means that the case (i) is possible only under this condition. \[In this respect the result $`\alpha 1.6`$  is of qualitative importance!\] Substituting $`A`$ from (4) into (9) and setting $`d_{\epsilon _0}1`$, we find

$$\epsilon _d\omega _0p^{1/2(2\alpha )},$$
(10)
$$p=\omega _0ma^2\left(N^{\alpha /3}\gamma ^{1\alpha /3}\right)^2.$$
(11)
Here $`p`$ is the main dimensionless parameter which we will encounter several times. By definition, in the case (i) the value of $`\epsilon _d`$ Eq.(10) should be greater than $`\omega _0`$, which implies $`p1`$. Below we will see that $`p1`$ is not only a necessary, but also a sufficient condition for the case (i) to occur: Under this condition the case (ii) is impossible.

To describe what happens after the front $`\epsilon _0(t)`$ reaches the energy scale $`\epsilon _d`$, we need to make an assumption that the discreteness of the levels leads to a slowing down of the evolution cascade at $`\epsilon <\epsilon _d`$. Such an assumption  seems to be quite natural in view of the fact that discrete harmonics of different frequencies are out of resonance, which should suppress the cascade in the wavenumber space. If this is the case, the next stage of evolution at the energies $`\epsilon <\epsilon _d`$ will be a fast equilibration \[with the characteristic time $`\tau _{\mathrm{coll}}(\epsilon =\epsilon _d,n_\epsilon =A/\epsilon _d^\alpha )`$; harmonics with energies $`\epsilon \epsilon _d`$ equilibrate due to the interaction with harmonics $`\epsilon \epsilon _d`$, rather than with each other\]. Further evolution should lead to the destruction of the non-equilibrium distribution $`n_\epsilon =A/\epsilon ^\alpha `$ at $`\epsilon \epsilon _d`$, the process having the form of a self-similar back cascade in the wavenumber space (back wave) :

$$n_\epsilon (t)=A\epsilon _0^\alpha (t)\stackrel{~}{f}(\epsilon /\epsilon _0(t)),\epsilon >0,t>t_{}$$
(12)
\[$`\stackrel{~}{f}(x)f(x)`$ at $`x\mathrm{}`$\], where $`\epsilon _0(t)`$ is given by the same formula (3). Clearly, Eq. (12) is applicable only if $`tt_{}`$ greater than the above-mentioned time of equilibration of the harmonics $`\epsilon <\epsilon _d`$. As follows from an estimate of corresponding terms of collision integral , the back wave should create a quasi-equilibrium distribution at $`\epsilon \epsilon _0(t)`$, that is $`\stackrel{~}{f}(x)1/x`$ at $`x1`$, with infinite at $`t=t_{}`$ and ever decreasing afterwards temperature $`\epsilon _0^\alpha (t)`$.

Since after the formation of the back wave all the harmonics with $`\epsilon \epsilon _0(t)`$, including the lowest one (!), are in a quasi-equilibrium, the process of long-range ordering in the case (i) is somewhat trivial: The state of the harmonics $`\epsilon \epsilon _0(t)`$ is determined by the instant temperature. In particular, the number of condensate particles $`N_0(t)`$ follows from Eq. (12) by conservation of the total number of particles:

$$N_0(t)=\frac{2A}{\omega _x\omega _y\omega _z}\epsilon _0^{3\alpha }(t)_0^{\mathrm{}}𝑑xx^2[x^\alpha \stackrel{~}{f}(x)].$$
(13)
This universal law of the condensate build-up should take place until the back wave cascade reaches the initial energy scale and the further growth of $`N_0`$ becomes sensitive to the details of preparation of the initial state.

We now turn to the most non-trivial case (ii), where the condition $`\xi 1`$ occurs when the front of the wave $`\epsilon _0(t)`$ reaches some energy $`\epsilon _\xi `$. From (1)-(3) and (6) we have the following estimate for the parameter $`\xi `$ at the front of the wave:

$$\xi _{\epsilon _0}\omega _0\epsilon _0^{2(\alpha 1)}/ma^2A^2.$$
(14)
Substituting $`A`$ from (4) into (14) and setting $`\xi _{\epsilon _0}1`$, we find

$$\epsilon _\xi \omega _0p^{1/2(\alpha 1)},$$
(15)
which, in particular, means that the case (ii) can only occur if $`p>1`$. From Eqs. (10) and (15) we conclude that, apart a cross-over region (corresponding to $`p1`$), the following alternative takes place. Either discreteness of the low-lying levels becomes relevant ($`p1`$), or there takes place a break-down of Knudsen regime ($`p1`$). Geometric size $`R_{}`$ of the area where the break-down of Knudsen regime takes place is obviously related to $`\epsilon _\xi `$ by $`m\omega _0^2R_{}^2\epsilon _\xi `$, which gives $`R_{}l_0p^{1/4(\alpha 1)}`$, where $`l_0=1/\sqrt{m\omega _0}`$ is the typical size of the single-particle groundstate wavefunction in the trap.

A natural question then is: What happens at distances $`r<R_{}`$ from the center of the trap after the Knudsen regime is broken down? A self-consistent scenario of future evolution can be constructed by noticing that further increase of the collisional rate due to the stimulated scattering naturally leads to the anti-Knudsen regime, when the free-path-length of the low-energy particles becomes much smaller than $`R_{}`$. This immediately leads to a quasi-homogeneous picture of evolution: the region $`r<R_{}`$ can be formally thought of as a number of independent homogeneous sells of the size much larger than the free-path-length, but much smaller than $`R_{}`$. Evolution in each sell follows homogeneous scenario, ultimately leading to the formation of the quasicondensate. Indeed, in the anti-Knudsen regime the evolution during the time period on the order of collision time is insensitive to the external potential (the criteria for the Knudsen regime and for the sensitivity to the potential within the collision time coincide). But this time is enough to form the wave (2)-(3) in the wavenumber space (with the exponent $`\alpha `$ now corresponding to the homogeneous case) and then to form quasicondensate. During the evolution in the wavenumber space, the free-path length of the particles with $`\epsilon \epsilon _0(t)`$ is getting progressively smaller, which renders the proposed scenario self-consistent.

A minor deviation from the pure homogeneous picture is that now the moment $`t_{}`$ depends on the distance from the center of the potential, so that the coherent regime first should start at $`r=0`$ (the point of maximal initial density) and then gradually occupy all the anti-Knudsen region up to $`rR_{}`$. Thus, the process of quasicondensation takes on a form of wave propagating in the coordinate space from $`r=0`$ towards $`rR_{}`$. When this wave reaches $`rR_{}`$, the process of formation of the quasicondensate correlation properties is finished, and within the region $`r<R_{}`$ we have a state of superfluid turbulence.

Note that the fact that the quasicondensate formation in the anti-Knudsen regime goes independently of the external potential should lead to certain qualitative effects. For example, after the anti-Knudsen regime is well formed, the trapping potential can, in principle, be turned off without any dramatic effect on the (short-range) ordering process. Clearly, experimental realization of such a regime would be interesting and instructive from the fundamental point of view. Another effect, coming from the same reason, is that the quasicondensate, just upon its formation, is out of equilibrium with respect to a global coherent motion in the potential, which means that a breathing mode for the quasicondensate motion should ultimately be excited in a trapped gas.

As far as the evolution in the anti-Knudsen region is of quasi-homogeneous character, all the quantities characterizing this process can be estimated from homogeneous relations of Refs. . As an initial condition to the corresponding homogeneous case one should take an extrapolation of the inhomogeneous solution (2)-(3) to the region of breaking down the Knudsen regime. Hence, the typical initial homogeneous single-particle energy and occupation number are nothing else than $`\epsilon _\xi `$ and $`A/\epsilon _\xi ^\alpha `$, respectively. Correspondingly, the number of particles involved in the quasi-homogeneous BEC process can be estimated as

$$\stackrel{~}{N}\gamma ^{1\alpha /3}N^{\alpha /3}p^{(3\alpha )/2(\alpha 1)}.$$
(16)
An important step of the qusi-homogeneous regime is associated with coming the homogeneous back wave (2)-(3) to the initial energy scale (which in our case corresponds to $`\epsilon _\xi `$). At this stage the number of qussicondensate particles $`\stackrel{~}{N}_0`$ is on the order of the total number of particles. Eq.(16) thus yields an estimate of the number of quasicondensate particles which will be created in the course of anti-Knudsen self-evolution in the system, even if the trapping potential is turned off upon the formation of the anti-Knudsen regime.

The state of superfluid turbulence, which arises together with the formation of the quasi-condensate local correlation properties, is characterized by the typical separation between the vortex lines. This separation $`r_0`$ is minimal just upon the quasicondensate formation, when $`r_0`$ is on the order of the quasicondensate healing length. From the estimate of the quasicondensate density upon its formation  we find $`r_0`$ and come to the following relation:

$$R_{}/r_0p^\nu ,\nu =\frac{4\alpha _03}{4(\alpha 1)(2\alpha _01)},$$
(17)
where $`\alpha _01.24`$ is the exponent $`\alpha `$ for the homogeneous case. From (17) we see once again that the formation of the superfluid turbulence is possible only if $`p1`$.

In conclusion, we note that in the experiment of Ref.  the parameter $`p`$ is, roughly speaking, of order unity. This means that the full-scale non-equilibrium BEC scenario cannot be seen with these initial conditions. Theoretically, the most direct way of achieving $`p1`$ is to radically increase $`\omega _0`$ immediately after the formation of the non-equilibrium initial state, during the time which is much less than the collisional time.


text:
# Lanczos study of the 𝑆=1/2 frustrated square-lattice antiferromagnet in a magnetic field

## I Introduction

Frustrated magnets are known to exhibit rich physics and diverse properties in zero external field (compare this volume and ). Exciting behavior in an external magnetic field was also observed recently in several cases: For instance, the $`S=7/2`$ three-dimensional antiferromagnet Gd3Ga5O12 is disordered (but gapless) in zero external field and orders only in finite magnetic fields . Another example is the recently discovered two-dimensional $`S=1/2`$ spin-gap material SrCu2(BO3)2 which exhibits several plateaux in the magnetization curve . The origin of these plateaux and the nature of the corresponding spin states are under intense debate  and it is suggested that some of the plateaux in SrCu2(BO3)2 might be described by certain types of ordered states.

A plateau in the magnetization curve of the triangular lattice at one third of the saturation magnetization is theoretically well known (see, e.g., ). However, general magnetization plateaus in two dimensions (including those of SrCu2(BO3)2) are still not well understood despite recent progress on the subject .

## II The model

In  we have studied the aforementioned issues in one of the simplest two-dimensional frustrated magnets, the frustrated square lattice antiferromagnet whose Hamiltonian is given by

$$H=J\underset{\stackrel{}{x},\stackrel{}{y}}{}\stackrel{}{S}_\stackrel{}{x}\stackrel{}{S}_\stackrel{}{y}+J^{}\underset{[\stackrel{}{x},\stackrel{}{y}]}{}\stackrel{}{S}_\stackrel{}{x}\stackrel{}{S}_\stackrel{}{y}h\underset{\stackrel{}{x}}{}S_\stackrel{}{x}^z$$
(1)
where the first sum is over nearest neighbor pairs $`\stackrel{}{x},\stackrel{}{y}`$ and the second one over diagonal neighbor pairs $`[\stackrel{}{x},\stackrel{}{y}]`$. Here we will present details of  and further results obtained from a numerical diagonalization of the $`S=1/2`$ version of the Hamiltonian (1) with the Lanczos method.

At zero field, the frustrated square lattice $`S=1/2`$ Heisenberg model has been studied extensively using also exact diagonalization . One reason for the popularity of this model is that it exhibits melting of long-range magnetic order by quantum fluctuations: For small frustration $`J^{}J`$, one finds antiferromagnetic Néel order while for $`J^{}J`$ a certain type of collinear order described by a wave vector $`(\pi ,0)`$ or $`(0,\pi )`$. Finally, in an intermediate region around $`J^{}=J/2`$, a spin-liquid groundstate with a gap to magnetic excitations is found.

In the presence of an external field, an important quantity is the magnetization

$$M=\frac{1}{SV}\underset{\stackrel{}{x}}{}S_\stackrel{}{x}^z$$
(2)
which we normalize to saturation values $`|M|=1`$. Magnetization curves for the model (1) have already been computed by exact diagonalization . However, the study of  was restricted to a special lattice of $`4\times 4`$ sites which can also be interpreted in terms of several other geometries. For instance, a strip of width 4 is equivalent to a frustrated four-leg ladder and then plateaux are expected in the magnetization curve at $`M=0`$, $`1/4`$, $`1/2`$ and $`3/4`$ (see, e.g., ). Precisely these values were also observed in , but it remains to be clarified whether they are an artifact of the special geometry. The other study  excluded the region $`J^{}>J/2`$ which we find to be the most interesting one and boundary conditions were used in  which frustrate the order which wants to develop.

## III Single magnon dispersion

It is instructive to look first at the dispersion of a single magnon above the ferromagnetic state with $`M=1`$ (all spins aligned along the field). When we formally set $`h=0`$, the one-magnon dispersion for the $`S=1/2`$ model is given by

$`_{1\mathrm{s}}(k_x,k_y)`$ $`=`$ $`2(J+J^{})+J\left(\mathrm{cos}k_x+\mathrm{cos}k_y\right)`$ (4)
$`+2J^{}\mathrm{cos}k_x\mathrm{cos}k_y.`$
For $`J^{}<J/2`$, this dispersion has a single minimum at $`k_x=k_y=\pi `$ corresponding to Néel order. On the other hand, for $`J^{}>J/2`$ two equivalent minima are found at $`k_x=0`$, $`k_y=\pi `$ or $`k_x=\pi `$, $`k_y=0`$ which signal collinear order. The case $`J^{}=J/2`$ is special: Here lines of minima are found for either $`k_x=\pi `$ or $`k_y=\pi `$.

## IV Magnetization curves

We now discuss the behavior in the entire field range using exact diagonalization. Guided by the dispersion (4), we have considered only rectangular clusters with an even number of sites in both directions. Finite-size magnetization curves are shown in Fig. 1 for $`J^{}/J=0.6`$. As we have pointed out in section II, finite-size effects must be examined carefully. Even if the $`M=3/4`$ and $`1/4`$ plateaux are quite pronounced on the $`4\times 4`$ and $`6\times 4`$ clusters, we believe that they are artifacts of the special system sizes. Indeed, the data for wider strips indicates a vanishing plateau at $`M=3/4`$ in the thermodynamic limit. In contrast, the plateau at $`M=1/2`$ is only a little narrower on the $`6\times 6`$ cluster than for the smaller clusters. In the extrapolation (bold line in Fig. 1), we have therefore drawn a plateau with $`M=1/2`$ but none at $`M=1/4`$ and $`3/4`$. The value of the spin gap (i.e. the boundary of the $`M=0`$ plateau) has been taken from . Otherwise the extrapolation was obtained by the standard procedure of connecting the mid-points of the steps of the magnetization curves at the largest available system size.

Fig. 2 shows the projection of the magnetization curves onto the horizontal axis, i.e. the locations of the finite-size jumps in the magnetization curves. From this one can then first read off how the width of the $`M=1/2`$ plateau varies with $`J^{}/J`$. An analysis of its finite-size behavior  indicates that the $`M=1/2`$ plateau exists in the region $`0.5J^{}/J0.65`$.

Fig. 2 further shows that the transition to saturation is special at $`J^{}=J/2`$, as is expected from (4): A jump of size $`\delta M=1/L`$ occurs just below saturation where $`L`$ is the length of the shorter edge of the cluster. This implies a smooth transition to saturation in the thermodynamic limit. The particularity of the transition to saturation for $`J^{}J/2`$ was already noted in  where, however, the jump was obscured by the choice of cluster geometries. We therefore believe the asymptotic form proposed in  just to be a manifestation of these boundary effects in combination with the questionable assumption that $`h(M)`$ has a power series expansion. It should be possible to obtain the correct asymptotic form from Bose condensation of magnons  into the lines of minima of the single magnon dispersion (4), but this computation has not yet been carried out.

For $`J^{}J/2`$, the magnetization curves further show a pronounced finite-size plateau at $`M=11/L`$, i.e. just below the jump. It was actually proposed in  that a plateau just below saturation should exist for $`J^{}J/2`$, but our candidate disappears as $`L\mathrm{}`$ and we suspect an error in the treatment  of the single-magnon dispersion (4).

## V Static structure factors

To gain further insight into the phase diagram and in particular in order to understand the nature of the $`M=1/2`$ plateau state better, we use the static structure factors which are defined as*** The symmetrized variant $`\left(S^{\alpha \beta }(k_x,k_y)+S^{\alpha \beta }(k_y,k_x)\right)/2`$ is used when the cluster is not square or the groundstate not symmetric under reflection across the diagonal of the cluster ($`xy`$).

$$S^{\alpha \beta }(k_x,k_y)=\frac{1}{V^2}\underset{x,y,r,s}{}\mathrm{e}^{i\left(k_xx+k_yy\right)}S_{r,s}^\alpha S_{r+x,s+y}^\beta .$$
(5)
The finite magnetization leads to a trivial peak in $`S^{zz}(0,0)`$. We observe additional peaks in $`S^{\alpha \beta }(k_x,k_y)`$ at $`(0,\pi )`$ or $`(\pi ,0)`$ and $`(\pi ,\pi )`$. The evolution of these peaks with $`J^{}/J`$ is shown in Figs. 3 and 4 for two values of $`M`$.

Néel and collinear order develops in the transverse components for small and large $`J^{}`$, respectively. Just at $`M=1/2`$ and for $`0.5J^{}/J0.65`$, the longitudinal components develop uuud order with three spins pointing up and one down on a four-spin plaquette. While it may be possible that some order exists in $`S^{zz}(0,\pi )`$ at $`M=2/3`$ for intermediate $`J^{}`$ (see Fig. 3), the structure factors at other magnetizations and Monte Carlo simulations of the classical model indicate disordered states for magnetizations sufficiently far above or below the uuud state at $`M=1/2`$ .

## VI Groundstate quantum numbers

For an odd number of magnons $`V/2S^z`$, one expects the groundstate to carry $`k_x=k_y=\pi `$ (corresponding to Néel order) for $`J^{}J`$ while for $`J^{}J`$ one expects collinear order and thus $`k_x=0`$, $`k_y=\pi `$ or $`k_x=\pi `$, $`k_y=0`$ – compare the discussion of the one-magnon sector in section III. When $`V/2S^z`$ is even, the groundstate carries momentum $`k_x=k_y=0`$. However, as was pointed out in , groundstate level crossings are still expected in the sectors with $`V/2S^z=2n`$ and $`n`$ odd. Then the groundstate is even under diagonal reflections for $`J^{}J`$ and odd for $`J^{}J`$.

The unique groundstate level crossings are shown in Fig. 5 (in some sectors more than one groundstate level crossing is observed). The crossing at $`M=0`$ on the $`6\times 6`$ cluster was taken from . Note that according to (4), the crossing point must move to $`J^{}=J/2`$ for $`M1`$ where a first-order transition between Néel order and collinear order occurs. The position of the crossing at $`M=1/2`$ agrees with our earlier estimates for the transition between the uuud plateau state and collinear order. Also the crossing at $`M=0`$ is close to recent estimates (see, e.g., ) for the location of the zero-field transition between the disordered spin-liquid state and the collinear state.

## VII Summary

We have studied the frustrated square lattice $`S=1/2`$ Heisenberg model in the presence of a strong magnetic field using the Lanczos method. A plateau with $`M=1/2`$ exists in the region $`0.5J^{}/J0.65`$ and corresponds to an up-up-up-down ordered state. The transverse spin components show Néel and collinear order for small and large $`J^{}`$, respectively. The intermediate region contains the ordered $`M=1/2`$ state and disordered regions at higher and lower fields . Finally, we have found that the transition to saturation is special at $`J^{}=J/2`$.

Acknowledgments: I am indebted to O.A. Petrenko and M.E. Zhitomirsky for collaboration on the subject , to S. Gluzman and D. Poilblanc for discussions and suggestions as well as to the Deutsche Forschungsgemeinschaft for funding participation in the conference ‘Highly Frustrated Magnetism 2000’.