text:
# Spin squeezing via quantum feedback

## Abstract

We propose a quantum feedback scheme for producing deterministically reproducible spin squeezing. The results of a continuous nondemolition atom number measurement are fed back to control the quantum state of the sample. For large samples and strong cavity coupling, the squeezing parameter minimum scales inversely with atom number, approaching the Heisenberg limit. Furthermore, ceasing the measurement and feedback when this minimum has been reached will leave the sample in the maximally squeezed spin state.

Squeezed spin systems KitUed93  of atoms and ions have attracted considerable attention in recent years due to the potential for practical applications, such as in the fields of quantum information QInf98  and high-precision spectroscopy Spect . Quantum correlations of squeezed spin states outperform classical states in analogy with squeezed optical fields. Moreover, squeezed spin sates are also multiparticle entangled states Soretal01 . Recent proposals for their generation include the absorption of squeezed light sqzlight , collisional interactions in Bose-Einstein condensates Soretal01 ; collBEC , and direct coupling to an entangled state through intermediate states such as collective motional modes for ions MolSor99  or molecular states for atoms HelmYou01 .

Other proposals create spin squeezing via quantum nondemolition (QND) measurements KuzBigMan98 ; KuzManBig00 ; Duanetal00 . A striking recent achievement of QND measurements is the entanglement of two macroscopic atomic samples JulKozPol01 . These QND schemes produce conditional squeezed states that are dependent on the measurement record. On the other hand, unconditional squeezing would ensure that the state is deterministically reproducible. Mølmer Mol99  has shown that alternating QND measurements and incoherent feedback can produce sub-Poissonian number correlations. However, that work does not treat the quantum effects of the measurement back action or the feedback on the mean spin (which is assumed to be zero). Hence it cannot predict the strength of the entanglement.

In this Rapid Communication, we suggest achieving spin squeezing via feedback that is coherent and continuous. We consider a continuous QND measurement of the total population difference of an atomic sample. The results of the measurement, which conditionally squeeze the atomic sample, are used to drive the system into the desired, deterministic, squeezed spin state. This involves amplitude modulation of a radio-frequency (rf) magnetic field, where the feedback strength varies in time such that the mean number difference is kept at zero.

An ensemble of $`N`$ two-level atoms can be described by a spin-$`J`$ system Dic54 , i.e., a collection of $`2J=N`$ spin-$`\frac{1}{2}`$ particles. The collective spin operators are given by $`J_\alpha =_{k=1}^N\frac{1}{2}\sigma _\alpha ^{(k)}(\alpha =x,y,z)`$, where $`\sigma _\alpha ^{(k)}`$ are the Pauli operators for each particle. Thus, $`J_z`$ represents half the total population difference. Coherent spin states (CSS) have variances normal to the mean spin direction equal to the standard quantum limit (SQL) of $`J/2`$. Introducing quantum correlations among the atoms reduces the variance below the SQL in one direction at the expense of the other KitUed93 . Such squeezed spin states can be characterized by the squeezing parameter Soretal01

$$\xi _{𝐧_1}^2=N(\mathrm{\Delta }J_{𝐧_1})^2/(J_{𝐧_2}^2+J_{𝐧_3}^2),$$
(1)
where $`𝐧_i(i=1,2,3)`$ are orthogonal unit vectors. Systems with $`\xi _𝐧^2<1`$ are spin squeezed in the direction $`𝐧`$ and also have multiparticle entanglement Soretal01 .

Let the internal states, $`|1`$ and $`|2`$, of each atom be the degenerate magnetic sublevels of a $`J=\frac{1}{2}`$ state, e.g., an alkali ground state. Each atom is prepared in an equal superposition of the two internal states, thus giving a CSS of length $`J`$ in the $`x`$-direction. The atomic sample is placed in a strongly driven, heavily damped, optical cavity, as shown in Fig. 1$`(a)`$. The cavity field is assumed to be far off resonance with respect to transitions probing state $`|2`$, see Fig. 1$`(b)`$. This dispersive interaction causes a phase shift of the cavity field proportional to the number of atoms in $`|2`$. Thus, the QND measurement of $`J_z`$ (since $`N`$ is conserved) is effected by the homodyne detection of the light exiting the cavity CorMil98 .

This interaction is defined by the Hamiltonian $`\mathrm{}\chi J_zb^{}b`$ where $`b`$,$`b^{}`$ are the cavity field operators and $`\chi =g^2/4\mathrm{\Delta }`$, with one-photon Rabi frequency $`g`$, and optical detuning $`\mathrm{\Delta }`$ CorMil98 . For strong coherent driving we can use the semiclassical approximation $`bi|\beta |+b`$, where now $`b`$ represents small quantum fluctuations around the classical amplitude $`i|\beta |`$. The interaction is thus

$$H_{\mathrm{int}}=\mathrm{}\chi |\beta |J_z(ib+ib^{}),$$
(2)
where we have chosen an initial $`J_z`$ splitting of $`\chi |\beta |^2`$.

Following the procedure of Sec. VII in Ref. WisMil94  we can adiabatically eliminate the cavity dynamics if the cavity decay rate $`\kappa \chi |\beta |J_z`$, which requires $`\kappa \chi |\beta |\sqrt{N}`$ (since the initial $`\mathrm{\Delta }J_z=\sqrt{J/2}`$). The evolution of the atomic system due the measurement is thus

$$\dot{\rho }=M𝒟[J_z]\rho ,$$
(3)
where $`M=4\chi ^2|\beta |^2/\kappa `$ is the measurement strength \[equivalent to $`2D`$ in Eq. (22) of Ref. CorMil98 \], and $`𝒟[r]\rho r\rho r^{}(r^{}r\rho +\rho r^{}r)/2`$. This equation represents decoherence of the atomic system due to photon number fluctuations in the cavity field, with the result of increased noise in the spin components normal to $`J_z`$.

The effect of Eq. (2) on the cavity field is a phase shift proportional to $`J_z`$, and thus the output Homodyne photocurrent is given by WisMil94

$$I_c(t)=2\sqrt{M}J_z_c+\zeta (t),$$
(4)
where $`\zeta (t)`$ is a white-noise term satisfying $`\mathrm{E}[\zeta (t)\zeta (t^{})]=\delta (tt^{})`$ and $`\mathrm{E}`$ is the ensemble average. The conditional master equation for the atomic system is then WisMil94

$$d\rho _c=dtM𝒟[J_z]\rho _c+\sqrt{M}dW(t)[J_z]\rho _c,$$
(5)
where $`dW(t)=\zeta (t)dt`$ is an infinitesimal Wiener increment and $`[r]\rho r\rho +\rho r^{}\mathrm{Tr}[(r+r^{})\rho ]\rho `$.

The effects of this evolution on the initial CSS are a decrease in the variance of $`J_z`$ with corresponding increases for $`J_y`$ and $`J_x`$ (i.e., spin squeezing), as well as a stochastic shift of the mean $`J_z`$ away from its initial value of zero. This shift, indicated by state 2 of Fig. 2, is equal to

$$dJ_z_c=2\sqrt{M}dW(t)(\mathrm{\Delta }J_z)_c^22\sqrt{M}I_c(t)dtJ_z^2_c.$$

Here the approximation assumes that $`J_z_c=0`$, which will be relevant when feedback is in place.

The average or unconditioned evolution, Eq. (3), is simply recovered by taking the ensemble average of all possible conditioned states, i.e., $`\rho (t)=\mathrm{E}[\rho _c(t)]`$. This leads to a spin state with $`J_z`$ variance equal to $`J/2`$. In other words, the unmonitored measurement does not affect $`J_z`$ and the squeezed character of the individual conditioned states is lost, indicated by state 3 in Fig. 2.

To retain the reduced fluctuations of $`J_z`$ in the average evolution, we employ a feedback mechanism that uses the measurement record to continuously drive the system into the $`J_z_c=0`$ squeezed state. The idea is to cancel the stochastic shift of $`J_z_c`$ due to the measurement. This simply requires a rotation of the mean spin about the $`y`$ axis equal and opposite to that caused by Eq. (Spin squeezing via quantum feedback), as illustrated in Fig. 2. The feedback Hamiltonian must therefore take the form

$$H_{\mathrm{fb}}(t)=\lambda (t)I_c(t)J_y/\sqrt{M},$$
(6)
where $`\lambda (t)`$ is a time-varying feedback strength. This feedback driving can be implemented by modulating an applied rf magnetic field Sangetal01 , as shown in Fig. 1.

Following again the methods of Ref. WisMil94  to find the total stochastic master equation, we can calculate the conditioned shift of the mean $`J_z`$ due to the feedback. Again using the assumption that $`J_z_c=0`$ we have

$$dJ_z_{\mathrm{fb}}\lambda (t)I_c(t)dtJ_x_c/\sqrt{M}.$$
(7)
Since the idea is to produce $`J_z_c=0`$ via the feedback, the approximations above and in Eq. (Spin squeezing via quantum feedback) apply and we find that the required feedback strength is

$$\lambda (t)=2MJ_z^2_c/J_x_c.$$
(8)
This type of feedback control is essentially a form of state-estimation-based feedback Dohetal00 . Although Eq. (6) looks like direct current feedback, the strength of this feedback (8) is determined by conditioned state expectation values. $`I_c(t)`$ only appears directly in $`H_{\mathrm{fb}}`$ due to the assumption that the feedback works and so $`J_z_c=0`$.

Being dependent on conditioned expectation values, which are computationally very expensive, Eq. (8) is not practical in an experimental sense. What is required is a predetermined series of data points or ideally an equation for $`\lambda (t)`$, like in Fig. 1. To find a suitable expression we begin by assuming the feedback is successful and replace the conditioned averages by ensemble averages,

$$\lambda (t)2MJ_z^2/J_x.$$
(9)
This approximation will be valid if the unconditioned state has high purity since then it must comprise of nearly identical highly pure conditioned states.

The evaluation of both the purity ($`\mathrm{Tr}[\rho ^2]`$) and the averages in Eq. (9) ($`A\mathrm{Tr}[A\rho ]`$) requires the unconditioned master equation (ME) WisMil94

$$\dot{\rho }=M𝒟[J_z]\rho i\lambda (t)[J_y,J_z\rho +\rho J_z]+\frac{\lambda (t)^2}{M}𝒟[J_y]\rho .$$
(10)
The terms in this equation describe, respectively, the noise due to the measurement back-action, the feedback optical driving, and the noise introduced by the feedback. The state determined by Eq. (10), with $`\lambda (t)`$ given by Eq. (9), has a purity very close to one (see below). Since the state is very close to a pure state, we are justified in applying the approximation of Eq. (9).

Note that Eq. (10) describes the exact unconditioned evolution of the atomic system where the feedback strength is arbitrarily defined by $`\lambda (t)`$. Equation (9) thus describes one particular feedback scheme, however, it can only (easily) be evaluated numerically. To find an approximate analytical expression we look at when $`𝐉𝒪(J)`$, for which the atomic sample remains near the minimum uncertainty state $`J_y^2J_z^2J^2/4`$. This is equivalent to a linear approximation represented by replacing $`J_x`$ with $`J`$ in the commutator $`[J_y,J_z]=iJ_x`$, which allows us to calculate $`J_y^2`$ directly from the ME (and hence $`J_z^2`$). The decrease of $`J_x`$ from $`J`$ is then related to the increase of $`J_y^2`$ from $`J/2`$ \[see Fig. 2$`(b)`$\] due to the measurement back action. Using these approximations we obtain

$$\lambda (t)Me^{Mt/2}(1+2JMt)^1.$$
(11)
We can analytically approximate the degree of squeezing produced by the particular feedback scheme represented by Eq. (11). For our model Eq. (1) becomes

$$\xi _z^2=\frac{NJ_z^2}{J_x^2}\frac{\lambda (t)J}{MJ_x}e^{Mt}(1+2JMt)^1.$$
(12)
This leads to a minimum at $`t_{}1/M`$ of

$$\xi _{\mathrm{min}}^2e/2J,J1.$$
(13)
Thus, the minimum attainable squeezing parameter asymptotically approaches an inverse dependence on the sample size, i.e., the Heisenberg limit AndLuk01 .

The approximations leading up to Eqs. (9) and (11) can be justified by numerically solving the ME (10) for a given $`\lambda (t)`$, for example, using the Matlab quantum optics toolbox Tan99 . The approximation of Eq. (9) can be tested by calculating the purity for $`\rho `$ described by the ME with this particular feedback. The expectation values in $`\lambda (t)`$ are found by iteratively solving the ME \[updating $`\lambda (t)`$ each time step\], and thus we also have

$$\xi _z^2=2JJ_z^2/J_x^2.$$
(14)
The results of this simulation are shown in Fig. 3$`(a)`$, where the purity is given by the dotted curve and $`\xi _z^2`$ is curve A. Clearly, the purity remains near unity for times of interest. This implies that the measurement and feedback scheme has worked to produce nearly identical, nearly pure, conditioned states \[for times $`𝒪(M^1)`$\], and we are therefore justified in using Eq. (9).

The further approximations to obtain the analytical Eq. (11) are also good, as shown by curve B in Fig. 3$`(a)`$, where the fit to curve A for the early evolution and minimum is nearly perfect. We are not interested in later times since the idea is to cease the measurement and feedback when the minimum $`\xi _z^2`$ is reached. The analytical expression for the squeezing parameter, Eq. (12), is also plotted as curve C in Fig. 3$`(a)`$. Although the minimum is not a perfect fit to the exact numerical results, it has the correct order of magnitude and so we expect the $`J^1`$ scaling of Eq. (13) to be correct.

The scaling of $`\xi _{\mathrm{min}}^2`$ is obtained numerically from solutions of the ME \[and thus Eq. (14)\] with feedback described by Eq. (11). This $`\lambda (t)`$, shown above to be a good approximation, is the suitable form for experimental realization. The numerical results, along with the analytical expression (13), are plotted in Fig. 3$`(b)`$ which clearly verifies the $`J^1`$ dependence. The analytical coefficient ($`e/2`$) represents an error of $`18\%`$ compared to the numerical fit, and the optimum time ($`t_{}`$) is also slightly out as shown by Fig. 3$`(a)`$. Nevertheless, these errors only apply to scaling coefficients, not to the scalings themselves.

Experimentally, the limit to squeezing will be dominated by spontaneous losses due to absorption of QND probe light. The rate of this loss is $`N\gamma \mathrm{\Omega }^2/4\mathrm{\Delta }^2`$, where $`\gamma `$ is the spontaneous emission rate and $`\mathrm{\Omega }=g|\beta |`$. To reach the Heisenberg limit (requiring a time $`t_{}=1/M`$) we want the total loss to be negligible, i.e., we require

$$N\gamma \frac{\chi |\beta |^2}{\mathrm{\Delta }}\frac{\kappa }{4\chi ^2|\beta |^2}1g^2N\kappa \gamma ,$$
(15)
which is the very strong coupling regime of cavity QED. Note we also required $`\kappa \chi |\beta |\sqrt{N}`$ for the adiabatic elimination of the cavity field, and to satisfy both we thus require $`\mathrm{\Delta }\gamma N^{3/2}`$. It is not surprising that Eq. (15) is the same requirement as for André and Lukin’s model AndLuk01  implementing the countertwisting Hamiltonian KitUed93 , since writing Eq. (10) in Lindblad form reveals such a term. Similarly, the condition for achieving some squeezing, i.e., $`\xi ^2<1`$, will be $`g^2N>\kappa \gamma `$. Further, we have calculated that a free space model \[also given by Eq. (10) but with $`M`$ equal to $`N`$ in Ref. WisTho01 \] will also produce some squeezing, although the Heisenberg limit cannot be reached since by $`t_{}=1/M`$ all atoms will be lost from the sample.

Figure 3$`(b)`$ also indicates that our continuous scheme is very robust to any experimental errors in the feedback strength, as opposed to a single-shot method. The latter approach consists of a single (integrated) measurement pulse (see e.g., Ref. Duanetal00 ), followed by a single feedback pulse. If there is a relative error of $`ϵ`$ in the feedback strength, this will induce an error term $`(\mathrm{\Delta }J_z^{\mathrm{err}})^2ϵ^2J/2`$, which will dominate the total variance for $`J1`$. Thus $`\xi _{\mathrm{min}}^2`$ will have a lower bound of $`ϵ^2`$, and will never be better than $`J^1`$. On the other hand, as shown in Fig. 3$`(b)`$, a large ($`20\%`$) error in $`\lambda (t)`$ for continuous feedback does not affect the $`J^1`$ scaling. We have also found this theoretical scaling to be unaffected by inefficient measurements. Finite feedback delay time will also have a limited effect as long as it is faster than $`(JM)^1`$.

This Rapid Communication has presented a scheme for producing a spin squeezed atomic sample via QND measurement and feedback. The advantage over previous QND schemes KuzBigMan98 ; KuzManBig00  is that it provides unconditional, or deterministically reproducible, squeezing. For very strong cavity coupling, the theoretical squeezing approaches the Heisenberg limit $`\xi ^21/N`$, while some squeezing will be produced at weaker coupling and even in free space (thus presenting a simple experimental test for quantum feedback). This indicates a stronger squeezing mechanism than collisional interactions in a Bose-Einstein condensate where the scaling is $`N^{2/3}`$ Soretal01 ; collBEC . Furthermore, by ceasing the measurement when this minimum is reached, the maximally squeezed state could be maintained indefinitely.

L.K.T and H.M.W wish to acknowledge inspiring discussions with Klaus Mølmer and Eugene Polzik.


format: text

text:
# Do Distinct Cosmological Models Predict Degenerate Halo Populations?

## 1 Introduction

In cosmological models dominated by cold dark matter (CDM) and an unclustered energy component (such as a cosmological constant), gravitational instability of primordial density fluctuations produces a population of dark matter halos, each in approximate virial equilibrium. Depending on their mass, these halos may host individual galaxies, galaxy groups, or rich galaxy clusters. Provided that one focuses on systems of overdensity $`\rho /\overline{\rho }200`$, and regards higher density structures as subsystems within their parent halos, the halo population itself is insensitive to the gas pressure forces that influence the sub-dominant baryon component. In this paper, we investigate whether changes to cosmological parameters — specifically the matter density parameter $`\mathrm{\Omega }_m`$ and parameters that describe the shape and amplitude of the primordial power spectrum — always produce measurable changes in the dark halo population, or whether two models with different combinations of these parameters can give rise to halo populations that are effectively indistinguishable.

Our interest in this question is spurred by recent developments in the theory of biased galaxy formation. The uncertain relation between galaxies and dark matter is the primary limitation in testing cosmological models against observations of galaxy clustering. The “halo occupation distribution” (HOD) characterizes this relation statistically in terms of the probability distribution $`P(N|M)`$ that a halo of virial mass $`M`$ contains $`N`$ galaxies of a specified type, together with prescriptions that specify the relative spatial and velocity distributions of galaxies and dark matter within these halos. Numerous recent papers have shown that the HOD framework is a powerful tool for analytic and numerical calculations of clustering statistics, for modeling observed clustering, and for characterizing the results of semi-analytic or numerical studies of galaxy formation (e.g., Kauffmann, Nusser, & Steinmetz 1997; Jing, Mo, & Börner 1998; Kauffman et al. 1999; Benson et al. 2000; Ma & Fry 2000; Peacock & Smith 2000; Seljak 2000; Berlind & Weinberg 2002; Bullock et al. 2001; Marinoni & Hudson 2002; Scoccimarro et al. 2001; Yoshikawa et al. 2001; White, Hernquist, & Springel 2001). In particular, Peacock & Smith (2000), Marinoni & Hudson (2002), and Berlind & Weinberg (2002) have argued that the HOD can be determined empirically from observed galaxy clustering, given an assumed cosmological model that determines the mass function and spatial and velocity clustering of the dark halo population. Empirical determinations of the HOD can provide insight into the physics of galaxy formation, and they may sharpen the ability of large scale structure studies to test cosmological models, since a model with an incorrect dark halo population may be unable to match the data for any choice of HOD (see Berlind & Weinberg 2002 for further discussion).

Suppose that we find a combination of cosmology and HOD that reproduces all aspects of observed galaxy clustering. Can we infer that the cosmology and the derived HOD are both correct? If there is another cosmological model that predicts an indistinguishable halo population, then the answer is clearly no, since the combination of this alternative cosmology with the same HOD would predict identical galaxy clustering. These cosmological models would be degenerate in the sense that they could not be distinguished by galaxy clustering data without relying on a predictive theory of galaxy formation (which might yield different HODs for the two cosmologies). Cosmological models that predict distinguishable halo populations could still be degenerate in this sense, if changes to the HOD can mask the differences in the halo populations; we will speculate on this point in §4, but we reserve a quantitative examination of it to future work.

We focus our analyses on CDM cosmological models with Gaussian initial density fluctuation fields. Motivated by cosmic microwave background measurements (Netterfield et al., 2002; Pryke et al., 2002), we restrict our attention to spatially flat models with a cosmological constant or $`\mathrm{\Omega }_m=1`$, though we demonstrate in passing that galaxy clustering data at $`z=0`$ probably cannot distinguish a flat model from an open model with the same $`\mathrm{\Omega }_m`$. The parameters that define our cosmological models are: the matter density parameter $`\mathrm{\Omega }_m`$; the normalization $`\sigma _8`$ of the power spectrum $`P(k)`$, which is the r.m.s. fluctuation of the linear density field filtered with a top-hat filter of radius 8$`h^1\mathrm{Mpc}`$ ; and the shape of $`P(k)`$, which we characterize by the spectral index $`n_s`$ of the inflationary power spectrum and the shape parameter $`\mathrm{\Gamma }`$ of the matter transfer function (see, e.g., Efstathiou, Bond, & White 1992). Although existing data impose constraints on these parameters, individually or in combination, here we allow each to vary independently over a fairly broad range, so that we can isolate the physical effects of the matter density, the amplitude of mass fluctuations, and the shape of the power spectrum on the resulting halo population. We concentrate entirely on the halo populations at $`z=0`$, since this is where galaxy clustering data will be good enough to allow empirical HOD determinations in the near future.

Using $`N`$-body simulations described in §2, we measure (in §3) mass functions, two-point correlation functions, mean pairwise radial velocities, and pairwise velocity dispersions of dark matter halos. First, we consider models that have the same initial power spectrum $`P(k)`$ but different matter density parameter $`\mathrm{\Omega }_m`$. We then move to cluster normalized models in which $`\sigma _8`$ is changed to compensate the change in $`\mathrm{\Omega }_m`$, so that the amplitude of the halo mass function is kept approximately fixed at a cluster scale. Finally, we investigate models in which both the amplitude and the shape of $`P(k)`$ are changed in order to match the amplitude and slope of the halo mass function at the cluster scale. We review our results and briefly discuss their implications for cosmological tests in §4. Our results overlap previous numerical studies of the halo mass function and halo clustering (recent examples include Jing 1998; Governato et al. 1999; Colberg et al. 2000; Jenkins et al. 2001), but they differ in the examination of controlled parameter sequences rather than specific cosmological models. Also, while some of these studies have focused on the large scale correlations of cluster mass halos, we devote considerable attention to the lower masses and smaller spatial scales that are important for non-linear galaxy clustering. Throughout the paper, we adopt a Hubble constant $`H_0=100h\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1=70\mathrm{km}\mathrm{s}^1\mathrm{Mpc}^1`$ (Freedman et al., 2001).

## 2 Numerical Methods and Tests

In this work, statistics are mainly based on simulations run with a Particle-Mesh (PM) $`N`$-body code (Park, 1990). We run PM simulations with different combinations of density and power spectrum parameters. Each simulation follows the evolution of 2003 particles in a periodic cube 200$`h^1\mathrm{Mpc}`$ on a side, with a 4003 mesh to compute the gravitational force starting at $`z=19`$ and advancing to $`z=0`$ in 40 equal steps of expansion factor $`a`$. The mass of each particle is $`2.78\times 10^{11}\mathrm{\Omega }_mh^1M_{}`$. We adopt the parameterization of Efstathiou, Bond, & White (1992) for the fluctuation power spectrum $`P(k)`$. We divide our simulations into three categories corresponding to parameter changes in the matter density and the normalization and the shape of the initial fluctuation power spectrum. For each category, four cosmological models are investigated with $`(\mathrm{\Omega }_m,\mathrm{\Omega }_\mathrm{\Lambda })=`$(0.2,0.8), (0.3,0.7), (0.4,0.6) and (1.0,0.0). We generally regard the (0.3,0.7) model as the central one for comparison. In the central model, the spectral index is $`n_s=1.0`$, the spectral shape parameter is $`\mathrm{\Gamma }=0.20`$, and the r.m.s. mass fluctuation within spheres of radius 8$`h^1\mathrm{Mpc}`$ is $`\sigma _8=0.9`$, which is consistent with the observed abundance of clusters (Eke, Cole, & Frenk, 1996). For each model, we generate four independent realizations, and we use the dispersion among these realizations to estimate the statistical error bars on our clustering measures associated with the finite simulation volume. The total volume is comparable to that expected for the Sloan Digital Sky Survey (York et al., 2000) out to its median redshift $`z_{\mathrm{med}}0.1`$. We use a FoF algorithm (Davis et al., 1985) to identify the dark matter halos, and, unless explicitly stated otherwise, the linking length is 0.2 times the mean inter-particle separation. With this linking length, the FoF algorithm picks out structures of typical mass overdensity $`\rho /\overline{\rho }200`$ (Davis et al., 1985).

The advantage of a PM code is that it is relatively cheap to run multiple simulations with a large dynamic range in mass. The limitation of a PM code is its moderate force resolution, but this limitation is not too serious for our purposes, since we only identify halos and do not attempt to measure their internal structure. We perform two tests to make sure that the PM resolution is adequate for our purpose. First, we make a comparison between PM simulations and simulations run with the GADGET tree code (Springel, Yoshida, & White, 2001). Second, we perform a self-similar scaling test for the PM code.

GADGET is a publicly available $`N`$-body code that calculates particle accelerations by the hierarchical tree method of Barnes & Hut (1986). Its advantages over the PM algorithm are its spatial and temporal adaptability. The gravitational softening length, $`ϵ`$, can be set significantly smaller than the initial particle grid spacing to increase force resolution, though if $`ϵ`$ is too small 2-body relaxation effects will begin to influence the halo population. Particles have individual timesteps, which can vary continuously to ensure accurate time integration given the adopted force resolution. GADGET has fully periodic boundary conditions for comoving integrations via the Ewald summation method (Hernquist, Bouchet, & Suto, 1991).

We have compared the results of two simulations with 1003 particles in a 100$`h^1\mathrm{Mpc}`$ box, one evolved with the PM code using a 2003 mesh (thus having equal resolution to the larger simulations on which our results will be based), and one evolved with GADGET using $`ϵ=0.3h^1`$Mpc. The initial conditions of the two simulations were identical to each other and used the parameters of the central model described above. Slices through these simulations are plotted for comparison in Figure 1. The structures formed in the two simulations are very similar at the level of detail discernible in this plot. The halo mass functions of the two simulations, also shown in Figure 1, trace each other nearly exactly. The PM simulation produces slightly ($`10\%`$) fewer halos below $`10^{13}h^1M_{}`$, which is to be expected given its lower force resolution; the diameter of an overdensity 200 sphere of $`10^{13}h^1M_{}`$ is 1.04$`h^1\mathrm{Mpc}`$, or about two PM grid cells. When the mesh size on this PM simulation was increased to 3003, the resulting mass function was consistent with the GADGET simulation even on these small mass scales.

In addition to the halo mass function, the quantities important to our analysis are the correlation functions for both the mass distributions and the halo populations, which are shown in Figure 2. As with the halo mass functions, the matter correlation functions for the PM and GADGET simulations match very well. They begin to deviate at separations smaller than $`0.5h^1`$Mpc, or under one mesh cell. The halo correlation functions produced by the two methods are also consistent. Larger simulations and multiple realizations of each model are efficiently produced with the PM code, which takes up to an order of magnitude less computing time at the force resolution adopted here.

We next perform a self-similar scaling test for the PM code. Assuming an Einstein-de Sitter universe ($`\mathrm{\Omega }_m=1`$), we use a pure power-law power spectrum $`P(k)k^{n_s}`$ with $`n_s=1`$ (and no transfer function modulation) for the initial density fluctuations of the scale-free model. A 4003 force mesh is used to follow the gravitational evolution of 2003 particles. The normalization of the power spectrum is set so that the nonlinear mass $`M_{}`$ at the final output corresponds to 400 particles. We define $`M_{}`$ by the condition $`\sigma (M_{})=\delta _c=1.69`$, where $`\sigma (M_{})`$ is the r.m.s. linear mass fluctuation on mass scale $`M_{}`$. The r.m.s. fluctuation is $`\sigma =1.0`$ at about 1/26 of the box size. If we identify this scale as 8$`h^1\mathrm{Mpc}`$ , the implied size of the periodic simulation cube is then 207.62$`h^1\mathrm{Mpc}`$ at $`z=0`$. Earlier outputs can be identified with higher redshifts or with larger volumes at $`z=0`$. In the language of the redshift description, the simulation is evolved for 40 time steps from $`z=19`$, with 4 outputs (at redshifts $`z=`$1.00, 0.60, 0.25, and 0.00), when the corresponding $`M_{}`$ values are 50, 100, 200, and 400 particles, respectively. The following results are averaged over four independent realizations of the simulation, and the error bars represent the dispersion among the four realizations divided by $`\sqrt{3}`$ to yield the error on the mean.

Efstathiou et al. (1988) performed the first study on scale-free cosmological models using $`N`$-body simulations and showed that many quantities in scale-free models have well-defined scaling behaviors. In an Einstein-de Sitter scale-free model, the cumulative halo mass functions at redshifts $`z_1`$ and $`z_2`$ have the scaling behavior

$$N(M>M_2;z_2)=\left(\frac{1+z_2}{1+z_1}\right)^{6/(3+n_s)}N(M>M_1;z_1),M_1=\left(\frac{1+z_2}{1+z_1}\right)^{6/(3+n_s)}M_2.$$
(1)
The left panel of Figure 3 shows the cumulative halo mass function at different redshifts in our scale-free simulations (solid curves). The dotted curves in Figure 3 are derived by scaling the mass function at $`z=0`$ according to the above scaling rule. The scaling properties shown in the figure are generally very good. The high mass end of the mass function drops systematically below the self-similar scaling as $`M_{}`$ increases, a likely result of the finite simulation volume, which supresses the amplitude of fluctuations on the largest scales. Comparison of the first and last outputs suggests that this effect becomes noticeable above $`M2\times 10^{15}h^1M_{}`$ at $`z=0`$ (for our adopted $`\sigma _8=1.0`$ scaling), or a mass scale $``$ 0.1% of the total simulation mass (10% in length scale).

The $`N`$-point correlation functions $`\xi (r)`$ are expected to have a similarity transformation with the variable $`s=r(1+z)^{2/(3+n_s)}`$ for an Einstein-de Sitter universe (Efstathiou et al., 1988). The scaling rule for the correlation functions at redshifts $`z_1`$ and $`z_2`$ is then

$$\xi (r_2;z_2)=\xi (r_1;z_1),r_1=\left(\frac{1+z_2}{1+z_1}\right)^{2/(3+n_s)}r_2.$$
(2)
From Figure 3 and Figure 4, it is clear that mass correlation functions and halo correlation functions at different values of $`M/M_{}`$ obey the scaling rule quite well. There are some noticeable departures from self-similar scaling in the $`M_{}/2`$ panel at the two earliest output times (especially the earliest), when $`M_{}/2`$ corresponds to 25 and 50 particles, respectively.

The clustering of halos is biased relative to that of the matter. The halo bias factor $`b_h`$ can be defined via $`\xi _h(r)=b_h^2\xi _m(r)`$, where $`\xi _h`$ and $`\xi _m`$ are the two-point correlation functions of halos and matter, respectively. The relation between the clustering of halos and that of the matter has been extensively studied based on analytical models and numerical simulations (e.g., Cole & Kaiser 1989; Kashlinsky 1991; Mo & White 1996; Jing 1998, 1999; Porciani, Catelan, & Lacey 1999; Sheth & Lemson 1999; Sheth & Tormen 1999; Sheth, Mo, & Tormen 2001; and references therein). Mo & White (1996) give an analytic formula for the bias factor of halos of a given mass based on an extended Press-Schechter analysis (Press & Schechter, 1974). This formula is quite accurate for halos with $`M>M_{}`$. Jing (1998) empirically modifies the formula to also fit the $`N`$-body simulation results for halos of lower mass,

$$b_h(M,z)=\left(\frac{1}{2\nu ^4}+1\right)^{(0.060.02n_s)}\left(1+\frac{\nu ^21}{\delta _{c,0}}\right),\nu =\delta _c(z)/\sigma (M),$$
(3)
where $`\sigma (M)`$ is the r.m.s. fluctuation of the mass density at a mass scale $`M`$, $`\delta _c(z)`$ is the threshold density contrast for collapse of a homogeneous spherical overdense region at redshift $`z`$, and $`\delta _{c,0}1.69`$ for an Einstein-de Sitter universe. For $`n_s=1`$ here, $`\nu =(M/M_{})^{2/3}`$. We derive the square of the bias factor numerically as a function of $`M/M_{}`$ at different redshifts by averaging $`\xi _h(r)/\xi _m(r)`$ over comoving separations between 5$`h^1\mathrm{Mpc}`$ and 30$`h^1\mathrm{Mpc}`$ . In this range, the ratio is almost constant. We compare the computed bias factor with that given by Jing’s fitting formula in Figure 5. Two conclusions can be drawn immediately from the comparison. First, since bias factors at different redshifts overlap with each other when plotted as a function of $`M/M_{}`$, the simulation displays the correct scaling behavior. Second, the bias factors agree well with Jing’s fitting formula. At higher mass, the computed bias is somewhat lower than that given by Jing’s formula. In fact, the same trend can be found in Jing’s plot (see his Figure 2). We attribute this slight deviation in part to a sample volume effect and in part to systematic errors in the approximate analytic formula given by Mo & White (1996) (see Sheth, Mo, & Tormen 2001).

Both the GADGET simulation test and the scale-free model test assure that the numerical artifacts in the PM simulations have almost negligible effect on the statistics we care about in this paper. For example, the mass scales that we are concerned with are generally higher than the level where the PM and GADGET halo mass functions begin to deviate slightly. Moreover, the correlation functions are measured at scales $``$ 0.5 mesh cells, where the PM and GADGET results are consistent. The PM simulations reproduce the expected scaling behaviors in the scale-free model, which also demonstrates that we can measure masses of halos and their spatial distribution correctly based on these simulations. The scaling tests and code comparison suggest that there are some inaccuracies when the halo mass is lower than 50 particles or more than 0.1% of the total mass in the simulation, but even in these regimes there should be little impact on our conclusions because we compare different cosmologies evolved with the same code and numerical parameters. The PM code is therefore an ideal tool for our purpose in this paper.

## 3 Comparison of Halo Populations

### 3.1 Changing $`\mathrm{\Omega }_m`$ with $`\sigma _8`$, $`n_s`$, and $`\mathrm{\Gamma }`$ Fixed

We start with a simple case in which we fix the power spectrum (as linearly evolved to $`z=0`$, with $`n_s=1.0`$, $`\sigma _8=0.9`$, $`\mathrm{\Gamma }=0.20`$) and change $`\mathrm{\Omega }_m`$. In linear theory, the evolved fluctuations in these models are identical, and this independence of $`\mathrm{\Omega }_m`$ holds in the Zel’dovich approximation (Zel’dovich, 1970) as well as its extension, the adhesion approximation, as emphasized by Weinberg & Gunn (1990). In fact, it holds remarkably well into the fully nonlinear regime, as shown in the top and middle panels of Figure 6, where slices of particle distributions at $`z=0`$ for models with different $`\mathrm{\Omega }_m`$ but identical initial conditions are compared. These slices are taken from GADGET simulations that evolve the particle distribution from $`z=19`$. Particle distributions of three spatially flat models with different values of $`\mathrm{\Omega }_m`$ and one open model with $`\mathrm{\Omega }_m=0.2`$ are extremely similar to each other (though particle velocities are higher for higher $`\mathrm{\Omega }_m`$). Nusser & Colberg (1998) extend the Zel’dovich/adhesion analytic explanation for this similarity of evolved structure by showing that the equations of motion of a collisionless gravitating system of particles in an expanding universe can be put in a form with almost no dependence on $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$, if the linear perturbation growth factor is used as the time variable. The evolutionary histories of the four models in Figure 6 are different, but this difference is captured almost entirely by the dependence of the linear growth factor on redshift. The bottom two panels demonstrate this point, showing the particle distributions of the flat $`\mathrm{\Omega }_m=1.0`$ and $`\mathrm{\Omega }_m=0.2`$ models at the redshifts when the linear growth factor is 0.5, redshifts $`z=1.0`$ and $`z=1.74`$, respectively (we define the growth factor to be unity at $`z=0`$.) In most comparisons of N-body simulations with different $`\mathrm{\Omega }_m`$, the amplitude and/or shape of $`P(k)`$ is adjusted in concert with $`\mathrm{\Omega }_m`$. Such adjustments are well motivated by observational and theoretical considerations, but they mask the fact that the influence of $`\mathrm{\Omega }_m`$ at fixed $`P(k)`$ is extremely simple. Earlier examples of matched comparisons like those in Figure 6 include Figure 1 of Davis et al. (1985) and Figure 5 of Nusser & Colberg (1998).

Figure 6 suggests that halo mass functions for different models in this case should be nearly the same when measured as a function of particle number, i.e., mass divided by $`\mathrm{\Omega }_m`$. This expectation is confirmed in Figure 7, which shows mass functions (from PM simulations) of four spatially flat models. In the left panel, all four (cumulative) mass functions have a similar shape and only shift horizontally relative to each other. When the halo masses are divided by $`\mathrm{\Omega }_m`$, the four curves become nearly identical and overlap with the $`\mathrm{\Omega }_m=1`$ curve. To allow a closer inspection, in the right panel we plot the ratio of the scaled mass functions to that calculated using Jenkins et al.’s (2001) fitting formula for the central model. The maximum relative difference between the four scaled mass functions is only about 10%.

The 10% residual differences reflect slight changes in the characteristic densities and profiles of halos in different cosmologies. For Figure 7, we have identified halos with our standard FoF linking length parameter $`b=0.2`$, which picks out structures with mean overdensity $`\rho /\overline{\rho }200`$ (Davis et al. 1985; Lacey & Cole 1993), close to the value $`\rho /\overline{\rho }178`$ predicted for the post-virialization overdensity in the spherical collapse model for an $`\mathrm{\Omega }_m=1`$ universe (Gunn & Gott 1972; Peacock 1999). If we instead use $`b=0.16`$, roughly doubling the overdensity threshold ($`0.2^3/0.16^31.95`$), then halo masses drop by $``$20% as previously linked particles are unlinked, but the level of agreement among different cosmological models stay the same. Another fairly common procedure is to scale $`b^3`$ in proportion to the virial overdensity predicted by the spherical collapse model, thus using a different $`b`$ for each cosmological model. If we adopt this approach, taking the virial overdensities from Eke, Cole, & Frenk (1996), then the discrepancy of the scaled mass functions reverses sign and rises in amplitude to a maximum $`20\%`$. We conclude that the effect of a pure $`\mathrm{\Omega }_m`$ change is described by the simple scaling of halo masses to an accuracy $`10\%`$ for halos defined at fixed overdensity, that the discrepancies reflect the expected cosmology dependence of the characteristic virial overdensity, but that the magnitude of the cosmology correction predicted by the spherical collapse model is too large. Jenkins et al. (2001) also find better agreement among halo mass functions of different cosmological models for halos defined at constant overdensity rather than overdensities scaled according to the spherical collapse model. Henceforth we will use the fixed linking length $`b=0.2`$ for all models, but our results would not be substantially different if we used $`b=0.16`$ or scaled $`b`$ with cosmology.

The mass correlation functions $`\xi _m(r)`$ for these four models are identical at large scales and gradually depart from each other at small scales ($`r<1`$$`h^1\mathrm{Mpc}`$ ), where the higher overdensities for lower $`\mathrm{\Omega }_m`$ boost $`\xi _m(r)`$. We calculate the two-point correlation function for halos of mass $`M`$ by counting pairs of halos with masses in the range ($`M/\sqrt{2}`$, $`\sqrt{2}M`$). Figure 8 shows halo correlation functions $`\xi _h(r)`$ for four values of $`M/M_{}`$. In each case, the four $`\mathrm{\Omega }_m`$ models are indistinguishable from each other, as we would expect, since $`M_{}`$ scales with $`\mathrm{\Omega }_m`$, and comparisons at fixed $`M/M_{}`$ therefore remove any dependence on $`\mathrm{\Omega }_m`$. The difference in $`\xi _m(r)`$ arises from differences in spatial structure within the halos themselves (which are difficult to discern in Figure 6). We also calculate the halo bias factor from the mass and halo correlation functions (not shown). When plotted as a function of $`M/M_{}`$, the bias factors of the four models are almost the same, as expected, and they agree with Jing’s (1998) fitting formula to a similar degree as seen in Figure 5. The halo correlation function plummets at small $`r`$ because two halo centers cannot be separated by less than the sum of their virial radii. Halos for different $`\mathrm{\Omega }_m`$ have the same virial radius at a given $`M/M_{}`$ (in simulation terms, the scaling of $`M_{}`$ with $`\mathrm{\Omega }_m`$ comes from the particle mass, not a change in halo size), so even this exclusion signature is the same in different models.

Although mass functions and spatial clustering of halos are nearly identical as a function of $`M/M_{}`$, different values of $`\mathrm{\Omega }_m`$ produce different velocity fields. Figure 9 shows the mean pairwise radial (inward) velocity, $`v_{12}(𝐯_\mathrm{𝟏}𝐯_\mathrm{𝟐})(𝐫_\mathrm{𝟏}𝐫_\mathrm{𝟐})/|𝐫_\mathrm{𝟏}𝐫_\mathrm{𝟐}|`$, and the pairwise velocity dispersion, $`\sigma _{12}v_{12}^2v_{12}^2`$, at $`M=M_{}`$ and $`M=8M_{}`$. Here $`𝐯_i`$ and $`𝐫_i`$ ($`i=1,2`$) are the velocities and positions of a halo pair, and the average is over all halo pairs with separations around $`r`$. Both $`v_{12}(r)`$ and $`\sigma _{12}(r)`$ drop as $`\mathrm{\Omega }_m`$ drops. The mean pairwise radial velocity $`v_{12}(r)`$ increases for massive halos and increases much faster for pairs with small separations, a sign of the non-linear infall velocities induced by high mass halos. The velocity dispersion $`\sigma _{12}(r)`$ continues rising out to large separations, and its amplitude at large scales seems to be nearly independent of halo mass. Linear perturbation theory predicts a relation between the peculiar velocity and matter density fields,

$$𝐯(𝐱)f(\mathrm{\Omega }_m)\delta (𝐱)\frac{(𝐱^{}𝐱)}{|𝐱^{}𝐱|^3}d^3x^{},$$
(4)
where $`\delta (𝐱)`$ is the mass density contrast and $`f(\mathrm{\Omega }_m)\mathrm{\Omega }_m^{0.6}`$ (see, e.g., Peebles 1993). In the four cosmological models considered here, the mass density contrasts $`\delta (𝐱)`$ should be the same, so linear velocity fields simply scale with $`\mathrm{\Omega }_m^{0.6}`$. The dotted curves in Figure 9 show the effect of dividing $`v_{12}(r)`$ and $`\sigma _{12}(r)`$ by $`\mathrm{\Omega }_m^{0.6}`$; the agreement of these curves with each other and with the $`\mathrm{\Omega }_m=1`$ curve demonstrates that the impact of a pure $`\mathrm{\Omega }_m`$ change on halo peculiar velocities is well described by the linear theory scaling even on non-linear scales.

Figure 6 suggests that the agreement of appropriately scaled halo mass functions, correlation functions, and pairwise velocity statistics should hold at other epochs for which the linear growth factors are equal. Though we do not show the results here, we have verified explicitly that the matching of scaled statistical measures in Figures 7, 8, and 9 holds equally well at the redshifts when the four models have growth factors of 0.5 ($`z=1.0`$, 1.36, 1.50, and 1.74 for $`\mathrm{\Omega }_m=1.0,`$ 0.4, 0.3, 0.2). Thus, for models with the same $`P(k)`$ shape and the same $`\sigma _8`$ at $`z=0`$, the influence of $`\mathrm{\Omega }_m`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }`$ on the evolutionary history of the halo population is entirely encoded in the linear growth factor.

### 3.2 Cluster-normalized Models with Fixed $`P(k)`$ Shape

Changing $`\mathrm{\Omega }_m`$ on its own preserves the shape of the halo mass function but shifts the mass scale in proportion to $`\mathrm{\Omega }_m`$. This scaling implies that the threshold volume for collapse is independent of $`\mathrm{\Omega }_m`$, given the same power spectrum. If we require that the four cosmological models that differ in $`\mathrm{\Omega }_m`$ have halo mass functions that match in amplitude at some physical mass scale (rather than $`M/M_{}`$), we need to adjust the power spectrum in a way that builds more massive halos in low-$`\mathrm{\Omega }_m`$ models and fewer massive halos in high-$`\mathrm{\Omega }_m`$ models. In other words, we need a larger threshold volume for collapse in a low-$`\mathrm{\Omega }_m`$ model, so that its nonlinear mass $`M_{}`$ increases. Therefore, if the shape of the power spectrum is kept unchanged, a higher $`\sigma _8`$ is necessary in a lower $`\mathrm{\Omega }_m`$ model in order to compensate for the lower mean mass density.

Putting this idea in quantitative form, White, Efstathiou, & Frenk (1993) argued that the observed abundances of massive clusters of galaxies impose the constraint $`\sigma _80.57\mathrm{\Omega }_m^{0.56}`$, with little dependence on the assumed shape of $`P(k)`$. Many authors have revisited, refined and recalibrated this “cluster normalization” constraint (see Pierpaoli, Scott, & White 2001 and references therein). One widely used formulation is that of Eke, Cole, & Frenk (1996), who found $`\sigma _8=(0.52\pm 0.04)\mathrm{\Omega }_m^{0.52+0.13\mathrm{\Omega }_m}`$ for a spatially flat universe. Inspired by these results, we now consider a sequence of cluster normalized models in which we keep the shape of $`P(k)`$ fixed but scale the amplitude as $`\sigma _8=0.9(\mathrm{\Omega }_m/0.3)^{0.5}`$.

Figure 10 shows cumulative halo mass functions for these cluster normalized models. Mass functions of models with $`\mathrm{\Omega }_m=`$0.2, 0.3, and 0.4 match in amplitude at a cluster scale ($`M5\times 10^{14}h^1M_{}`$). The mass function of the $`\mathrm{\Omega }_m=1`$ model crosses the others at a somewhat lower mass scale, $`M2\times 10^{14}h^1M_{}`$. The discrepancy in matched mass scale just shows that the $`\sigma _8\mathrm{\Omega }_m^{0.5}`$ scaling does not hold all the way to $`\mathrm{\Omega }_m=1`$, and the Eke, Cole, & Frenk (1996) formula indeed implies that we should adopt a somewhat higher $`\sigma _8`$ for $`\mathrm{\Omega }_m=1`$. However, the crucial result for our purposes is that the mass functions of these cluster normalized models have systematically different shapes, and there is no value of $`\sigma _8`$ we could choose that would make the $`\mathrm{\Omega }_m=1`$ model match the others at more than a single mass scale. Over the smaller range $`\mathrm{\Omega }_m=0.20.4`$, the cluster normalization condition keeps the halo mass functions reasonably well aligned over about a decade in mass, but by $`M10^{13}h^1M_{}`$ the halo space densities of the $`\mathrm{\Omega }_m=0.2`$ and 0.4 models differ from those of the central model by $`30\%`$, and from each other by nearly a factor of two.

The reason that cluster normalization does not preserve the shape of the halo mass function is that raising $`\sigma _8`$ increases $`M_{}`$ but simultaneously drives down the space density of $`M_{}`$ clusters by increasing the scale of non-linearity. For models with the same $`P(k)`$ shape and the same value of $`M_{}`$, the space density of halos at fixed $`M/M_{}`$ is proportional to $`\mathrm{\Omega }_m`$ (see equation (7) below). The mass functions of cluster normalized models with different $`\mathrm{\Omega }_m`$ must therefore cross at different values of $`M/M_{}`$, and because the shapes remain the same as a function of $`M/M_{}`$, they cannot stay aligned as a function of $`M`$.

In these cluster normalized models, the amplitude of $`\xi _m(r)`$ increases as $`\mathrm{\Omega }_m`$ decreases, as shown in the upper left panel of Figure 11. At large scales, where the mass density field is still linear, $`\xi _m(r)`$ is just proportional to the amplitude of its Fourier transform, the linear power spectrum, and hence to $`\sigma _8^2`$. The roughly constant logarithmic offset in Figure 11 shows that this multiplicative scaling holds approximately even into the strongly non-linear regime. However, the remaining panels of Figure 11 show that the halo correlation functions are very similar at a fixed physical mass. Qualitatively, this similarity makes sense, since a given mass scale corresponds to a higher value of $`M/M_{}`$ when $`\mathrm{\Omega }_m`$ is higher, and the correspondingly higher halo bias tends to compensate for the lower mass clustering amplitude. However, it is impressive just how well this cancellation between mass clustering and halo bias works quantitatively. The halo exclusion signature does shift to larger scale for lower $`\mathrm{\Omega }_m`$ because of the larger virial radii needed to compensate the lower mass density.

Figure 12 shows that the mean pairwise velocity and pairwise dispersion increase with $`\mathrm{\Omega }_m`$ in this cluster normalized sequence, but the dependence is weak, especially at large separations. As in Figure 9, more massive halos have higher pairwise velocities, especially at small separations. The weak dependence on $`\mathrm{\Omega }_m`$ at large scales can be understood with the help of equation (4). In linear theory, velocity is proportional to $`f(\mathrm{\Omega }_m)\sigma _8`$, which, when combined with the cluster normalization condition, reduces to approximately $`\mathrm{\Omega }_m^{0.1}`$. Non-linear evolution at small separations amplifies the difference between models, but only slightly.

### 3.3 Matching Mass Functions

By varying the amplitude of the power spectrum with $`\mathrm{\Omega }_m`$, we can match the amplitudes of the cumulative halo mass functions at a given mass scale. However, this combination of parameter changes does not guarantee a match between the shapes of the mass functions at that mass scale. In order to match both the amplitude and the shape of the halo mass functions, we must adjust the relative amplitudes of fluctuations at different scales, or, in other words, the shape of $`P(k)`$. We therefore need to determine the appropriate combinations of parameters ($`\sigma _8`$ and parameters related to the shape of the power spectrum) for each $`\mathrm{\Omega }_m`$ model so that all models produce halo mass functions that agree both in amplitude and in shape.

#### 3.3.1 Analytic Solution to Parameter Combinations

Instead of searching the parameter space to obtain the right power spectrum combinations, we make use of the analytic form of the halo mass function. There are three commonly used analytic formulae for halo mass functions. Press & Schechter (1974) developed the first theoretical framework based on the spherical collapse model (also see Bond et al. 1991) and provided the first theoretical formula (PS formula) for the halo mass function, which is still widely used. The PS formula is known to underestimate the halo abundance in the high mass tail (see Jenkins et al. 2001 and references therein). Sheth & Tormen (1999) give an analytic formula (ST formula) obtained by fitting to the results of $`N`$-body simulations, which turns out to be consistent with the ellipsoidal collapse model (Sheth, Mo, & Tormen, 2001). Based on more simulation results, Jenkins et al. (2001) provide the most accurate fitting formula for halo mass functions. However, this formula includes an absolute value of a function, which makes it difficult to use for our present purpose: obtaining an analytic solution of parameter combinations that lead to a good match in halo mass functions. We therefore use the slightly less accurate ST formula.

In the notation of Jenkins et al. (2001), the mass function is defined in terms of the dimensionless function

$$f(\nu )\frac{M}{\rho _0}\frac{\mathrm{d}n(M)}{\mathrm{d}\mathrm{ln}\sigma ^1},$$
(5)
and the ST formula for the halo mass function can be expressed as (Sheth & Tormen, 1999)

$$f(\nu )=A\sqrt{\frac{2a}{\pi }}\left[1+(a\nu ^2)^p\right]\mathrm{exp}(a\nu ^2/2),\nu =\delta _c/\sigma (M),$$
(6)
where $`n(M)`$ is the abundance of halos with mass less than $`M`$, $`\delta _c1.69`$ is the threshold density contrast for collapse, $`\sigma (M)`$ is the r.m.s. fluctuation of the mass density at a mass scale $`M`$, $`\rho _0`$ is the mean density of the current universe, $`A=0.3222`$, $`a=0.707`$, and $`p=0.3`$. The mass function can be written as

$$\frac{\mathrm{d}n}{\mathrm{d}\mathrm{ln}M}=\frac{\mathrm{d}n}{\mathrm{d}\mathrm{ln}\sigma ^1}\frac{\mathrm{d}\mathrm{ln}\sigma ^1}{\mathrm{d}\mathrm{ln}M}=\alpha \frac{\rho _0}{M}f(\nu ),$$
(7)
where $`\alpha =\mathrm{d}\mathrm{ln}\sigma ^1/\mathrm{d}\mathrm{ln}M=(3+n_{\mathrm{eff}})/6`$ characterizes the local index $`n_{\mathrm{eff}}`$ of the power spectrum and is determined by the shape of the power spectrum. The slope of the mass function is

$$\frac{\mathrm{d}}{\mathrm{d}\mathrm{ln}M}\left(\frac{\mathrm{d}n}{\mathrm{d}\mathrm{ln}M}\right)=\frac{\mathrm{d}n}{\mathrm{d}\mathrm{ln}M}\left(\alpha \frac{\mathrm{d}\mathrm{ln}f}{\mathrm{d}\mathrm{ln}\nu }1+\frac{\mathrm{d}\mathrm{ln}\alpha }{\mathrm{d}\mathrm{ln}M}\right).$$
(8)
We use a power-law to locally approximate the power spectrum (so that the last term of the right-hand side of equation (8) reduces to zero).

We require that the halo mass function of a given cosmological model match as closely as possible that of the central model at some mass scale. In other words, at this mass scale, both the amplitudes (eq. 7) and the slopes (eq. 8) of the two halo mass functions should be equal. We thus solve these two equations for two unknowns ($`\alpha `$,$`\nu `$) for each cosmological model. All the information contained in the power spectrum is fully embedded in ($`\alpha `$,$`\nu `$). Therefore, by solving for these two parameters at a given mass scale, we can also determine the normalization $`\sigma _8`$ and shape of the power spectrum.

The power spectrum can be expressed as $`P(k)k^{n_s}T^2(k;\mathrm{\Gamma })`$, where $`n_s`$ is the spectral index of the inflationary power spectrum and $`T(k;\mathrm{\Gamma })`$ is the transfer function with shape parameter $`\mathrm{\Gamma }`$. A change in either $`n_s`$ or $`\mathrm{\Gamma }`$ leads to a change in the shape of the power spectrum. We first keep $`\mathrm{\Gamma }`$ fixed and only tilt $`P(k)`$. Since $`\alpha =(3+n_{\mathrm{eff}})/6`$ and $`n_{\mathrm{eff}}`$ is simply the sum of $`n_s`$ and an index given by the transfer function, it is straightforward to obtain the change in $`n_s`$ with respect to that of the central model: $`\mathrm{\Delta }n_s=6\mathrm{\Delta }\alpha `$. Once the shape of $`P(k)`$ is determined, the normalization $`\sigma _8`$ can be obtained from $`\nu `$. Assuming that mass functions are matched against that of the central model at a cluster mass scale $`M=5\times 10^{14}h^1M_{}`$, we calculate $`n_s`$ and $`\sigma _8`$ as a function of $`\mathrm{\Omega }_m`$. The left panel of Figure 13 shows that both $`n_s`$ and $`\sigma _8`$ increase sharply towards low $`\mathrm{\Omega }_m`$ and drop slowly towards high $`\mathrm{\Omega }_m`$. In order to match with the halo mass function of the central model, we need $`n_s=`$1.33, 0.82, 0.42 and $`\sigma _8=`$1.16, 0.78, 0.55 for models with $`\mathrm{\Omega }_m=`$0.2, 0.4, and 1.0, respectively.

Alternatively, we can keep the spectral index $`n_s`$ fixed ($`n_s`$=1 in our calculation) and change the shape parameter $`\mathrm{\Gamma }`$. The right panel of Figure 13 shows $`\mathrm{\Gamma }`$ as a function of $`\mathrm{\Omega }_m`$ when halo mass functions are matched at $`M=5\times 10^{14}h^1M_{}`$. The shape parameter $`\mathrm{\Gamma }`$ also has a sharp increase towards low $`\mathrm{\Omega }_m`$ and a slow drop towards high $`\mathrm{\Omega }_m`$. For models with $`\mathrm{\Omega }_m=`$0.2, 0.4, and 1.0, solutions to $`\mathrm{\Gamma }`$ are 0.35, 0.14, and 0.06, respectively. In this case, the normalization $`\sigma _8`$ of each model is almost identical to that of the tilted model, which indicates that the two ways of changing the shape of the power spectrum are equivalent at this mass scale. Although pure-$`n_s`$ and pure-$`\mathrm{\Gamma }`$ changes alter the power spectrum in different ways, they produce power spectra (and thus mass functions) that match closely over a fairly large dynamic range, covering most of the regime in which $`\sigma `$ is of order unity. The near degeneracy can be broken by examining the low mass end of the halo mass functions. More generally, we could allow both $`n_s`$ and $`\mathrm{\Gamma }`$ to change. In that case, for each value of $`\mathrm{\Omega }_m`$, there would be a 1-dimensional locus in the $`n_s\mathrm{\Gamma }`$ plane along which models would satisfy the requirement of matching the central model’s halo mass function at some mass scale. However, we expect that results from this more general case would be similar to those of the two simple cases.

#### 3.3.2 Simulation Results

Using the power spectrum parameters given by the above analytic solutions, we run simulations and make comparisons with the central model. Here we only show results for different tilt models, since models with different values of $`\mathrm{\Gamma }`$ give similar results.

Figure 14 shows cumulative halo mass functions for the four cosmological models. The mass functions agree with each other quite well over a large mass range (from $`10^{13}h^1M_{}`$ to $`10^{15}h^1M_{}`$), with only $`10\%`$ relative differences. We do not expect the cumulative mass functions to match exactly both in amplitude and in slope at $`M=5\times 10^{14}h^1M_{}`$, for three reasons. First, the ST formula is not a perfect fit to the halo mass functions in the simulations. Second, we locally approximate the power spectrum as a power-law in order to obtain analytic solutions for the power spectrum parameters. Third, we require a match for differential mass functions instead of cumulative mass functions. Nevertheless, the cumulative mass functions trace each other remarkably well. The non-linear mass scales are quite different for these four models: $`M_{}=3.10\times 10^{13}h^1M_{}`$, $`1.03\times 10^{13}h^1M_{}`$, $`3.66\times 10^{12}h^1M_{}`$, and $`7.50\times 10^8h^1M_{}`$, for $`\mathrm{\Omega }_m=0.2`$, 0.3, 0.4, 1.0, respectively.

Figure 15 shows two-point correlation functions of the mass and halos. Mass is more strongly clustered in low $`\mathrm{\Omega }_m`$ models, having a correlation function that is steeper and higher in amplitude than in higher $`\mathrm{\Omega }_m`$ models. However, the clustering of halos exhibits the opposite trend. Halos in the $`\mathrm{\Omega }_m=1`$ model are highly biased with respect to those of other models, as one might guess from the low nonlinear mass of this model, and this bias more than compensates for the weaker mass clustering. Conversely, at the mass scale $`M=1.03\times 10^{13}h^1M_{}`$, shown in the upper right panel of Figure 15, halos in the $`\mathrm{\Omega }_m=0.2`$ model are actually anti-biased with respect to the mass. Figure 11 shows that for cluster normalized models with fixed $`P(k)`$ shape, the effects of bias almost exactly cancel the change in $`\sigma _8`$, leaving $`\xi _h(r)`$ nearly independent of $`\mathrm{\Omega }_m`$. Since our mass function matching approach imposes a redder $`P(k)`$ shape (lower $`n_s`$ or $`\mathrm{\Gamma }`$) for higher $`\mathrm{\Omega }_m`$, and increased large scale power amplifies the importance of bias, it is not surprising that bias over compensates mass clustering changes in these models, leaving $`\xi _h(r)`$ steadily dependent on $`\mathrm{\Omega }_m`$.

Figure 16 shows the mean pairwise radial velocity and pairwise velocity dispersion for the four cosmological models. As in the case of cluster normalized models, models with higher $`\mathrm{\Omega }_m`$ have larger mean pairwise velocities and velocity dispersions than lower $`\mathrm{\Omega }_m`$ models. The dependence on $`\mathrm{\Omega }_m`$ is stronger and more systematic than that in Figure 12, again a sign of the redder power spectrum of higher $`\mathrm{\Omega }_m`$ models.

## 4 Summary and Discussion

We have examined how changes to the matter density parameter $`\mathrm{\Omega }_m`$ and the shape and amplitude of the linear power spectrum $`P(k)`$ affect the halo mass function, the two-point halo correlation function, and the first and second moments of the halo pairwise velocity distribution. A change in $`\mathrm{\Omega }_m`$, at fixed $`P(k)`$, simply shifts the halo mass scale. Therefore, if halo masses are scaled in proportion to $`\mathrm{\Omega }_m`$, halo populations of different $`\mathrm{\Omega }_m`$ models have identical mass functions and clustering properties. However, mean pairwise velocities and pairwise velocity dispersions, which scale as $`\mathrm{\Omega }_m^{0.6}`$ on large scales (and fairly far into the non-linear regime) break this degeneracy. A change of the power spectrum amplitude ($`\sigma _8`$) along with $`\mathrm{\Omega }_m`$, maintaining a cluster normalization condition $`\sigma _8\mathrm{\Omega }_m^{0.5}`$, produces halo populations that have nearly identical spatial clustering at a fixed physical mass scale and systematic but small differences in the mean pairwise velocities and pairwise velocity dispersions. However, while these cluster normalized models have the same halo space density at $`M`$ several $`\times 10^{14}M_{}`$, the shapes of their mass functions are systematically different. If we further allow the shape of the power spectrum to change (either by tilting it or by changing its shape parameter $`\mathrm{\Gamma }`$), we can produce halo populations whose mass functions match very well both in amplitude and in shape, over a large mass range. However, the changes to the power spectrum shape cause correlation functions, mean pairwise velocities, and velocity dispersions for halos in the same mass range to differ. We conclude that the halo populations produced by distinct cosmological models are not degenerate. If they are indistinguishable by one statistic, they can be told apart using another statistic.

Our results imply that a perfect knowledge of the halo population and its properties would allow us to pin down the underlying cosmological parameters – specifically the value of $`\mathrm{\Omega }_m`$ and the shape and amplitude of $`P(k)`$ – even before bringing in constraints from other cosmological measurements. However, a galaxy redshift survey detects galaxies rather than halos. Could changes to the halo occupation distribution (HOD) mask the differences in the halo populations of different cosmological models? Models that differ only in $`\mathrm{\Omega }_m`$ lead to the same halo mass function and halo clustering in terms of $`M/M_{}`$, but with $`M_{}\mathrm{\Omega }_m`$. Their halos could be thus be populated the same way as a function of $`M/M_{}`$ to produce indistinguishable galaxy spatial clustering. However, any dynamically sensitive statistics – virial masses of clusters and groups, pairwise velocity dispersions, the parameter combination $`\mathrm{\Omega }_m^{0.6}/b`$ inferred from redshift-space distortions, the galaxy-mass correlation function from galaxy-galaxy lensing – would distinguish models with different values of $`\mathrm{\Omega }_m`$ immediately. Velocity bias within halos could mask some of these changes, but not all of them (Berlind & Weinberg, 2002). Cluster-normalized models with different $`\mathrm{\Omega }_m`$ and the same $`P(k)`$ shape have different halo mass functions, and they would therefore require a different $`P(N|M)`$ in order to keep the galaxy density fixed. This change would likely cause differences in the galaxy clustering even though the halo clustering is similar, since galaxy clustering is highly sensitive to $`P(N|M)`$ (Berlind & Weinberg 2002 and references therein). Changing the shape of $`P(k)`$ by the amount required to match mass function shapes makes substantial differences to the halo clustering and velocities, which seem unlikely to be masked by a change in the HOD.

These speculations suggest that distinct cosmological models, which produce nondegenerate halo populations, cannot have galaxy populations that are indistinguishable in every spatial and velocity statistic, even if one allows complete freedom in the way that galaxies occupy these halos. We reserve a detailed investigation of this issue for future work. If our optimistic speculation holds true, then high precision measurements of galaxy clustering and galaxy-galaxy lensing should impose strong constrains on cosmological models without reliance on a priori models of galaxy bias.

We thank Volker Springel for advice on GADGET and Changbom Park for the use of his PM code. DW thanks the Institute for Advanced Study for its hospitality and the Ambrose Monell Foundation for its financial support during the final stages of this work. This work was supported by NSF grant AST-0098584.


text:
# Bose-enhanced chemistry: Amplification of selectivity in the dissociation of molecular Bose-Einstein condensates

A significant effort is currently underway with the goal of producing a stable Bose-Einstein condensate in a vapor of molecules. Several groups are pursuing schemes in which atomic BECs are converted to molecular BECs via photoassociation  and/or Feshbach resonance . Progress is also being made towards forming a molecular BEC by direct evaporation of molecules in a magnetic trap. In this case sympathetic buffer-gas cooling is used to load the trap , as opposed to the optical techniques used in the alkali-metal atomic vapors.

Once the technical problems are solved, the realization of molecular condensates will open a new field of quantum-degenerate coherent chemistry. The most striking new feature of the resulting molecular dynamics would be the unprecedented role of nonlinear collective effects, i.e. Bose enhancement and Pauli blocking, induced by the $`N`$-body quantum statistics. The probability of any single boson to undergo a transition is enhanced by the number of identical particles occupying the final state, whereas fermions would be subject to Pauli blocking of the output channels. Thus any process which retains a critical phase-space density would greatly differ from that of the isolated single-particle process which dominates ’thermal’ gas-phase chemistry .

The role of Bose stimulation in the association and dissociation of a molecular BEC has already been studied theoretically in the regime where a single atomic mode is coupled to a single molecular mode . In this situation the entire condensate is theoretically predicted to undergo large amplitude coherent oscillations between atoms and molecules . The Bose-enhanced oscillation frequency is $`\sqrt{N}\mathrm{\Omega }`$, where $`\mathrm{\Omega }`$ is the single-particle atom-molecule coupling frequency and $`N`$ is the total number of condensate particles. While this prediction is significantly modified by the inclusion of quantum-field effects , the collective evolution remains strikingly different from the single-pair (non-degenerate) dynamics.

The purpose of this Letter is to extend these ideas to the case where more complicated molecules are dissociated into a multiplicity of output channels. In such multi-channel quantum-degenerate processes there are two important competing effects: Bose stimulation, which leads to exponential growth of the output channel populations, and condensate depletion, in which the finite number of initial molecules leads to competition between output channels. If certain channels have a slight advantage in coupling strength, they will deplete the condensate before other channels attain a macroscopic population. Related effects have been experimentally observed in the phenomenon of superradiant Rayleigh scattering of laser light from an elongated BEC , indicating that there are no overwhelming obstacles to obtaining winner-takes-all scenarios in BEC systems.

Studying the simplest example of the dissociation of a triatomic ABC molecule into AB+C or A+BC pairs, we show that the combined effects of Bose stimulation and competition between output channels for a finite initial number of molecules leads to an amplification of selectivity which can dramatically enhance any non-degenerate control scheme. Relatively small controllability of the single-particle process translates to near complete controllability of the Bose-stimulated $`N`$-body process. Thus, collective amplification and mode competition effects can be employed to achieve the central goal of the fields of optimal- and coherent-control .

We consider a BEC of triatomic $`ABC`$ molecules. The ground-state molecules are resonantly coupled by a stimulated Raman transition to two open dissociation channels: $`AB+C`$ ($`AB`$ molecules and $`C`$ atoms) and $`A+BC`$ ($`A`$ atoms and $`BC`$ molecules), as depicted in Fig. 1. After adiabatic elimination of the intermediate $`ABC^{}`$ state the $`N`$-body second-quantized Hamiltonian in the interaction representation and the rotating wave approximation, reads

$$\widehat{H}=\frac{\mathrm{}}{2}\widehat{\mathrm{\Phi }}(𝐫)\left[\mathrm{\Omega }_1\widehat{\mathrm{\Psi }}_{AB}^{}(𝐫)\widehat{\mathrm{\Psi }}_C^{}(𝐫)+\mathrm{\Omega }_2\widehat{\mathrm{\Psi }}_A^{}(𝐫)\widehat{\mathrm{\Psi }}_{BC}^{}(𝐫)\right]+h.c.$$
(1)
where $`\mathrm{\Omega }_1,\mathrm{\Omega }_2`$ are the coupling frequencies into the $`AB+C`$ and $`A+BC`$ dissociation channels, respectively and $`\widehat{\mathrm{\Phi }}`$, $`\widehat{\mathrm{\Psi }}_{AB}`$, $`\widehat{\mathrm{\Psi }}_C`$, $`\widehat{\mathrm{\Psi }}_A`$, $`\widehat{\mathrm{\Psi }}_{BC}`$ are the bosonic field operators for $`ABC`$, $`AB`$, $`C`$, $`A`$, $`BC`$, respectively.

Using the Hamiltonian (1), the Heisenberg equations of motion for these operators are

$`i{\displaystyle \frac{d}{dt}}\widehat{\mathrm{\Phi }}`$ $`=`$ $`\mathrm{\Omega }_1\widehat{\mathrm{\Psi }}_{AB}\widehat{\mathrm{\Psi }}_C+\mathrm{\Omega }_2\widehat{\mathrm{\Psi }}_A\widehat{\mathrm{\Psi }}_{BC}`$ (2a)
$`i{\displaystyle \frac{d}{dt}}\widehat{\mathrm{\Psi }}_{AB}`$ $`=`$ $`\mathrm{\Omega }_1\widehat{\mathrm{\Phi }}\widehat{\mathrm{\Psi }}_C^{}`$ (2b)
$`i{\displaystyle \frac{d}{dt}}\widehat{\mathrm{\Psi }}_C^{}`$ $`=`$ $`\mathrm{\Omega }_1^{}\widehat{\mathrm{\Phi }}^{}\widehat{\mathrm{\Psi }}_{AB}`$ (2c)
$`i{\displaystyle \frac{d}{dt}}\widehat{\mathrm{\Psi }}_A`$ $`=`$ $`\mathrm{\Omega }_2\widehat{\mathrm{\Phi }}\widehat{\mathrm{\Psi }}_{BC}^{}`$ (2d)
$`i{\displaystyle \frac{d}{dt}}\widehat{\mathrm{\Psi }}_{BC}^{}`$ $`=`$ $`\mathrm{\Omega }_2^{}\widehat{\mathrm{\Phi }}^{}\widehat{\mathrm{\Psi }}_A`$ (2e)
Equations (2) are closely related to the quantum equations of motion of a multimode laser , in that the $`AB+C`$ dissociation channel (Eqs. (2b,2c)) is only coupled to the $`A+BC`$ channel (Eqs. (2d,2e)) via the molecular reservoir $`\widehat{\mathrm{\Phi }}`$. In order to demonstrate the collective amplification of selectivity, we use a five-mode model, neglecting the spatial variation of the various fields. This approximation, which assumes there is only one available state per specie, is valid for a zero temperature $`ABC`$ condensate at exact resonance with the dissociation threshold of the $`AB+C`$ and $`A+BC`$ channels. The field operators of Eqs. (2) are thus taken to represent the particle creation and annihilation operators of the respective modes. In the undepleted pump approximation , the intense $`\widehat{\mathrm{\Phi }}`$ field is replaced by a constant c-number $`\mathrm{\Phi }\sqrt{N}`$. Assuming without loss of generality that $`\mathrm{\Omega }_1=\mathrm{\Omega }_1^{}`$ and $`\mathrm{\Omega }_2=\mathrm{\Omega }_2^{}`$, the initial time evolution of the operators $`\widehat{\mathrm{\Psi }}_{AB},\widehat{\mathrm{\Psi }}_C,\widehat{\mathrm{\Psi }}_A`$, and $`\widehat{\mathrm{\Psi }}_{BC}`$ is thus given as

$`\widehat{\mathrm{\Psi }}_{AB}(t)`$ $`=`$ $`\widehat{\mathrm{\Psi }}_{AB}(0)\mathrm{cosh}\omega _1ti\widehat{\mathrm{\Psi }}_C^{}(0)\mathrm{sinh}\omega _1t,`$ (3a)
$`\widehat{\mathrm{\Psi }}_C^{}(t)`$ $`=`$ $`\widehat{\mathrm{\Psi }}_C^{}(0)\mathrm{cosh}\omega _1t+i\widehat{\mathrm{\Psi }}_{AB}(0)\mathrm{sinh}\omega _1t,`$ (3b)
$`\widehat{\mathrm{\Psi }}_A(t)`$ $`=`$ $`\widehat{\mathrm{\Psi }}_A(0)\mathrm{cosh}\omega _2ti\widehat{\mathrm{\Psi }}_{BC}^{}(0)\mathrm{sinh}\omega _2t,`$ (3c)
$`\widehat{\mathrm{\Psi }}_{BC}^{}(t)`$ $`=`$ $`\widehat{\mathrm{\Psi }}_{BC}^{}(0)\mathrm{cosh}\omega _2t+i\widehat{\mathrm{\Psi }}_A(0)\mathrm{sinh}\omega _2t,`$ (3d)
where $`\omega _1=\sqrt{N}\mathrm{\Omega }_1`$ and $`\omega _2=\sqrt{N}\mathrm{\Omega }_2`$. Assuming $`\widehat{\mathrm{\Psi }}_{AB}(0)=\widehat{\mathrm{\Psi }}_C(0)=\widehat{\mathrm{\Psi }}_A(0)=\widehat{\mathrm{\Psi }}_{BC}(0)=0`$ (no initial product coherence), the average numbers of particles $`n_j=\mathrm{\Psi }_j^{}\mathrm{\Psi }_j`$ in the modes $`j=AB,C,A,BC`$, are given, respectively, as

$`n_{AB}(t)`$ $`=`$ $`n_{AB}(0)\mathrm{cosh}^2\omega _1t`$ (4a)
$`+(1+n_C(0))\mathrm{sinh}^2\omega _1t,`$ (4b)
$`n_C(t)`$ $`=`$ $`n_C(0)\mathrm{cosh}^2\omega _1t`$ (4c)
$`+(1+n_{AB}(0))\mathrm{sinh}^2\omega _1t,`$ (4d)
$`n_A(t)`$ $`=`$ $`n_A(0)\mathrm{cosh}^2\omega _2t`$ (4e)
$`+(1+n_{BC}(0))\mathrm{sinh}^2\omega _2t,`$ (4f)
$`n_{BC}(t)`$ $`=`$ $`n_{BC}(0)\mathrm{cosh}^2\omega _2t`$ (4g)
$`+(1+n_A(0))\mathrm{sinh}^2\omega _2t`$ (4h)
Starting with a pure ABC condensate $`n_{AB}(0)=n_C(0)=n_A(0)=n_{BC}(0)=0`$, the average product populations grow as

$`n_{AB}(t)=n_C(t)=\mathrm{sinh}^2\sqrt{N}\mathrm{\Omega }_1t,`$ (5a)
$`n_A(t)=n_{BC}(t)=\mathrm{sinh}^2\sqrt{N}\mathrm{\Omega }_2t.`$ (5b)
The initial time evolution of Eqs. (5) corresponds to an initial nonexponential (quadratic rather than linear in time) spontaneous dissociation process, followed by an exponential stimulated process . We note that exactly the same equations have been routinely used for over three decades, to describe the phenomenon of parametric superfluorescence , i.e. the parametric amplification of noise photons, in quantum optics. The feature which produces the striking selectivity of Bose-enhanced chemistry is that any small difference between the two rate coefficients, $`\mathrm{\Omega }_2\mathrm{\Omega }_1`$, is exponentially amplified by the stimulated dissociation, so that

$$N_2(t)/N_1(t)e^{2\sqrt{N}(\mathrm{\Omega }_2\mathrm{\Omega }_1)t}.$$
(6)
The linearized model of Eqs. (3) is only valid as long as the population of the triatomic mode $`\widehat{\mathrm{\Phi }}^{}\widehat{\mathrm{\Phi }}`$ is large and the effect of its depletion on the products population growth is negligible. In order to go beyond this approximation we solve the full $`N`$-body problem numerically. Starting with no products and fixing the total number of particles $`N`$, the accessible Hilbert space is restricted to Fock states of the type $`|n_{ABC},n_{AB},n_C,n_A,n_{BC}`$ with $`n_{AB}=n_C`$, $`n_A=n_{BC}`$ and $`n_{ABC}=N(n_A+n_C)`$. Using these states to represent the Hamiltonian (1) we solve numerically the $`N`$-body Schrödinger equation,

$`i{\displaystyle \frac{d}{dt}}c_{n_1,n_2}=`$ $`\mathrm{\Omega }_1\sqrt{N(n_1+n_2)}(n_1+1)c_{n_1+1,n_2}`$ (7)
$`+\mathrm{\Omega }_1^{}\sqrt{N(n_1+n_2)+1}n_1c_{n_11,n_2}^{}`$ (8)
$`+\mathrm{\Omega }_2\sqrt{N(n_1+n_2)}(n_2+1)c_{n_1,n_2+1}`$ (9)
$`+\mathrm{\Omega }_2^{}\sqrt{N(n_1+n_2)+1}n_1c_{n_1,n_21}^{},`$ (10)
for the dynamics of the probability amplitudes

$$c_{n_1,n_2}(t)=N(n_1+n_2),n_1,n_1,n_2,n_2|\mathrm{\Psi }.$$
(11)
In Fig. 2, we plot the expectation values $`n_C`$, $`n_A`$, corresponding to the populations in the $`AB+C`$ channel and the $`A+BC`$ channel, respectively, as a function of time, for $`\mathrm{\Omega }_2=2\mathrm{\Omega }_1`$. When $`N=1`$ (Fig. 2(a)) the ratio between the two channel populations is, as expected for the single-particle process $`n_A/n_C=|\mathrm{\Omega }_2/\mathrm{\Omega }_1|^2=4`$. However, for $`N=1000`$ (Fig. 2(b)), the Bose stimulated population ratio is $`n_A/n_C>40`$. Thus, the reaction outcome is highly dependent on the total number of particles, becoming more selective as $`N`$ increases.

Our numerical calculations are limited by available computation power to $`N1000`$. However, the number of particles in current BEC experiments in trapped dilute alkali gases is sometimes as high as $`10^7`$, suggesting a far greater selectivity enhancement. In Fig. 3 we plot the channel population ratio when the leading channel population reaches a maximum, as a function of $`N`$ for various values of $`\mathrm{\Omega }_2/\mathrm{\Omega }_1`$. The curves are linear fits to the numerical data, demonstrating a power-law dependence of the selectivity as $`N^\alpha `$. The power $`\alpha `$ increases with the rate coefficient ratio $`\mathrm{\Omega }_2/\mathrm{\Omega }_1`$, as shown in Fig. 4.

In order to analytically estimate the dependence of the selectivity on the total particle number, we take a perturbative approach where we first obtain the time $`t_d`$ at which the leading channel depletes the condensate, and then evaluate the channel population ratio at $`t_d`$. Assuming a purely exponential growth, the depletion time is $`t_d=\mathrm{log}N/2\sqrt{N}\mathrm{\Omega }_2`$. Substituting into Eq. (6), we find that the selectivity is enhanced as $`N^\alpha `$, where

$$\alpha =1\mathrm{\Omega }_1/\mathrm{\Omega }_2,$$
(12)
so that $`\mathrm{log}[N_2(t_d)/N_1(t_d)]\alpha \mathrm{log}N`$. This result is confirmed by the numerical calculations in Fig. 3, giving $`\alpha `$ values of 0.17, 0.32, and 0.49 for $`\mathrm{\Omega }_2/\mathrm{\Omega }_1=1.25`$, $`\mathrm{\Omega }_2/\mathrm{\Omega }_1=1.5`$, and $`\mathrm{\Omega }_2/\mathrm{\Omega }_1=2`$, in good agreement with the predicted values of 0.2, 0.33, and 0.5, respectively. Extrapolating to $`N=10^6`$, we find population ratios that are greater than $`13`$, $`100`$, and $`1200`$ for $`\mathrm{\Omega }_2/\mathrm{\Omega }_1=1.25`$, $`\mathrm{\Omega }_2/\mathrm{\Omega }_1=1.5`$, and $`\mathrm{\Omega }_2/\mathrm{\Omega }_1=2`$, respectively (as compared to single particle ratios of 1.56, 2.25, and 4, respectively).

The dependence of the selectivity amplification constant $`\alpha `$ on the rate coefficient ratio $`\mathrm{\Omega }_2/\mathrm{\Omega }_1`$ is plotted in Fig. 4. The solid line depicts the estimate of Eq. (12). At the limit of a high rate coefficient ratio the selectivity is enhanced as $`N`$, implying that all the particles would dissociate into the preferred channel. We note that since the linearized model underestimates the condensate depletion time, the numerical values of $`\alpha `$ are consistently lower than the theoretical estimate of Eq. (12).

While the five-mode model is sufficient to produce Bose-enhanced gain, it neglects various spatial effects. In particular, channel losses due to the relative translational motion of the products are not accounted for under the assumption that energy conservation forces the reaction products to be at rest. In fact, by imparting additional momentum via the laser fields, such losses could be used to increase selectivity by setting a threshold for amplification. The population in each rearrangement channel would only be stimulated if the respective pump rate $`\sqrt{N}\mathrm{\Omega }_j`$ is faster than the channel depletion. Thus, if the motion in the $`AB+C`$ channel is faster than in the $`A+BC`$ channel it will be possible to tune the laser intensity so that only one pathway is Bose stimulated. In particular, if $`A`$ and $`B`$ are atoms of similar mass $`m_Am_B`$, and $`C`$ is a much lighter atom $`m_Cm_A,m_B`$, the $`A+BC`$ channel will be enhanced at much lower pump intensities than the $`AB+C`$ channel. Consequently, selectivity could be attained even for $`\mathrm{\Omega }_1=\mathrm{\Omega }_2`$, an issue currently under active investigation.

By drawing an analogy from the phenomena of BEC superradiance  to multi-channel reactive scattering in molecular condensates, we have shown that the selectivity of reactive processes carried out in a quantum-degenerate gas, can be greatly amplified by collective enhancement effects. This selectivity is the result of interplay between bose enhancement, which leads to exponential growth, and competition between modes for a finite number of initial molecules. The effect was demonstrated using the simple model system of two-channel photodissociation, where due to a small difference between the two channel rate coefficients, population was shown to be directed by the $`N`$-body stimulated process into the favored channel, approaching a nearly complete selectivity in the limit of large $`N`$. In the quantum degenerate regime, it is therefore only necessary to exert slight control over the reaction rates in order to achieve near-total control over the reaction output. Future work would address the question of complementary Pauli-blocking effects when a molecular BEC is dissociated into fermionic constituents.

We are grateful to James R. Anglin for valuable discussions. This work was supported by the National Science Foundation through a grant from the Institute for Theoretical Atomic and Molecular Physics at Harvard University and the Smithsonian Astrophysical Observatory.