text:
# Entanglement and spin squeezing of Bose condensed atoms

## I introduction

Entanglement lies at the heart of the difference between the quantum and classical multi-particle world. It is the phenomenon that enables quantum information and quantum computing with many qubits. Recently, several interesting developments have occurred in studies of massively entangled atomic states. Based on the proposals of Sorensen and Molmer Sørensen and Mølmer (1999a), a controlled, entangled state of $`4`$ions was successfully created by the NIST ion trap group Sackett et al. (2000). Zeilinger and coworkers prepared three entangled photon or GHZ states by selecting from two beams of entangled photon pairs Bouwmeester et al. (1999). Entanglement between two atoms and a microwave photon were also detected in a ”step-by-step” process Rauschenbeutel et al. (2000).

Of these and other related developments, the idea of Molmer and Sorensen Sørensen and Mølmer (1999a) is especially interesting. In their scheme a direct coupling to the multiparticle, entangled final state was made possible through a virtual, intermediate state which was a common (quantum) mode of the motion of all the ions. Similar type interactions were also obtained by Milburn milburn99 . Both proposals allow for creation of massive entangled states by unitary evolution, starting from certain pure initial states.

A Bose-Einstein condensate is a convenient source of atoms in initially pure and separable states. Zoller and coworkers suggested creating massive entangled spin squeezed states from a two component condensate using the inherent atom-atom interactions Sørensen et al. (2001) and further investigated schemes to dynamically create many particle entangled states of a two component BEC Micheli et al. . Spin exchange interaction of a spinor condensate was also proposed as a candidate for creating entangled pairs of atoms Pu and Meystre (2000); Duan et al. (2000).

Most of these proposals work in the two mode approximation where one motional state is assumed for each spinor component of condensed atoms. Recently, Sorensen demonstrated the validity of this method within Bogoliubov theory Sørensen (2002). Similar studies were previously performed, extensively, for condensate atoms in a double well potential Javanainen (1999); Zapata et al. (1998); Spekkens and Sipe (1999); Imamoglu et al. (1997); Milburn et al. (1997); Ho and Yip (2000).

Recently, two of us proposed a new type of atom-atom coupling that achieves optimal spin squeezing in two mode Bose condensed atoms Helmerson and You (2001). In this study we provide further details of our proposal and compare its potential for both the creation of massive entangled states Sørensen and Mølmer (1999a) and improved spin squeezing Sørensen et al. (2001); Pu and Meystre (2000); Duan et al. (2000); Kitagawa and Ueda (1993). This paper is organized as follows. First, we briefly review the Raman process with photon momentum transfer for a single three-level $`\mathrm{\Lambda }`$-type atom. We then extend the analysis to the case of two atoms, using intermediate molecular states. The resulting effective interaction is then compared with models in several other recent studies of creating entanglement with Bose condensed atoms. We then provide an explicit calculation of the time reversed dynamics to confirm the generation of maximally entangled states without requiring single atom resolved measurements Micheli et al. ; Vogels et al. . We conclude with a careful analysis of the limitations and advantages of our model.

## II Single atom, two-photon process through an atomic intermediate state

We consider a three level $`\mathrm{\Lambda }`$-type atom described by the Hamiltonian

$`H={\displaystyle \underset{\mu =g,e,g^{}}{}}\left[{\displaystyle \frac{p^2}{2M}}+V_\mu (\stackrel{}{r})+\mathrm{}\omega _{\mu g}\right]|\mu \mu |+H_{\mathrm{int}},`$ (1)
where the dipole interaction term is

$`H_{\mathrm{int}}`$ $`=`$ $`\stackrel{}{d}_{eg}\stackrel{}{E}_1(\stackrel{}{r})e^{i\stackrel{}{k}_1\stackrel{}{r}i\omega _1t}|ge|+h.c.`$ (2)
$`\stackrel{}{d}_{eg^{}}\stackrel{}{E}_2(\stackrel{}{r})e^{i\stackrel{}{k}_2\stackrel{}{r}i\omega _2t}|g^{}e|+h.c.`$
for the general Raman process from initial state $`|g`$ to final state $`|g^{}`$ through an intermediate state $`|e`$. $`V_\mu (\stackrel{}{r})`$ denotes the, possibly different, trapping potentials for the different electronic states $`|\mu `$ with internal energy $`\omega _{\mu g}=\omega _{\mu \mu }\omega _{gg}`$. $`\stackrel{}{E}_j(\stackrel{}{r})e^{i\stackrel{}{k}_j\stackrel{}{r}i\omega _jt}`$ denotes the running wave amplitude of the laser field.

We define $`\mathrm{\Delta }=\omega _1\omega _{eg}`$ and $`\delta _j=\omega _j\omega _1`$. Without loss of generality, assuming $`\stackrel{}{d}_{eg}=\stackrel{}{d}_{eg^{}}=\stackrel{}{d}`$, $`\mathrm{\Omega }_j=2\stackrel{}{d}\stackrel{}{E}_j/\mathrm{}`$ is the Rabi frequency for the corresponding dipole connected transition. In the limit of a resonant two photon process with large detuning from the intermediate atomic excited state $`|e`$, we adiabatically eliminate $`|e`$ to obtain an effective Hamiltonian

$`_{\mathrm{eff}}={\displaystyle \underset{\mu =g,g^{}}{}}\left[{\displaystyle \frac{P^2}{2M}}+V_\mu (\stackrel{}{r})\right]|\mu \mu |+V_{B/R},`$ (3)
with,

$`V_B`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{4\mathrm{\Delta }}}{\displaystyle \underset{j}{}}\mathrm{\Omega }_je^{i\stackrel{}{k}_j\stackrel{}{r}i\delta _jt}{\displaystyle \underset{j}{}}\mathrm{\Omega }_j^{}e^{i\stackrel{}{k}_j\stackrel{}{r}+i\delta _jt}|gg|,`$ (4)
$`V_R`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{4\mathrm{\Delta }}}\mathrm{\Omega }_1e^{i\stackrel{}{k}_1\stackrel{}{r}i\omega _1t}\mathrm{\Omega }_2^{}e^{i\stackrel{}{k}_2\stackrel{}{r}+i\omega _2t}|gg^{}|+h.c.,`$ (5)
where distinguish the two types Raman process according to whether the final internal state $`|g^{}`$ is the same as or different from the initial state $`|g`$, respectively. When $`|g^{}=|g`$, we identify the general Raman process as Bragg diffraction Kozuma et al. (1999).

In momentum space, the above two photon process involves the simultaneous absorption of a photon $`(\stackrel{}{k}_1,\omega _1)`$ by an atom in state $`|g;\stackrel{}{p}`$ and stimulated transition to state $`|g^{};\stackrel{}{p}\mathrm{}\stackrel{}{k}_2+\mathrm{}\stackrel{}{k}_1`$ with emission of a photon $`(\stackrel{}{k}_2,\omega _2)`$. For the case of Bragg diffraction ($`g^{}=g`$) with two counter-propagating waves

$`\stackrel{}{k}_1`$ $`=`$ $`\stackrel{}{k},`$
$`\stackrel{}{k}_2`$ $``$ $`\stackrel{}{k},`$
$`\stackrel{}{K}`$ $`=`$ $`\stackrel{}{k}_1\stackrel{}{k}_22\stackrel{}{k},`$ (6)
we obtain

$`V_B(\stackrel{}{r},t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{4\mathrm{\Delta }}}\left[\right|\mathrm{\Omega }_1|^2+|\mathrm{\Omega }_2|^2`$ (7)
$`+\mathrm{\Omega }_1\mathrm{\Omega }_2^{}e^{i\delta _2t}e^{i\stackrel{}{K}\stackrel{}{r}}+h.c.]|gg|.`$
The first two terms are AC Stark shifts that can be neglected. For an atom initially in state $`|g,\stackrel{}{p}_g`$, the last two terms of Eq. (7) couples to momentum states $`\stackrel{}{p}_g\pm 2\mathrm{}\stackrel{}{k}`$, respectively. By appropriately choosing $`\delta _2`$ one can selectively enhance coupling to only one momentum state or side mode. The resonance condition is defined by energy conservation in the two photon process

$`\mathrm{}\omega _1+{\displaystyle \frac{\stackrel{}{p}_g^2}{2M}}\mathrm{}\omega _2+{\displaystyle \frac{\left(\stackrel{}{p}_g+\mathrm{}\stackrel{}{K}\right)^2}{2M}},`$ (8)
which gives

$`\mathrm{}\delta _2={\displaystyle \frac{\stackrel{}{p}_g^2}{2M}}{\displaystyle \frac{\left(\stackrel{}{p}_g+\mathrm{}\stackrel{}{K}\right)^2}{2M}}{\displaystyle \frac{\mathrm{}^2\stackrel{}{K}^2}{2M}},`$ (9)
for an atom with $`\stackrel{}{p}_g0`$. Hence under conditions resonant for Bragg diffraction, the effective coupling (7) creates superposition states like,

$`|\psi =\alpha |g|\stackrel{}{p}_g+\beta |g|\stackrel{}{p}_g\pm 2\mathrm{}\stackrel{}{k}.`$ (10)
For Bose-condensed atoms, the above discussion still applies to any of the identical single atoms. When a condensate of N atoms is involved (initial motional state $`|p_g0`$), one simply creates a condensate in state (10), which is not an entangled state because it is simply putting N atoms into the same state (10), i.e.

$`\left(a_{|\psi }^{}\right)^N|\mathrm{vac}\psi (\stackrel{}{r}_1)\psi (\stackrel{}{r}_2)\mathrm{}\psi (\stackrel{}{r}_N)=|\psi ^N.`$ (11)
In practice this can only be done if atoms are noninteracting, or that the effect of interactions is small during the time of atomic Bragg diffraction. Bragg diffraction of condensate atoms was demonstrated first by the NIST-Gaithersburg group Kozuma et al. (1999). Higher motional states $`\stackrel{}{K}=\pm 2m\stackrel{}{k}`$ ($`m`$ integer) was also obtained through, higher order, $`2m`$-photon couplings.

For the Raman process with two co-propagating waves $`\stackrel{}{k}_1\stackrel{}{k}_2`$ involving nearly degenerate ground states $`|g`$ and $`|g^{}`$, we have

$`V_R={\displaystyle \frac{\mathrm{}}{4\mathrm{\Delta }}}\mathrm{\Omega }_1\mathrm{\Omega }_2^{}e^{i\delta _2t}|gg^{}|+h.c.,`$ (12)
where AC Stark shifts from second order processes involving the same laser fields have been neglected. The Raman resonance condition is then simply $`\omega _1\omega _2=\omega _{g^{}g}`$. Note that co-propagating Raman results in an effective coupling (12) which is constant, independent of the position $`\stackrel{}{r}`$. The two photon process then is essentially insensitive to motional effects $`(\stackrel{}{k}_2\stackrel{}{k}_1)`$.

## III Two-atom, two-photon process through an intermediate molecular state

We now consider a model involving two $`\mathrm{\Lambda }`$-type atoms whose initial and final states are described by the same non-interacting atomic states. The intermediate excited state, however, is now chosen to be a bound, molecular, excited state, similar to those utilized in recent photo-association experiments Wynar et al. (2000). Neglecting configurations not directly involved in the two photon process, the two atom Hamiltonian can be written as

$``$ $`={\displaystyle \underset{\mu =g,g^{}}{}}\left[{\displaystyle \underset{i=1,2}{}}\left({\displaystyle \frac{\stackrel{}{p}_i^2}{2M}}+V_\mu (\stackrel{}{r}_i)\right)+2\mathrm{}\omega _{\mu g}\right]|\mu ,\mu \mu ,\mu |+{\displaystyle \underset{\mu =g,g^{}}{}}U_{\mu \mu }(\left|\stackrel{}{r}_1\stackrel{}{r}_2\right|)|\mu ,\mu \mu ,\mu |`$ (13)
$`+\left[{\displaystyle \underset{i=1,2}{}}\left({\displaystyle \frac{\stackrel{}{p}_i^2}{2M}}+V_{g^{}g}(\stackrel{}{r}_i)\right)+U_{g^{}g}(\left|\stackrel{}{r}_1\stackrel{}{r}_2\right|)+\mathrm{}\omega _{g^{}g}\right]|g^{},gg^{},g|`$
$`+\left[{\displaystyle \underset{i=1,2}{}}\left({\displaystyle \frac{\stackrel{}{p}_i^2}{2M}}+V_{e_b}(\stackrel{}{r}_i)\right)+U_{e_b}(\left|\stackrel{}{r}_1\stackrel{}{r}_2\right|)+\mathrm{}\omega _{eg}\right]|e_be_b|+H_{\mathrm{int}},`$
where states $`|\mu ,\mu `$ for $`\mu =g,g^{}`$ denote symmetrized electronic states of two atoms, each in $`|\mu `$. $`|g^{},g`$ denote the symmetrized state with one in $`|g`$ and the other in $`|g^{}`$. $`V_{\mu =g,g^{},g^{}g,e_b}`$ denote external trapping potentials in state $`|\mu `$, $`U_{\mu \mu }`$ and $`U_{g^{}g}`$ are atom-atom interactions which at large inter-atomic distances are described by the usual van der Waals type terms, and at short range modified by coulomb effects. The intermediate molecular state family, denoted by $`|e_b`$, can be accessed through the direct dipole coupling term $`H_{\mathrm{int}}`$. $`U_{e_b}`$ is the Born-Oppenheimer molecular (of atom 1 and 2) binding potential for internal state manifold $`|e_b`$, which contain bound levels to be used as intermediate states.

Introducing center of mass $`\stackrel{}{R}=(\stackrel{}{r}_1+\stackrel{}{r}_2)/2`$ and relative coordinate $`\stackrel{}{r}=\stackrel{}{r}_2\stackrel{}{r}_1`$ for the two atoms, we can express

$`\stackrel{}{r}_1=\stackrel{}{R}{\displaystyle \frac{\stackrel{}{r}}{2}},`$ $`\stackrel{}{r}_2=\stackrel{}{R}+{\displaystyle \frac{\stackrel{}{r}}{2}},`$ (14)
and

$`{\displaystyle \frac{\stackrel{}{p}_1^2}{2M}}+{\displaystyle \frac{\stackrel{}{p}_2^2}{2M}}={\displaystyle \frac{\stackrel{}{P}^2}{2(2M)}}+{\displaystyle \frac{\stackrel{}{p}^2}{2(M/2)}}.`$ (15)
When $`|e_b`$ is asymptotically connected with $`|e,\mu =\left(|e|\mu +|\mu |e\right)/\sqrt{2}`$ (for $`\mu =g,g^{}`$), we can express the dipole coupling as

$`H_{\mathrm{int}}={\displaystyle \underset{i=1,2}{}}\stackrel{}{d}_i\stackrel{}{E}_1e^{i\stackrel{}{k}_1\stackrel{}{r}_ii\omega _1t}\left(|ge|\right)_i+h.c.`$
$`{\displaystyle \underset{i=1,2}{}}\stackrel{}{d}_i\stackrel{}{E}_2e^{i\stackrel{}{k}_2\stackrel{}{r}_ii\omega _2t}\left(|g^{}e|\right)_i+h.c.`$
$`={\displaystyle \frac{\mathrm{}\mathrm{\Omega }_1}{\sqrt{2}}}\mathrm{cos}(\stackrel{}{k}_1\stackrel{}{r}/2)e^{i\stackrel{}{k}_1\stackrel{}{R}i\omega _1t}|g,ge_b|+h.c.`$
$`{\displaystyle \frac{\mathrm{}\mathrm{\Omega }_2}{\sqrt{2}}}\mathrm{cos}(\stackrel{}{k}_2\stackrel{}{r}/2)e^{i\stackrel{}{k}_2\stackrel{}{R}i\omega _2t}|g^{},g^{}e_b|+h.c.,`$ (16)
where we have assumed $`\stackrel{}{d}_i=\stackrel{}{d}_i^{}=\stackrel{}{d}`$ for the $`|\mu |e`$ transition of atom $`i`$. We note that the electronic dipole coupling $`\mathrm{\Omega }_j`$ will vary considerably with $`\stackrel{}{r}`$ for small values of $`r`$, due to the relatively short range of the molecular interactions.

To express the intermediate state in terms of the particular molecular resonance or bound state $`|e_b`$, we now take a look at the eigen-structure within the $`|e_b`$ family. We will essentially base our discussion on some kind of molecular state asymptotically connected to the limit of one ground $`|g`$ and one excited atom $`|e`$, just as the $`O_g^{}`$ state, extensively discussed by D. Heinzen and others in their photo association work Boesten et al. (1999). The recent experiment by Heinzen and coworkers on the production of ground state molecules from an atomic condensate by a two-photon Raman process Wynar et al. (2000) provides additional motivation to explore our ideas experimentally. The photo-association process Boesten et al. (1999); Bohn and Julienne (1996); Kostrun et al. (2000) used by Heinzen et al. relies on the transition strength of going from a (two atom) ’free’ to a (molecular) ‘bound’ state via an intermediate, excited (molecular) ‘bound’ state. What we desire, on the other hand, is a transition from a (two atom) ‘free’ to ’free’ (relative motion) state via an intermediate, excited (molecular) ’bound’ state. More detailed discussions about feasibility of such transitions can be given provided all molecular potential curves are available.

The question of whether this is possible or not does not seem to depend on the sample density (in the weakly interacting limit), but will depend, as we shall show below, to a large degree on the trap strength and the excited bound state structure. Assuming all trapping potentials $`V_g`$, $`V_g^{}`$, and $`V_{e_b}`$ to be harmonic, the separation of center-of-mass and relative coordinates for the trapping potentials is

$`{\displaystyle \frac{1}{2}}M{\displaystyle \underset{j=1,2,3}{}}\nu _j^2x_{1j}^2+{\displaystyle \frac{1}{2}}M{\displaystyle \underset{j=1,2,3}{}}\nu _j^2x_{2j}^2`$ (17)
$`=`$ $`{\displaystyle \frac{1}{2}}\left(2M\right){\displaystyle \underset{j=1,2,3}{}}\nu _j^2X_j^2+{\displaystyle \frac{1}{2}}\left(M/2\right){\displaystyle \underset{j=1,2,3}{}}\nu _j^2x_j^2`$
$`=`$ $`V_{tR}+V_{tr},`$
where $`x_{ij}=\widehat{e}_j\stackrel{}{r}_i`$, $`X_{ij}=\widehat{e}_j\stackrel{}{R}`$, and $`x_j=\widehat{e}_j\stackrel{}{r}`$ with $`\widehat{e}_j`$ the unit vector in the $`j`$ direction.

When short range interactions $`U_{\mu \mu }(\stackrel{}{r})`$ are approximated by their optical contact forms, Wilkens and coworkers provide analytic solutions for two interacting atoms inside a harmonic trap Busch et al. (1998). In principle, from

$`\left[{\displaystyle \frac{p^2}{2(M/2)}}+V_{tr}^\mu (\stackrel{}{r})+U_{\mu \mu }(\stackrel{}{r})\right]\varphi _n^\mu (\stackrel{}{r})=E_n^\mu \varphi _n^\mu (\stackrel{}{r}),`$ (18)
we can find all bound states $`\varphi _n^\mu (\stackrel{}{r})`$ due to the external trapping potential $`U_{\mu \mu }(\stackrel{}{r})`$. We note that if the range of the inter-atomic interaction $`r`$ size of the ground state of the harmonic trap, then these bound states should resemble the standard low energy scattering solutions.

Similarly, the bound state of the intermediate, excited molecular level satisfies

$`\left[{\displaystyle \frac{p^2}{2(M/2)}}+V_{tr}^{e_b}(\stackrel{}{r})+U_{e_b}(\stackrel{}{r})\right]\varphi _m^{e_b}(\stackrel{}{r})=E_m^{e_b}\varphi _m^{e_b}(\stackrel{}{r}),`$ (19)
where the energy $`E_m^{e_b}`$ is measured from the asymptotic energy $`\mathrm{}\omega _{eg}`$. The index $`m`$ represents the rotational and vibrational quantum numbers. For all bound states with $`E_m^{e_b}<0`$, the effect of external trapping potential $`V_{tr}^{e_b}`$ is much less than for the free atoms, since the corresponding bound state wave function is localized inside the molecular potential well $`U_{e_b}`$.

Similarly, complete bound states $`|l`$ can also be introduced for $`|g^{},g`$ if needed. They satisfy

$`\left[{\displaystyle \frac{p^2}{2(M/2)}}+V_{tr}^{g^{}g}(\stackrel{}{r})+U_{g^{}g}(\stackrel{}{r})\right]\psi _l(\stackrel{}{r})=E_l^{g^{}g}\psi _l(\stackrel{}{r}).`$ (20)
In the $`(\stackrel{}{R},\stackrel{}{r})`$ basis, keeping only one intermediate bound state $`m_b`$, we can then write our previous Hamiltonian, Eq. (13), as

$``$ $`=`$ $`{\displaystyle \underset{n_\mu }{}}{\displaystyle \underset{\mu =g,g^{}}{}}\left[{\displaystyle \frac{\stackrel{}{P}^2}{2(2M)}}+V_{tR}^\mu (\stackrel{}{R})+2\mathrm{}\omega _{\mu g}+E_n^\mu \right]|\mu ,\mu ;n_\mu \mu ,\mu ;n_\mu |`$ (21)
$`+{\displaystyle \underset{l}{}}\left[{\displaystyle \frac{\stackrel{}{P}^2}{2(2M)}}+V_{tR}^{g^{}g}(\stackrel{}{R})+\mathrm{}\omega _{g^{}g}+E_l^{g^{}g}\right]|g^{},g;lg^{},g;l|`$
$`+\left[{\displaystyle \frac{\stackrel{}{P}^2}{2(2M)}}+V_{tR}^\mu (\stackrel{}{R})+\mathrm{}\omega _{eg}+E_{m_b}^{e_b}\right]|e_b;m_be_b;m_b|+H_{\mathrm{int}},`$
and the interaction term Eq. (III) now becomes

$`H_{\mathrm{int}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}\mathrm{\Omega }_1}{\sqrt{2}}}{\displaystyle \underset{n_g}{}}\eta _{n_gm_b}e^{i\stackrel{}{k}_1\stackrel{}{R}i\omega _1t}|g,g;n_ge_b;m_b|+h.c.`$ (22)
$`{\displaystyle \frac{\mathrm{}\mathrm{\Omega }_2}{\sqrt{2}}}{\displaystyle \underset{n_g^{}}{}}\eta _{n_g^{}m_b}e^{i\stackrel{}{k}_2\stackrel{}{R}/2i\omega _2t}|g^{},g^{};n_g^{}e_b;m_b|+h.c.,`$
where

$`\eta _{n_\mu m_b}`$ $`=`$ $`{\displaystyle \frac{1}{|\stackrel{}{d}|}}n_\mu \left|d(\stackrel{}{r})\mathrm{cos}(\stackrel{}{k}\stackrel{}{r}/2)\right|m_b`$ (23)
$``$ $`{\displaystyle \frac{1}{d}}{\displaystyle 𝑑\stackrel{}{r}\varphi _n^\mu (\stackrel{}{r})d(\stackrel{}{r})\mathrm{cos}(\stackrel{}{k}\stackrel{}{r}/2)\varphi _{m_b}^{e_b}(\stackrel{}{r})},`$
which is a measure of the free-bound transition amplitude. we note that in contrast to other expressions for the free-bound transition amplitude Boesten et al. (1999); Bohn and Julienne (1996); Kostrun et al. (2000), the motional wave functions $`|n_\mu `$ and $`|m_b`$ are now normalized in the trap rather than energy or momentum normalized. Although the modification of the inner part of the trap potential for $`r100`$ ($`a_0`$) due to the inter-atomic potentials is to be included, these values can be determined for selected molecular states from the results of photo-association experiments. They can also be computed directly if accurate potentials are available. Detailed discussions are given by Verhaar Boesten et al. (1999), Julienne Bohn and Julienne (1996) and Javaninen Kostrun et al. (2000).

In general, the coupling coefficient will vary depending on intermediate state. It is therefore important to pick the largest one of them. We then perform an adiabatic elimination of the excited state by appropriate choice of laser detunings, and obtain the effective Hamiltonian,

$``$ $`=`$ $`{\displaystyle \underset{n_\mu }{}}{\displaystyle \underset{\mu =g,g^{}}{}}\left[{\displaystyle \frac{\stackrel{}{P}^2}{2(2M)}}+V_{tR}^\mu (\stackrel{}{R})+2\mathrm{}\omega _{\mu g}+E_n^\mu \right]|\mu ,\mu ;n_\mu \mu ,\mu ;n_\mu |`$ (24)
$`+{\displaystyle \underset{l}{}}\left[{\displaystyle \frac{\stackrel{}{P}^2}{2(2M)}}+V_{tR}^{g^{}g}(\stackrel{}{R})+\mathrm{}\omega _{g^{}g}+E_l^{g^{}g}\right]|g^{},g;lg^{},g;l|+V_{B/R},`$
with

$`V_B`$ $`={\displaystyle \frac{\mathrm{}}{4\mathrm{\Delta }}}2{\displaystyle \underset{n,l}{}}\left[{\displaystyle \underset{j}{}}\mathrm{\Omega }_j\eta _{nm_b}e^{i\stackrel{}{k}_j\stackrel{}{R}i\omega _jt}\right]\left[{\displaystyle \underset{j}{}}\mathrm{\Omega }_j\eta _{lm_b}e^{i\stackrel{}{k}_j\stackrel{}{R}i\omega _jt}\right]^{}|g,g;ng,g;l|,`$
$`V_R`$ $`={\displaystyle \frac{\mathrm{}}{4\mathrm{\Delta }}}2{\displaystyle \underset{n,n^{}}{}}\left[\mathrm{\Omega }_1\eta _{nm_b}e^{i\stackrel{}{k}_j\stackrel{}{R}i\omega _1t}\right]\left[\mathrm{\Omega }_2\eta _{n^{}m_b}e^{i\stackrel{}{k}_2\stackrel{}{R}i\omega _2t}\right]^{}|g^{},g^{};n^{}g,g;n|+h.c.,`$ (25)
where we have defined $`\mathrm{\Delta }=\omega _1\omega _2E_{n_b}^{e_b}/\mathrm{}`$. Note the $`(\sqrt{2})^2`$ enhancement due to excitation to an asymptotic, symmetric state.

More generally, one can include all motional states and still be able to find the effective coupling between the selected electronic states, to obtain similar results You (2001).

Our analysis now follows along the same lines as in the previous case of a single atom. For example, if we consider only the lowest motional state $`n_\mu =0`$ (a situation well approximated by Bose condensed atoms), then our Hamiltonian becomes

$`={\displaystyle \underset{\mu =g,g^{}}{}}\left[{\displaystyle \frac{\stackrel{}{P}^2}{2(2M)}}+V_{tR}^\mu (\stackrel{}{R})+2\mathrm{}\omega _{\mu g}+E_m^\mu \right]|\mu ,\mu ;0\mu ,\mu ;0|+V_{B/R},`$ (26)
with

$`V_B`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{2\mathrm{\Delta }}}\eta _{0m_b}\eta _{0^{}m_b}^{}\left[{\displaystyle \underset{j}{}}\mathrm{\Omega }_je^{i\stackrel{}{k}_j\stackrel{}{R}i\omega _jt}\right]\left[{\displaystyle \underset{j}{}}\mathrm{\Omega }_je^{i\stackrel{}{k}_j\stackrel{}{R}i\omega _jt}\right]^{}|g,g;0g,g;0^{}|,`$
$`V_R`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{2\mathrm{\Delta }}}\eta _{0m_b}\eta _{0^{}m_b}^{}\mathrm{\Omega }_1\mathrm{\Omega }_2^{}\left[e^{i\stackrel{}{k}_j\stackrel{}{R}i\omega _1t}\right]\left[e^{i\stackrel{}{k}_2\stackrel{}{R}i\omega _2t}\right]^{}|g^{},g^{};0g,g;0^{}|+h.c.`$ (27)
For Bragg diffraction by two counter-propagating waves, the above result reduces to

$`V_B(\stackrel{}{R},t)={\displaystyle \frac{\mathrm{}}{2\mathrm{\Delta }}}\left|\eta _{0m_b}\right|^2\left[\left|\mathrm{\Omega }_1\right|^2+\left|\mathrm{\Omega }_2\right|^2+\mathrm{\Omega }_1\mathrm{\Omega }_2^{}e^{i\delta _2t}e^{i\stackrel{}{K}\stackrel{}{R}}+\mathrm{\Omega }_1^{}\mathrm{\Omega }_2e^{i\delta _2t}e^{i\stackrel{}{K}\stackrel{}{R}}\right]|g,g;0g,g;0|.`$ (28)
Similar to the single atom case, Bragg diffraction involves the simultaneous absorption and stimulated emission of photons, i.e. the elementary Bragg process involves a total of two photons. In contrast, however, to the single atom case, a pair of atoms are now diffracted as illustrated in Fig. 1. Hence for Bose condensed atoms with $`p_i0`$, Bragg diffraction produces atoms in the state

$`|\psi _M`$ $`=`$ $`\alpha |g,g;0|p_10,p_20+\beta |g,g;0e^{\pm i\stackrel{}{K}\stackrel{}{R}}`$ (29)
$`=`$ $`\alpha |g,g|p_10,p_20`$
$`+\beta |g,g|\stackrel{}{p}_1=\mathrm{}\stackrel{}{K}/2,\stackrel{}{p}_2=\mathrm{}\stackrel{}{K}/2.`$
Not surprisingly, the momentum shift is only half the value of atomic Bragg diffraction, due to the fact that only two photons are involved for the pair of atoms. The resonance condition occurs at

$`\mathrm{}\omega _1+{\displaystyle \frac{P_i^2}{2(2M)}}+{\displaystyle \frac{p_i^2}{2(M/2)}}`$ $`=`$ $`\mathrm{}\omega _2+{\displaystyle \frac{P_f^2}{2(2M)}}+{\displaystyle \frac{p_f^2}{2(M/2)}},`$
$`p_i`$ $``$ $`p_f0,`$
$`\stackrel{}{P}_f`$ $`=`$ $`\stackrel{}{P}_i\pm \mathrm{}\stackrel{}{K},`$ (30)
which gives the resonance condition

$`\mathrm{}\delta _2={\displaystyle \frac{P_i^2}{2(2M)}}{\displaystyle \frac{(\stackrel{}{P}_i\pm \mathrm{}\stackrel{}{K})^2}{2(2M)}}{\displaystyle \frac{\mathrm{}^2\stackrel{}{K}^2}{2(2M)}},`$ (31)
i.e. half the atomic Bragg resonance (9).

It is worth mentioning that for a condensate with many atoms, the state obtained from diffraction via a two atom intermediate state will be complicated. It is generally not of the form $`(a_{|\psi _M}^{})^{N/2}|\mathrm{vac}`$.

For a nearly degenerate Raman process with two co-propagating waves, we obtain for atoms initially at rest,

$`V_R(\stackrel{}{R},t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{2}}{\displaystyle \frac{\mathrm{\Omega }_1\mathrm{\Omega }_2^{}}{\mathrm{\Delta }}}\eta _{0m_b}\eta _{0^{}m_b}^{}`$ (32)
$`e^{i\delta _2t}|g,g;0g^{},g^{};0^{}|+h.c.`$
The resonance condition now becomes $`\omega _1\omega _2=2\omega _{g^{}g}`$, i.e. the Raman process changes each atom’s state ($`gg^{}`$) and so each atom acquires the energy deficit $`\mathrm{}\omega _{g^{}g}`$. Note that there is no $`\stackrel{}{R}`$ dependence of the coupling coefficient in this case.

If the Raman process is arranged to allow for simultaneous scattering to two different final states $`|g^{}`$ and $`|g^{\prime \prime }`$, the coherent coupling becomes,

$`V_R^{[2]}(\stackrel{}{R},t)`$ $`=`$ $`{\displaystyle \frac{\mathrm{}\mathrm{\Omega }_R}{2}}e^{i\delta _2t}|g,g;0g^{},g^{};0^{}|+h.c.`$ (33)
$`+{\displaystyle \frac{\mathrm{}\mathrm{\Omega }_R^{}}{2}}e^{i\delta _2^{}t}|g,g;0g^{\prime \prime },g^{\prime \prime };0^{\prime \prime }|+h.c.`$
Such a coupling can be obtained among different Zeeman levels in the atomic ground state.

Finally, we note that our discussions above can also be applied to a pair of different species of atoms, e.g. using molecular intermediate state from Li-Cs dimer could create entangled pairs of Li and Cs. We also note that to be formally correct, the second order perturbation process also needs to be considered for the state $`|g^{},g`$, which should have some type of effective couplings with itself and as well as with the states $`|g,g`$ and $`|g^{},g^{}`$. we have neglected those couplings, assuming that resonance conditions for their survival in the interaction picture are not satisfied.

## IV Spin-spin interaction models and two mode approximations

In this section, we will be primarily concerned with process which produce entanglement of atomic spin states. We therefore consider a spin 1/2 system and adopt the notation $`|`$ for spin up and $`|`$ for spin down, with the following operators and identities:

$`\sigma _z`$ $`=`$ $`||||,`$
$`\sigma _+`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\sigma _x+i\sigma _y)=||,`$
$`\sigma _{}`$ $`=`$ $`{\displaystyle \frac{1}{2}}(\sigma _xi\sigma _y)=||,`$
$`\sigma _+^2`$ $`=`$ $`\sigma _{}^2=0,`$
$`\sigma _+\sigma _{}+\sigma _{}\sigma _+`$ $`=`$ $`1,`$
$`\sigma _x^2`$ $`=`$ $`\sigma _y^2=\sigma _z^2=1.`$ (34)
The Bragg/Raman coupling as discussed here is of the form

$`{\displaystyle \frac{\mathrm{}\mathrm{\Omega }_R}{2}}[(||)_1(||)_2+(||)_1(||)_2]`$
$`={\displaystyle \frac{\mathrm{}\mathrm{\Omega }_R}{2}}(\sigma _+^{(1)}\sigma _+^{(2)}+\sigma _{}^{(1)}\sigma _{}^{(2)}),`$ (35)
if we designate $`|g`$ and $`|g^{}`$ as $`|`$ and $`|`$, respectively. The Rabi frequency $`\mathrm{\Omega }_R`$ is the two-photon Rabi frequency for transitions between $`|`$ and $`|`$ via the intermediate molecular state. Hence it contains the single photon Rabi frequencies $`\mathrm{\Omega }_1`$ and $`\mathrm{\Omega }_2`$, as well as the free-bound transition amplitude, $`\eta _{n_\mu m_b}`$. This is a new type of coupling not widely discussed before. Its ability, however, to generate entanglement should be obvious as two atoms perform conditional evolution at all times, similar to the coupling for photon down conversion. Its prospects for creating massive GHZ states will now be studied.

We first compare our coupling scheme with other relevant models. In the original scheme of Sorensen and Molmer Sørensen and Mølmer (1999a), the two atom coupling takes the following form,

$`{\displaystyle \frac{\mathrm{}\mathrm{\Omega }_R}{2}}\sigma _x^{(1)}\sigma _x^{(2)}`$
$`={\displaystyle \frac{\mathrm{}\mathrm{\Omega }_R}{2}}\left(\right||+|\left|\right)_1\left(\right||+|\left|\right)_2,`$ (36)
a form different from ours. It is convenient to analyze such models for the case of many atoms in terms of collective spin operator $`J_{\mu =x,y,z}=\frac{1}{2}_i\sigma _\mu ^{(i)}`$. One can then show that the Hamiltonian for the Sorensen and Molmer scheme becomes

$`{\displaystyle \underset{i<j}{}}{\displaystyle \frac{\mathrm{}\mathrm{\Omega }_R}{2}}\sigma _x^{(i)}\sigma _x^{(j)}\mathrm{}\mathrm{\Omega }_R\left(J_x^2{\displaystyle \frac{N}{4}}\right),`$ (37)
while our model gives

$`{\displaystyle \underset{i<j}{}}{\displaystyle \frac{\mathrm{}\mathrm{\Omega }_R}{2}}\left(\sigma _+^{(i)}\sigma _+^{(j)}+\sigma _{}^{(i)}\sigma _{}^{(j)}\right)`$
$`={\displaystyle \frac{\mathrm{}\mathrm{\Omega }_R}{4}}\left(J_+^2+J_{}^2\right)`$
$`={\displaystyle \frac{\mathrm{}\mathrm{\Omega }_R}{2}}\left(J_x^2J_y^2\right).`$ (38)
Recently, a many body, two mode coupling was proposed by Zoller and coworker Sørensen et al. (2001). They considered an interacting, two component (i.e. $`|`$ and $`|`$) condensate whose many body Hamiltonian becomes

$`H`$ $`=`$ $`{\displaystyle \underset{j=,}{}}ϵ_{0,j}{\displaystyle 𝑑\stackrel{}{r}\left|\varphi _0(\stackrel{}{r})\right|^2a_j^{}a_j}`$ (39)
$`+{\displaystyle \frac{1}{2}}{\displaystyle \underset{j=,}{}}U_{jj}{\displaystyle 𝑑\stackrel{}{r}\left|\varphi _0(\stackrel{}{r})\right|^4a_j^{}a_j^{}a_ja_j}`$
$`+U_{}{\displaystyle 𝑑\stackrel{}{r}\left|\varphi _0(\stackrel{}{r})\right|^4a_{}^{}a_{}^{}a_{}a_{}}.`$
under the two mode approximation, where each component has the same spatial mode,

$`\mathrm{\Psi }_j(\stackrel{}{r})=\varphi _0(\stackrel{}{r})a_j.`$ (40)
Introducing the collective spin operators

$`J_x`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(a_0^{}a_1+a_1^{}a_0\right),`$
$`J_y`$ $`=`$ $`{\displaystyle \frac{i}{2}}\left(a_0^{}a_1a_1^{}a_0\right),`$
$`J_z`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left(a_1^{}a_1a_0^{}a_0\right),`$ (41)
and taking the index $`0`$ for $`|`$ and $`1`$ for $`|`$ Spekkens and Sipe (1999); Milburn et al. (1997); Elgaroy and Pethick (1999); Yurke (1986), one can show that $`J_{x,y,z}=\frac{1}{2}_i\sigma _{x,y,z}^{(i)}`$ as was used earlier. We also note that $`N=a_1^{}a_1+a_0^{}a_0`$ and the Casimir relation

$`J_x^2+J_y^2+J_z^2={\displaystyle \frac{N}{2}}({\displaystyle \frac{N}{2}}+1).`$ (42)
The interaction term Eq. (39) then becomes

$`{\displaystyle \frac{1}{2}}J_z^2\left(U_{00}+U_{11}2U_{01}\right){\displaystyle 𝑑\stackrel{}{r}\left|\varphi _0(\stackrel{}{r})\right|^4}.`$ (43)
In the discussions of the numerical solutions to follow, we will use the second quantized notation, and consider a pure initial state with a fixed, total number of atoms $`N`$. Specifically, we compare the three types of coupling,

$`V_M=\mathrm{}\mathrm{\Omega }_RJ_x^2,`$ (44)
considered by Sorensen and Molmer Sørensen and Mølmer (1999a), the spin squeezing model

$`V_S=\mathrm{}\mathrm{\Omega }_RJ_z^2,`$ (45)
recently considered by Zoller and coworkers Sørensen et al. (2001), and our proposed coupling

$`V_{B/R}={\displaystyle \frac{\mathrm{}\mathrm{\Omega }_R}{2}}\left(J_x^2J_y^2\right).`$ (46)
Given an initial state $`|\psi (0)=a_1^N|\mathrm{vac}/\sqrt{N!}`$, i.e. with all atoms condensed in state $`|`$, the time evolution operator for $`V_M`$ is analytically known and takes a simple form at $`\mathrm{\Omega }_Rt=\pi /2`$ for $`N`$ an even integer Sørensen and Mølmer (1999a),

$`U_M=e^{iJ_x^2(\mathrm{\Omega }_Rt)}=e^{i\pi N/2+i\pi /4}{\displaystyle \frac{1}{N!}}a_0^Na_1^N+e^{i\pi /4}.`$ (47)
It creates the maximally entangled N-GHZ state

$`|\mathrm{GHZ}_N`$ $``$ $`{\displaystyle \frac{1}{\sqrt{2}}}\left[|0_0|N_1+i^{N+1}|N_0|0_1\right].`$ (48)
In the context of Bose-Einstein condensation, state $`|\mathrm{GHZ}_N`$ is an example of the interesting fragmented condensate Nozieres and James (1982). The Sorensen and Molmer model (Eq. (44)) produces perfect GHZ-type states at selected times, while our coupling (Eq. (46)) doesn’t, in general, produce exact GHZ-type state. From numerical simulations, however, we find that our model, on average, produces more than 50% overlap with GHZ-type state (Eq. (48)) at selected times. Because of the geometric equivalence between $`J_x^2`$ and $`J_z^2`$, the spin squeezing model (Eq. (45)) also produces a perfect GHZ at selected times, except now the entanglement is between single particle states $`|\pm =(|\pm |)/\sqrt{2}`$ You  and the initial state has to become $`|\pm ^N`$.

For our model, we were not able to find a general analytic form of the time evolution operator, even at specific values of $`t`$. To numerically calculate the time evolution, we expand the wave function as

$`|\psi (t)={\displaystyle \underset{m=0}{\overset{N}{}}}c_m(t)|Nm_0|m_1,`$ (49)
where

$`|m_i={\displaystyle \frac{a_i^m}{\sqrt{m!}}}|\mathrm{vac},`$ (50)
with the initial conditions given by the $`c_m(0)`$s.

In Figs. 2 and 3 we show, for comparison, the generation of the maximum entangled GHZ state Eq. (48) as determined by model $`V_M`$ (Eq. (44)) and our “molecular diffraction” model (Eq. (46)). We see that while the Sorensen and Molmer model creates perfect massive GHZ states at selected times, our model can also creates significant overlap with the massive GHZ at selected times.

We can also compare the achievable spin squeezing between our model and that of $`V_S`$ (Eq. (45)), using the squeezing parameter

$`\xi ^2={\displaystyle \frac{N(\mathrm{\Delta }J_{\stackrel{}{n}_1})^2}{J_{\stackrel{}{n}_2}^2+J_{\stackrel{}{n}_3}^2}},`$ (51)
where $`\stackrel{}{n}_i,i=1,2,3`$ are mutually orthogonal unit vectors. Other discussions of spin squeezed states can be found in Kitagawa and Ueda (1993); Sørensen and Mølmer (1999b); Hald et al. (1999).

For the initial state

$`|\psi (0)=|+^N={\displaystyle \frac{1}{2^{N/2}}}{\displaystyle \frac{1}{\sqrt{N!}}}\left(a_0^{}+a_1^{}\right)^N|\mathrm{vac},`$ (52)
it has been shown Sørensen et al. (2001) that $`\xi ^2(t>0)<1`$ for some set of $`\stackrel{}{n}_i`$’s. This model of $`V_S`$ (Eq. (45)) is in fact the one-axis twisting model considered by Kitagawa and Ueda earlier Kitagawa and Ueda (1993). In this case the problem can be solved analytically, with the result

$`(\mathrm{\Delta }J_{\stackrel{}{n}_1})_{\mathrm{min}}^2N^{1/3}.`$ (53)
On the other hand, our model resembles the two-axis twisting model of Kitagawa and Ueda Kitagawa and Ueda (1993), and has to be solved numerically. In the limit of large $`N`$ and with the condensate initially in $`|^N`$, one can show that for an initial $`J_z=N/2`$ ($`\stackrel{}{J}`$-spin pointing along the positive z-direction),

$`(\mathrm{\Delta }J_{\stackrel{}{n}_1})_{\mathrm{min}}^2{\displaystyle \frac{1}{2}}.`$ (54)
The optimal squeezing in this case occurs with $`J_x+J_y`$. This result can be easily verified by taking a semiclassical approximation in the dynamical equation for $`J_x`$ and $`J_y`$. We find that the time scale of reaching maximum squeezing is $`1/(N\mathrm{\Omega }_R)`$ (see also Wineland et al. (1994)). For condensates containing $`10^6`$ atoms, even with a very weak coupling $`\mathrm{\Omega }_R1`$ (Hz), the maximum squeezing is reached within a micro-second.

In Figure 4 the squeezing parameters and the corresponding angle $`\theta _{\mathrm{min}}`$ are computed numerically for the model $`V_S`$ (Eq. (45)) with the initial state given by Eq. (52). Fig. 5 is a similar calculation for our model (Eq. (46)) starting with all atoms in state $`|`$. Our model achieves a better squeezing at an earlier time. In addition, the phase angle $`\theta _{\mathrm{min}}`$ is fixed at the value $`\pi /4`$. More detailed discussions of the spin squeezing properties of our model can be performed, including dissipation and finite system size Vardi and Anglin (2001); Law et al. (2000); André and Lukin (2002).

Finally we also note that our coupling scheme for simultaneous transitions to two different final states gives

$`V_{B/R}^{[2]}={\displaystyle \frac{\mathrm{}\mathrm{\Omega }_R}{4}}[(a_0^{}a_0^{}a_1a_1+a_0^{}a_0^{}a_1^{}a_1^{})+h.c.].`$ (55)
In the undepleted pump limit when dissipation exists, or when the coupling efficiency is small, this coupling can be approximated by

$`V_{B/R}^{[2]}\left(a_1^{}a_1^{}+a_1^{}^{}a_1^{}^{}\right).`$ (56)
This is similar to the recent ideas of Refs. Pu and Meystre (2000); Duan et al. (2000); Vogels et al.  which creates entangled EPR pairs every time an event occurs, just like in a parametric down conversion process Shih and Alley (1988); Hong et al. (1987).

## V Detection of the maximally entangled states

While it seems that there are promising schemes to generating massive entanglement of condensate atoms, the detection of such entanglement represents a significant challenge by itself as it is very difficult to perform individual atom resolved measurements Sackett et al. (2000); Wineland et al. (1994) in a condensate. Recently, it was suggested that performing a time reversed dynamic evolution can detect the coherence of the created N-GHZ state Micheli et al. ; Vogels et al. .

In this section, we detail the working mechanism of the proposed scheme by performing explicit calculations. We will use the $`V_M`$ model as an example. To simplify the algebra, we will use $`\tau =\mathrm{\Omega }_Rt`$, and denote the basis states according to

$`|n_1`$ $``$ $`|n_0,n_1=|n_0_0|n_1_1.`$ (57)
Assume we start with $`|^N`$, then at $`\tau =\pi /2`$, a maximally entangled state

$`|\psi (\tau ={\displaystyle \frac{\pi }{2}}){\displaystyle \frac{1}{\sqrt{2}}}(|N+i^{N+1}|0),`$ (58)
Immediately reversing the time dynamics ($`V_MV_M`$), it is easy to show that after another $`\tau =\pi /2`$, we recover the initial state $`|^N`$, i.e. all atoms to be detected are in state $`|`$. If instead, the created state at $`\tau =\pi /2`$ is a completely incoherent mixture of $`|0`$ and $`|N`$, given by the (N-body) density matrix

$`\rho (\tau ={\displaystyle \frac{\pi }{2}})={\displaystyle \frac{1}{2}}|00|+{\displaystyle \frac{1}{2}}|NN|,`$ (59)
then the time reversed dynamics will lead to observations of variable number of atoms in state $`|`$ as well.

For the general case, let’s assume the restricted form

$`\rho (\tau ={\displaystyle \frac{\pi }{2}})`$ $`=`$ $`p_0|00|+p_N|NN|`$ (60)
$`+\beta |0N|+\beta ^{}|N0|,`$
with $`0p_{0/N}1`$, $`p_0+p_N=1`$, and $`|\beta |^2p_0p_N`$. Time reversed evolution of $`V_M`$ for $`\tau `$ from $`\tau =\pi /2`$ will lead to

$`\rho (\tau +\pi /2)`$ $`=`$ $`U_M(\tau )\rho (\pi /2)U_M^{}(\tau )`$ (61)
$`=`$ $`e^{iJ_x^2\tau }\rho (\pi /2)e^{iJ_x^2\tau },`$
which can be explicitly evaluated with the use of rotation operation

$`e^{\pm iJ_x^2\tau }`$ $`=`$ $`e^{i\frac{\pi }{2}J_y}e^{\pm iJ_z^2\tau }e^{i\frac{\pi }{2}J_y},`$ (62)
and the disentangling relation Arecchi et al. (1972)

$`e^{i\frac{\pi }{2}J_y}`$ $`=`$ $`e^{J_+}e^{J_z\mathrm{ln}2}e^J_{}=e^J_{}e^{J_z\mathrm{ln}2}e^{J_+},`$ (63)
$`e^{i\frac{\pi }{2}J_y}`$ $`=`$ $`e^{J_+}e^{J_z\mathrm{ln}2}e^J_{}=e^J_{}e^{J_z\mathrm{ln}2}e^{J_+}.`$ (64)
Denote

$`A`$ $`=`$ $`N\left|e^{J_+}e^{iJ_z^2\tau }e^{J_+}\right|0,`$ (65)
$`B`$ $`=`$ $`N\left|e^{J_+}e^{iJ_z^2\tau }e^J_{}\right|N,`$ (66)
then we find

$`\rho _{NN}(\tau +\pi /2)=N|\rho (\tau +\pi /2)|N`$
$`={\displaystyle \frac{1}{2^{2N}}}[|A|^2p_0+|B|^2p_N+(AB^{}\beta +h.c.)]`$ (67)
with

$`A`$ $`=`$ $`{\displaystyle \underset{m=0}{\overset{N}{}}}{\displaystyle \frac{N!}{m!(Nm)!}}(1)^me^{i(m\frac{N}{2})^2\tau },`$ (68)
$`B`$ $`=`$ $`{\displaystyle \underset{m=0}{\overset{N}{}}}{\displaystyle \frac{N!}{m!(Nm)!}}e^{i(m\frac{N}{2})^2\tau }.`$ (69)
The following angular momentum algebra has been used,

$`J_+|m`$ $`=`$ $`\sqrt{(m+1)(Nm)}|m+1,`$
$`J_{}|m`$ $`=`$ $`\sqrt{m(Nm+1)}|m1,`$
$`J_z|m`$ $`=`$ $`\left(m{\displaystyle \frac{N}{2}}\right)|m,`$
$`m|J_+`$ $`=`$ $`(J_{}|m)^{}.`$ (70)
At $`\tau =\frac{\pi }{2}`$, we find

$`\rho _{NN}(\pi )={\displaystyle \frac{1}{2}}+(1)^{\frac{N}{2}}\mathrm{𝐈𝐦}(\beta )1,`$ (71)
because $`|\mathrm{𝐈𝐦}(\beta )|\left|\beta \right|\sqrt{p_0p_N}1/2`$, for any state. Thus we conclude that the probability of detecting atoms in state $`|`$ is nonzero for mixed states. In fact, for the completely incoherent mixture Eq. (59) with $`p_0=p_N=\frac{1}{2}`$ and $`\beta =0`$, $`\rho _{00}(\pi )=\rho _{NN}(\pi )=1/2`$.

It’s interesting to compare the average populations and fluctuations of spin up $`|`$ number operator ($`n_1=a_1^{}a_1`$) after the time reversed evolution for $`\tau =\frac{\pi }{2}`$ starting with the following different initial conditions,

* A pure state,
$`|\psi (\pi /2)=\frac{1}{\sqrt{2}}(|N+i^{N+1}|0)`$

$`n_1(\pi )`$ $`=`$ $`\mathrm{𝐭𝐫}(\rho n_1)=N,`$
$`n_1^2(\pi )`$ $`=`$ $`\mathrm{𝐭𝐫}(\rho n_1^2)=N^2,`$
$`\mathrm{\Delta }n_1(\pi )`$ $`=`$ $`\sqrt{n_1^2n_1^2}=0;`$ (72)
* An incoherent mixture,
$`\rho (\pi /2)=\frac{1}{2}(|NN|+|00|)`$

$`n_1(\pi )`$ $`=`$ $`\mathrm{𝐭𝐫}(\rho n_1)={\displaystyle \frac{N}{2}},`$
$`n_1^2(\pi )`$ $`=`$ $`\mathrm{𝐭𝐫}(\rho n_1^2)={\displaystyle \frac{N^2}{2}},`$
$`\mathrm{\Delta }n_1(\pi )`$ $`=`$ $`{\displaystyle \frac{N}{2}}.`$ (73)

Thus, these two cases become easily distinguishable. In Fig. 6, we compare the time dependent results of averaged population in state $`|`$ and its variance, from the time reversed evolution starting with the coherent GHZ state Eq (58) and the incoherent mixture Eq. (59). It is interesting to note that the average population of spin up $`|`$ number operator ($`n_1=a_1^{}a_1`$) for the completely incoherent mixture is in fact independent of $`\tau `$. We find that the width of the Gaussian like features near $`\mathrm{\Omega }_Rt\pi `$ is inversely proportional to $`/\sqrt{N}`$, but the plateau for the variance $`\mathrm{\Delta }n_1`$ is at around $`N/\sqrt{3}`$.

For factorizable N-atom states, whether a pure state $`|\psi (\pi /2)=\frac{1}{2^{N/2}}(|+|)^N`$ or a mixture $`\rho (\pi /2)=\frac{1}{2^N}(||+||)^N`$, the probability distribution after the time reversed evolution gives identical results that are different from the two multi-atom correlated states discussed earlier. In fact, the final results simply reflects the binomial distributions of $`n_1`$ in a equal superposition separable state and is given by

$`n_1(\pi )`$ $`=`$ $`\mathrm{𝐭𝐫}(\rho n_1)={\displaystyle \frac{N}{2}},`$
$`n_1^2(\pi )`$ $`=`$ $`\mathrm{𝐭𝐫}(\rho n_1^2)={\displaystyle \frac{N(N+1)}{4}},`$
$`\mathrm{\Delta }n_1(\pi )`$ $`=`$ $`{\displaystyle \frac{\sqrt{N}}{2}}.`$ (74)
In this particular example with $`V_M`$ as the interaction, the above results actually is independent of $`\tau `$. Thus the measurement of total single particle state population fluctuations after the time reversed evolution provides a strong indication for the coherence and correlation of the created entangled states.

## VI Limitations and advantages

There are several limitations to our model that needs to be discussed carefully. First, the molecular coupling needs to dominate over the atomic coupling; i.e., we need to achieve a molecular coupling

$`\mathrm{\Omega }_R^M`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Omega }_1\mathrm{\Omega }_2^{}}{\mathrm{\Delta }_M}}\eta _{0m_b}\eta _{0^{}m_b}^{}`$ (75)
$``$ $`\mathrm{\Omega }_R^A={\displaystyle \frac{\mathrm{\Omega }_1\mathrm{\Omega }_2^{}}{\mathrm{\Delta }_A}},`$
where $`\mathrm{\Delta }_{M/A}`$ are detunings from the molecular and atomic intermediate states, respectively. We can estimate analytically the free-bound transition amplitude using the interacting ground state wave function from Busch et al. (1998), and excited state wave function from the bound state solutions in a $`C_3/r^3`$ potential modified by an inner repulsive core Gao (1999). Most likely, however, the deepest molecular bound state with significant amplitude should be chosen in order to maximize the detuning from the atomic transition. In addition, we should have $`\mathrm{\Delta }_M\gamma _M`$ to minimize spontaneous emission which would lead to decoherence and loss of atoms.

A clear advantage of our scheme is the signature for diffraction via a molecular intermediate state. The diffracted atoms would move at half the speed of atoms that have undergone atomic Bragg diffraction. Even as the number of atoms is increased, our scheme achieves the same level of squeezing and the same high value of overlap with the massive GHZ state.

Our scheme is based on an engineered interaction that can be turned on and off like the scheme of Sorensen and Molmer (Eq. (44)). In particular, our scheme works for non-interacting atoms which can decrease the noise due to atom-atom interactions in a $`U(1)`$ symmetry breaking condensate state. We have also provided explicit calculations that demonstrate that the time reversed dynamics can be used to verify the coherence of the intended maximally entangled states.

Among all models discussed, there exists an even/odd number problem: The form of the final state depends on whether the total number of condensed atoms is even or odd number. If, however, the condensate resembles a coherent superposition of number states (grand canonical ensembles), we should average the results over the number distribution $`p(N)`$ as well. This will be explored in a future study.

This work benefited from our participation of the 2000 workshop on quantum degenerate gases at the Lorentz Center, University of Leiden. We thank H. Stoof for the hospitality. We thank Dr. M. S. Chapman and S. L. Rolston for helpful discussions. This work is supported by a grant from NSA, ARDA, and DARPA under ARO Contract No. DAAD19-01-1-0667, and by a grant from the NSF PHY-0113831.


format: text

text:
# Quenching of High 𝑝_⟂ Hadron Spectra by Hadronic Interactions in Heavy Ion Collisions at RHIC

## I Motivation and estimates

One of the major goals of the ongoing experiments at the Relativistic Heavy Ion Collider (RHIC) in Brookhaven National Laboratory is to find and probe for a temporarily occurring macroscopical state of deconfined quark–gluon matter. So far, especially hadronic abundancies have been experimentally thoroughly studied. Total hadron multiplicity measurements, however, do reflect mainly the chemical equilibration and possible hadronization processes of the late stages of the system, and do not deliver direct information on the early stage with maximum energy density and temperature, where a deconfined and highly excited state is generally expected to be formed. Still, one of the very interesting first results of the RHIC experiments is, that one has established a significant suppression of moderately high $`p_{}`$ hadrons produced in central $`A+A`$ collisions compared to rescaled peripheral collisions or rescaled (and extrapolated) $`p+p`$ collisions: The results for the lower center–of–mass energy of $`130\mathrm{AGeV}`$ PHENIX ; STAR  and also the preliminary data for $`200\mathrm{AGeV}`$ QM2002  show a suppression factor $`R(p_{})`$ of about $`1/5`$ for pions with $`p_{}5\mathrm{GeV}`$, and hence state a clear hint for (probably various) nuclear medium effects at work.

The most popular explanation for this phenomenon is the onset of the occurrence of so called ‘jet quenching’, anticipated in Wang . The idea is, that a high energetic parton moving through a dense coloured (and deconfined) medium will loose considerable energy due to collisions and induced gluon radiation, and its final fragmentation will give rise to particles with considerable lower energies Wang1 . Hence, measurements of jets seem to offer a direct access to probe the early stage when the deconfined matter is very dense. As a present theoretical consensus, the dominant mechanism for energy loss is thought to be the nonabelian radiation of gluons based on pQCD by the energetic partons BDMS ; Baier . As a result, the momenta of the jet partons are attenuated before hadronizing. One should internalize the word of caution, that this appealing picture of applying perturbative QCD calculations is strictly valid only for such high energetic jets with $`E>10\mathrm{}20\mathrm{GeV}`$ Baier . And indeed, it has been found that the suppression pattern of the hard jets (with $`E>10\mathrm{GeV}`$) can show some significant dependence on the initial properties of the deconfined state Kai .

A recent calculation based on the rather involved GLV formalism for a finite number of coloured collisions Levai  does indeed provide a good agreement with the data. Still, though, a phenomenological opacity parameter has to be adjusted Levai2 . In this study already various other effects (the Cronin effect as well as a slight modification of the gluon distribution due to shadowing) have been phenomenologically incorporated. These effects do partly counteract and compete with the pure parton energy loss and thus do make a detailed analysis already more delicate JMS02 ; VG02 . Contrary, the proposal was made that the observed spectra for central collisions show a significant $`m_{}`$ scaling, being a manifestation of a direct remnant of the initial gluons which were liberated from gluon saturated nuclear distribution functions SKMV01 . Here any possible later interactions of the jets with the surrounding medium which might or should alter the distributions are completly discarded.

It has also been raised very recently, that in an ideal strong quenching picture the $`p_{}`$ spectrum of the high momentum hadrons is simply given by early ‘surface emission’ of the outer regions, where roughly half of the there produced jets can escape into the vacuum without passing through any medium Muller . It has to be verified, whether such an idealised scenario of complete absorption of energetic partons in matter can actually be achieved. The proposal bears resemblance to a similar geometrical picture for the production of another hard hadronic probe, the $`J/\mathrm{\Psi }`$, being invoked much earlier for a possible explanation of the so called ‘anomalous’ suppression pattern observed in $`Pb+Pb`$ collisions at CERN–SPS blaizot .

Typically, it is assumed in all the above refered descriptions that the partons do exit the collision region as pointlike particles before finally fragmenting into a ‘jet’ of hadrons. We now will point out that this prejudice should be of major concern, being invalid for present settings at RHIC. The potential magnitude of the hadronization time (or, to be more precise, the time to build up the hadronic wavefunction) is based on a relativistic and simple quantum mechanical estimate DKMT  for light quark hadronic states

$`t_q^{\mathrm{hadr}}`$ $`ER^2`$ (1a)
or heavy quark meson states
$`t_Q^{\mathrm{hadr}}`$ $`{\displaystyle \frac{E}{m_Q}}R.`$ (1b)
Here $`E`$ denotes the energy of the leading ‘parton’ or jet, $`R`$ states the transversal size of the to be formed hadrons and $`m_Q`$ is the mass of a heavy quark. Taking (1a) with the average radius of the pion $`R_\pi 0.5\mathrm{fm}`$ or (1b) substituting $`m_Q`$ by $`m_\rho `$ for a $`\rho `$ meson and taking $`R_\rho 0.8\mathrm{fm}`$, one has for the formation time a simple understanding given as

$$t_F\mathrm{\hspace{0.17em}1}\mathrm{}1.2(E/\mathrm{GeV})\mathrm{fm}/\mathrm{c}.$$
(2)
Hence, for leading hadrons with moderately high $`p_{}10\mathrm{GeV}`$ original point–like jet–partons have established already a complete nonperturbative, transversal wavefunction after traveling a distance in the vacuum of length being smaller or equal than $`10\mathrm{fm}`$. Accordingly, the jets should, to a large fraction, materialize into hadrons still inside the expanding fireball, which has a transversal size of roughly

$$R(\tau )\mathrm{\hspace{0.17em}8}\mathrm{fm}+0.5\tau ,$$

where $`\tau `$ denotes the local proper time of the longitudinally (Bjorken like) and transversally expanding system with velocity $`v_t0.5c`$.

It is clear, though, that the parton cannot materialize into hadrons if it still resides in a deconfined phase. The dressing with a potential wavefunction and its subsequent fragmentation can only occur outside of deconfined matter, i.e. either in the hadronic system or in vacuum. The physics of hadronization out of a deconfined plasma for fast partons is not understood. Thus the estimate (2) should really be seen as a crude guess. The major fact we want to stress is that the final energetic hadrons (with a moderately high transverse momentum $`p_{}`$) observed at RHIC could well be realized as nearly fully established hadrons already inside the late stage hadronic fireball. If this is the case, these (pre–)hadrons will interact (by collisions) with the bulk hadronic matter making up the fireball GGX02 . In addition, refering to the phenomenology of colour transparency, the colourless objects becoming leading hadrons might interact as pre–mesons with the surrounding low momentum mesons with a cross section of roughly $`\sigma (t)=\sigma _0(t/t_F)`$ as long as their wavefuction is not completely being built up DKMT . This picture of potential late hadronic final state interaction bears again a close similarity to the one applied to explain the $`J/\mathrm{\Psi }`$ suppression pattern by either hadronic comovers gavin  or by early hadronic string–like excitations as early comovers jochen .

The typical total (elastic plus inelastic) cross section of two mesons is given by $`\sigma _010\mathrm{}20\mathrm{mb}`$. Indeed, at lower energies for the collisions, there might be strong resonance contributions, which give rise to a much larger cross section. The density of hadrons in the late fireball changes from about $`1\mathrm{fm}^3`$ down to $`0.1\mathrm{fm}^3`$. The mean free path of the ’fast’ hadron is then estimated to be $`\lambda 1\mathrm{}10\mathrm{fm}`$. Accordingly, a few collisions $`L/\lambda =0,1,2,3\mathrm{}`$ should take place. The system is potentially rather opaque! In return, a strong degradation of the momenta of the leading hadrons will result. This states the major motivation for the present study: What is the potential influence of a specified finite number of undergone collisions $`L/\lambda =0,1,2,3\mathrm{}`$ of the high momentum particles on their spectrum?

Looking at the available center of mass energy for an individual collision for mesons at the same rapidity, even for a value of $`p_{}=10\mathrm{GeV}`$, one gets a $`\sqrt{s}<2\mathrm{GeV}`$, if the target is a pion (assumed being locally) at rest. For a $`\rho `$ as target, one has an invariant mass of above 2.5$`\mathrm{GeV}`$ for particles with a transverse momentum larger than 3$`\mathrm{GeV}`$. This behaviour is depicted in Fig. 1. In any case, at such lower energies the dominant part of the collisions of the mesons or premesons with hadrons of the bulk system is nonperturbative and cannot be described by pQCD methods.

For the considered transverse momentum region, one has either elastic scattering, resonance scattering or also inelastic scattering resulting in a few final hadrons. In section II we explain the description of the collision and introduce the notations of a folding matrix. Since the surrounding hadronic gas is expected to consist of approximately 40% (direct) pions, 10% kaons and 40% rhos, we treat in the following a typical inelastic collision within the Fritiof scheme Fritiof  with a $`\rho `$–meson as characteristic target hadron being at rest. To also stress some model dependence, we compare with elastic scattering on a $`\pi `$ as target particle, where the scattering is taken simply as isotropic. In section III we show various results for the transverse momentum spectra at midrapidity for various hadronic particles by assuming a definite number of undergone collisions. The initial spectra are hereby generated by the Pythia code. We end our discussion in section IV with a summary and some conclusions.

## II Energy loss by (multiple) final state hadronic interactions

Matter of interest are the transverse momentum spectra $`dN_i/dyd\stackrel{}{p}_{}`$, where the index $`i`$ stands for the various hadron species and the transverse momentum vector $`\stackrel{}{p}_{}`$ has to be taken as two dimensional. In all what follows, we will assume a Bjorken–like space–time picture with azimuthal symmetry, so that the only changes of the spectra from the two dimensional transverse momentum vector are given by its length $`p_{}`$ and we therefore will skip for reasons of simplicity all dependencies on the azimuthal angle.

In the framework of this assumption we are interested in particles lying in the same rapidity window, introducing the short hand notation $`f_i(p_{})dN_i/dp_{}`$ for this distribution function. (Our notation is different from the more common form $`dN_i/p_{}dp_{}`$, which is due to practical reasons in the descriptions presented below.) In order to separate these distributions of (hard) initial hadronic particles of type $`i`$ from the dominant distributions of hadronic particles of type $`X`$ with moderate or low (i.e. ‘soft’) momenta, say with $`p_{}2`$ GeV, we also do introduce the distribution functions $`\stackrel{~}{f}_X(p_{},\varphi )`$. Here, $`\varphi `$ does not reflect the azimuthal angle in the cylindrical shaped Bjorken picture, but denotes the angle between the radial transverse direction and its actual momentum orientation. These latter distributions do describe the motion of the target particles. They are assumed to be steeply falling functions in the variable $`p_{}`$, corresponding to a rather low temperature in a picture of ‘thermalized’ matter, superimposed by some collective, transversal (‘hydrodynamical’) flow.

Within these definitions, one can express the resulting spectrum of particles of the kind $`j`$ after one collision of a (hard) particle of kind $`i`$ with a (soft) target hadron of kind $`X`$ as

$`f_j(p_{})`$ $`=`$ $`{\displaystyle \underset{i}{}}{\displaystyle 𝑑p__0f_i^0(p__0)𝑑p__X𝑑\varphi __X}`$
$`\times {\displaystyle \frac{1}{\stackrel{~}{N}_X}}\stackrel{~}{f}_X(p__X,\varphi __X)g_{ij}^X(p__0;p__X,\varphi __X;p_{}).`$
Here $`f_i^0(p__0)`$ has to be interpreted as the initial distribution of particle kind $`i`$, i.e. the distribution one would assume without any modifications. $`\stackrel{~}{N}_X`$ denotes the number of soft particles of kind $`X`$. The “folding matrix” $`g_{ij}^X(p__0;p__X,\varphi __X;p_{})`$ expresses the probability to get a hadron of species $`j`$ with momentum $`p_{}`$, if a hadron of the species $`i`$ and the transverse momentum $`p__0`$ collides with a target hadron $`X`$, which is described by the transverse momentum component $`p__X`$ and the orientation $`\varphi __X`$. Given this general structure with four independent variables, $`g_{ij}^X`$ is a quite complicated object. It depends on the model employed for the collision (see below), and which is connected with the target, as indicated by the upper index.

Since the non–transversal component of the target momentum, which is encoded in the variable $`\varphi __X`$, also leads to scattering in the longitudinal (beam) direction, equation (II) is, strictly speaking, only true if one assumes that scattering–out of a specified (tiny) region of rapidity is compensated by scattering–in from the next rapidity regions. This is the situation for the here assumed boost invariant Bjorken picture of 1–dimensional longitudinal expansion.

Setting in the following $`\stackrel{~}{f}_X(p__X,\varphi __X)/\stackrel{~}{N}_X=\delta (p__X)\delta (\varphi __X)`$, we assume for a general estimate that the low transverse momentum target is approximated by $`p_{}0`$, which is reasonable, as long as the hard momentum scale is clearly separated from the low momenta of the bulk hadronic particles. This leads to the more simple and self–explanatory folding equation for momentum degradation

$$f_j(p_{})=\underset{i}{}𝑑p__0f_i^0(p__0)g_{ij}^X(p__0,p_{}).$$
(4)
The basic problem is then effectively given by one–dimensional collisions. Here the reduced folding matrix $`g_{ij}^X(p__0,p_{})`$ has a straightforward interpretation: It indicates the probability that out of a given particle $`i`$ with transverse momentum $`p__0`$ one gets a particle $`j`$ with transverse momentum $`p_{}`$, if the collision is taken with a target $`X`$ at rest, where all this happens in the same rapidity interval. E.g. for elastic collisions $`i+Xi+X`$ ($`iX`$), one has the immediate normalisation $`_0^{\mathrm{}}𝑑p_{_{}}g(p__0,p_{})=1`$.

Multiple, subsequent collisions are described by applying iteratively eq. (4), where in a new iteration step the distributions $`f_i^0`$ have to be replaced by the $`f_i`$ resulting from the former step.

We do remark some words of caution:
(1) The separation of the transversal distribution into a “hard” contribution $`f_i`$ and a “soft” part $`\stackrel{~}{f}_X`$ and setting the latter to a delta distribution implies some drawback, as this separation is artificial at low $`p_{}`$ (in the detector one would see $`f_i+\stackrel{~}{f}_i`$) and thus leads to a “particle production” effect at low values of $`p_{}`$. This is connected with the normalization of the folding matrix $`g_{ij}^X`$, what one can see easily, if one for example looks at the simple situation of elastic scattering $`\pi ^+`$ on a $`\pi `$, i.e. $`i=\pi ^+`$ and $`X=\pi \frac{1}{3}(\pi ^++\pi ^0+\pi ^{}`$), leading to two pions with non vanishing transverse momenta. We remark here that one has to clearly stay in the discussion above a certain moderately large transversal momentum scale for an assumed clear separation of hard and low momentum scales.
(2) All our considerations assume, that a fixed number of collisions takes place, disregarding the fact, that the cross section could vary dramatically with the transverse momentum. For special collision partners, for example $`p+\pi `$, it is known PDG , that the cross section $`\sigma _{\mathrm{tot}}\left(p+\left[\pi ^++\pi ^{}\right]/2\right)`$ changes from $`30\mathrm{mb}`$ for $`p_{\mathrm{proton}}=3\mathrm{}6\mathrm{GeV}`$ to nearly five times this value at $`p_{\mathrm{proton}}1.8\mathrm{GeV}`$, if the pion is considerd to be at rest. Similarly, the elastic reaction $`\pi +\pi \rho \pi +\pi `$ also exceeds $`100\mathrm{mb}`$ for the projectile momentum $`p_\pi 2\mathrm{}2.5\mathrm{GeV}`$. This implies that, while assuming that one collision appears, at an other value of the transverse momentum (in these special cases) five collisons are to be awaited. We will comment on this point a little further in the conclusions.

### II.1 Initial distribution

For the initial distribution $`f_i^0`$ of individual hadrons $`i`$ at moderately high $`p_{}`$ when entering the final hadronic stage we do not have direct information. All various ideas, as mentioned in the introduction, might contribute and already steepen the spectrum. In order to see clearly the potential effect of the late hadronic interactions, we use individual distributions generated by Pythia v6.2 PYTHIA  gauged to individual $`pp`$–collisions.

It is essential to carefully adjust the parameters for the model calculations to data provided by the UA1 collaboration for $`p\overline{p}`$ collisions in the energy range $`\sqrt{s}200\mathrm{GeV}`$ and $`|y|<2.5`$ Albajar90 . Fitting the calculations to the data, one gets a best description of the experimental results by just tuning one parameter, the intrinsic transverse momentum distribution of the partons inside the nucleon, described by its averaged value $`k_{}(\sqrt{s})`$. The resulting spectra are displayed in Fig. 3.

In a second step one has to compare the calculated spectra with the parametrization $`(1+p_{}/\overline{p}_0)^n`$ used to describe the experimental data Albajar90  and extract from this fitted parmetrization ($`\overline{p}_0(\sqrt{s})`$ and $`n(\sqrt{s})`$) the characteristic value $`p_{}(\sqrt{s})`$. With this knowledge one can extrapolate the calculations of $`p\overline{p}`$ collisions down to the energy $`\sqrt{s}=130\mathrm{GeV}`$ by adjusting the intrinsic parameter $`k_{}`$ to finally meet the experimentally extrapolated value $`p_{}(\sqrt{s}=130\mathrm{GeV})=0.385\mathrm{GeV}`$ Albajar90 . Please note that this value is completely disregarded by the assumed pointwise interpolation of fits to data from different values of $`\sqrt{s}`$ for given transverse momentum $`p_{}`$, as used eg. in PHENIX ; Ex130 .

After having adjusted the microscopic Monte–Carlo calculations for $`p\overline{p}`$ collisions at $`\sqrt{s}=130`$ and $`200\mathrm{GeV}`$ and $`|y|<2.5`$, theoretical predictions for $`Au`$$`Au`$ collisions at these energies for midrapidity are given by respecting the correct isospin average corresponding to $`0.16pp+0.36nn+0.48pn`$. It is worth noting, that these distributions do not follow a true power law behaviour, as can be seen in Fig. 3, where at $`p_{}2\mathrm{GeV}`$ and $`6\mathrm{GeV}`$ significant changes in the slope occur.

Some special care has to be taken (for example mixing of calculations with non–weighted and some with reweighted events) to expand the Pythia Monte–Carlo calculations up to values of $`p_{}`$ larger than $`10\mathrm{GeV}`$ with best statistics. Naive Monte–Carlo implies, that at $`p_{}6\mathrm{GeV}`$ reasonable computation time ends and that one would have to increase the number of events by approximately a factor of 10 every time one wants to go further in the $`p_{}`$ spectrum by an other step of 2$`\mathrm{GeV}`$, as one can see in Fig. 3.

### II.2 Folding Matrices

The collision of a high $`p_{}`$ particle with particles of the medium is implemented according eq. (4) via a folding matrix $`g_{ij}(p__0,p_{})`$. Since a unique description of scatterings at these low energies (cf. Fig. 1) is not available, model descriptions have to be choosen.

In our first description A we calculate the folding matrices with Monte–Carlo calculations using the package Fritiof Fritiof  to simulate collisions with a $`\rho `$-meson. Since it is known, that within Fritiof the relative weighting between elastic and inelastic scatterings is not in full agreement with data, we therefore choose a constant weighting according the ratio 15:85. (This is given by the averaging of $`p+\rho `$–, $`\pi +\pi `$–, $`K+\pi `$ cross sections in a pomeron/reggeon exchange picture for $`\sqrt{s}>2\mathrm{GeV}`$ and also following simple estimates in a constituent quark model.)

As a second description B we use purely elastic scattering on a pion, where the angular distribution of the outgoing particles in the cm system is assumed to be isotropic, which seems to be a at a first glance a rather strong assumption.

In the description A of the collisions we are forced to extrapolate the Fritiof prescription downwards to energies, where it is strictly not fully applicable anymore (for values below an initial transverse momentum of $`p__03\mathrm{GeV}`$ we are in the resonance region, cf. Fig. 1). Therefore (or nevertheless) we stop the calculations at a minimal initial transverse momentum of about $`p__0^{\mathrm{min}}=1\mathrm{GeV}`$. Boundary problems applying the folding procedure acccording eq. (4) due to this cut could be easily cured, if one implements a more detailed weighting of inelastic and elastic cross sections, which could already allow to bend the calculations downwards to zero transverse momentum, or if one would be able to introduce a consistent picture of scattering of hadrons from very low (resonance region) to very high momenta (hard QCD processes).

The two different model descriptions show some significant differences, if one looks at the dependence of the final transverse momentum $`p_{}`$ for a given initial transverse momentum $`p__0`$. But not only the different models, but also different final particles behave differently. This is shown in Fig. 5, where as an example the convolution matrizes $`g_{i\pi ^+/p^+}^{(\pi )}`$ and $`g_{i\pi ^+/p^+}^{(\rho )}`$ are shown for the model descritpions A and B for the fixed value of the initial transverse momentum $`p__0=5\mathrm{GeV}`$.

While all the elastic distributions (Fig. 5 (bottom)) do follow a flat behaviour, where the threshold values are dictated by $`p__0`$ and the corresponding masses of the particles, the inelastic distributions (Fig. 5 (top)) do possess a much richer structure. First they show, as expected, a strong enhancement at low $`p_{}`$ values. In addition, these curves also depict a double peak structure for the largest $`p_{}`$–values in the forward scattering region. The lower peak at $`p_{}=p__0`$ is the simple one due to elastic scattering, while the second and rather astonishing peak for $`p_{}>p__0`$ is due to a large partial cross section for $`X+\rho X+\pi `$ implemented in Fritiof ($`X`$ denotes here either the $`\pi ^+`$ in the upper left figure or the proton in the upper right one of Fig. 5). This peak is located at $`(p_{}p__0)(m_\rho ^2m_\pi ^2)(2m_\rho )m_\rho /2`$. Physically, for such collision events the leading particles do become slightly accelerated. One might argue, that the likelihood of such contributions are completely overestimated, so that Fritiof should also fail in this part of the description. On the other hand, by including them, their net effect in our later discussion would be that one actually underestimates the emerging depletion of the transverse momentum spectrum. We have decided to omit these contributions from the inelastic part by hand compared to the basic Fritiof concept.

Considering eq. (4), one realizes the fact that the behaviour of the folding matrices for a fixed value of the resulting transverse momentum $`p_{}`$ as a function of the initial transverse momentum $`p__0`$ is of crucial importance. The corresponding figures for the examples from above are depicted in Fig. 5.

While the curves for the elastic scattering on a pion are smoothly falling for $`p__0p_{}`$ (or the corresponding threshold value), the curves for the scattering on a $`\rho `$ are nearly constant, but also showing the double peak structure at $`p_{}=p__0`$ and $`p_{}>p__0`$, as mentioned above. For larger $`p__0`$ the behaviour of all these folding curves gets less and less important for the final spectrum, since one has to fold them with the ininital distribution, whose steep slope is also indicated in Fig. 5. This means that the correct way of looking at the stopping is the function $`g_{ij}(p_{}=const)`$ instead of $`g_{ij}(p__0=const)`$. Looking at the potential size of the energy loss (at a partonic or a hadronic level), the value $`p__0p_{}`$ of $`g_{ij}(p__0,p_{})`$ is strongly skewed toward small values by the steeply falling distribution $`f_i^0(p__0)`$. Only the least stopped particles (i.e. particles with $`p__0`$ only somewhat larger than $`p_{}`$) in $`g_{ij}`$ contribute to the folding integral eq. (4). Simple estimates give, that values of $`g_{ij}`$ with $`p_{}<\frac{1}{2}p__0`$ contribute only at a one percent level. Although appealing from the various shapes of the folding matrix $`g_{ij}`$, assumptions like $`\mathrm{\Delta }E=\mathrm{\Delta }p_{}=p_{}p__0constp__0`$ for an averaged energy loss are much too oversimplistic.

These findings are in complete agreement with the considerations of Baier , where the energy loss of (very energetic) partons in partonic matter due to gluon radiation, calculated via perturbative QCD, is expressed via a function $`D(ϵ)`$, which has to be integrated in a sensible way with $`ϵ`$ as a measure of the energy shift.

Inspecting further Fig. 5, the $`p_{}>p__0`$ peak of the Fritiof results will lead to some (maybe unphysical) “unstopping”, as already emphasized. As those contributions are indeed then weighted with the much larger distribution $`f(p__0)`$ at $`p__0<p_{}`$, they would indeed be of more crucial importance than naively expected. Including those contributions would roughly weaken the deacceleration effect in the spectrum by a factor of 2 for $`p_{}3\mathrm{GeV}`$ and lead to a zero net effect at $`p_{}2\mathrm{GeV}`$. As already discussed, we will exclude these contributions with $`p__0<p_{}`$.

## III Quenching by multiple collisions and comparison with data

Employing the appropriate scaling of binary collisions ($`N_{coll}=945,\mathrm{\hspace{0.17em}905},\mathrm{\hspace{0.17em}19},\mathrm{\hspace{0.17em}3.7}`$ for the centralities 0…5%, 0…10%, 60…80% and 80…92% respectively), the Pythia calculations are in perfect agreement with the (at that time preliminary) PHENIX data PHENIX2  for charged hadrons and for neutral pions in the case of peripheral collisions, as can be seen in Figs 7 and 7. We can thus conclude that our Pythia calculation with the correct extrapolation provides an accurate basis for the experimental comparison.

Looking at the most central data, one recognizes the celebrated significant discrepancy between our scaled Pythia results and the data. The data resides below the theoretical curve for $`p_{}<1.5\mathrm{GeV}`$ and for $`p_{}>2.5\mathrm{GeV}`$. Only the latter case for $`p_{}2\mathrm{GeV}`$ we now want to adress by means of our scenario of possible (multiple) collisions of the jet–like (pre–)hadrons with the ‘soft’ hadrons sitting in the bulk matter of the fireball. There is some general agreement that the hadron spectra below $`p_{}2\mathrm{GeV}`$ should be dominated by soft physics and hydrodynamical expansion with some appreciable collective outward flow. In addition, in this lower momentum regime our basic assumption of a hadronic target being approximately at rest becomes no more valid. Here not only energy loss, but probably also some collisional energy gain (on a partonic and/or hadronic level) and also hadronic particle annihilation and (re–)absorption occurs which influences the shape of the spectrum.

Refering to eq. (4) and inspecting the two Figs. 5 and 5, one recognizes that the depletion in the spectrum by one collision should be stronger for description A with (in)elastic collision on a $`\rho `$ as the soft target particle, as one would also intuitively expect. This is so as the foldings matrices $`g_{ij}`$ for $`p__0>p_{}`$ are somewhat smaller for this first collision picture, though there is the same strength at $`p__0p_{}`$. In Figs 7 and 7 the three solid lines depicts the spectrum (for $`p_{}>2\mathrm{GeV}`$) for the most central region, when the particles have suffered on average 1,2,3 (in)elastic hadronic collisions on a $`\rho `$. The dotted lines in these figures indicate the effect for elastic scattering on a $`\pi `$. Indeed, for the inelastic case the overall depletion is somewhat stronger. Refering to our earlier remarks, we only do show the modifications for $`p_{}>2\mathrm{GeV}`$. One clearly recognizes the potential importance of these possible final state interactions. For $`p_{}>2\mathrm{GeV}`$ an average loss with $`L/\lambda =1\mathrm{}2`$ collisions are appropriate to explain the present experimental data both for charged hadrons and neutral pions.

We remark, that the potential energy loss by final state hadronic reactions is in the same range as considered in the various discussions for the jets in possible deconfined matter.

In Figs. 9 and 9 we now show calculations for various number of collisions at a $`\sqrt{s_{NN}}=200`$ GeV up to $`p_{}=10`$ GeV. Our results might serve also as a guideline and prediction concerning the present analysis of the 200 GeV runs at RHIC QM2002 .

In principle, the results are completely similar to the one carried out at the lower RHIC energies. Fig. 9 shows the resulting $`p_{}`$ spectra after 1,2,3 collisions according to the two collision pictures. As expected, the modification of the spectrum is similar to the one discussed for the situation at lower energies of $`\sqrt{s}=130\mathrm{AGeV}`$. Thus we can conclude that hadronic final state collisions are efficient in degrading the transverse energy of jet–like hadronic particles.

Fig. 9 gives the resulting suppression factor $`R(p_{})=\frac{dN/dp_{}|_{\text{modified}}}{dN/dp_{}|_{\text{initial}}}`$ for collisions from Fig. 9. Here the effect of the different collision schemes is more visible. With one to two inelastic collisions one will have a suppression of 50 to 80 percent for momenta $`p_{}4\mathrm{GeV}`$.

## IV Summary

We have motivated that most of the (pre–)hadrons stemming from a jet should still materialize in the dense system for transverse momenta up to 10 GeV. The late hadronic final state interactions with the bulk of comovers have a clear and nonvanishing effect in suppressing the spectrum. This is so for (in)elastic reactions on a $`\rho `$, which we have modelled via the FRITIOF scheme, and also for (isotropic) elastic scattering on a $`\pi `$. On the average one up to two such interactions should already be enough to explain quantitatively the RHIC results.

Although this comparison of the depletion of the spectra with real data is intriguing, one should be aware that our calculations are at best semi–quantitative. This is so because of a couple of simplifying assumptions invoked in the present study. We have neglected all possible dependencies of the cross sections of the colliding “hard” and “soft” (target) particles, but do assume that a fixed number of collisions has taken place, following a simple–minded opacity expansion. In particular, possible resonance contribution for the elastic scattering could enhance the effect of momentum degradation considerably. Second we made the further assumption, that the bulk hadronic targets are taken at rest. This is certainly not the case especially if one wants to quantitatively discuss the modifications in the spectra around $`p_{}=1\mathrm{}3\mathrm{GeV}`$, where both collective and ‘thermal’ dynamics should play a considerable role. Energy gain and loss mechanisms should both be at work and are probably competing in this momentum range. On a partonic level such a picture has been discussed in Wang2 , but this should be rather obvious if one thinks of a kinetic cascade–type description either on a partonic or hadronic description of the dynamics happening inside the fireball. Hence, the marriage of a transport description for the ‘soft’ matter by some standard ultrarelativistic heavy ion cascade together with the present ideas carried out in this study is one possible way to further proceed in order to obtain more quantitative results. We will leave this as a future option. A first and rather exploratory investigation in this direction has already been mentioned some time ago Bass .

Our emphasis is, however, at this stage to point out that, in order to draw some deductions for possible QCD effects of a deconfined QGP phase on the materializing jets, one has to disentangle these from the here investigated final state interactions, before definite conclusions on the importance of a potential dense partonic phase, or any other effect, can really quantitatively be drawn. This is certainly true for the case of the jets to be observed at RHIC. (At future LHC experiments observable jets with much higher $`p_{}`$ could actually materialize outside the fireball.)

Combining specific energy loss fomulae with fragmentation functions of quarks and gluons into hadrons can therefore only be understood as an extremely “effective” approach, since no conclusive answer about the relative strength of hadronic final state interactions can be extracted. The possible disentanglement, though, is very challenging and delicate. One possible way could be the precise measurements of $`p_{}`$–distributions of charmed particles like the $`D`$–mesons referee . According to a recent investigation DokKha01  heavy quarks should loose considerable less energy than light quarks or gluons (although an earlier study comes to other conclusions Munshi ). On the other hand the (elastic and inelastic) $`D`$–meson hadronic cross section should not be much smaller than that of light quark mesons. Hence, if the $`D`$–mesons should not show any significant $`p_{}`$–suppression at moderate $`p_{}3\mathrm{}10\mathrm{GeV}`$, this would then favour the partonic ’pQCD’ jet quenching scenario. At present, the PHENIX data PhenixEl  up to $`p_{}`$ values of $`3\mathrm{GeV}`$ still suffer from too large statistical errors and are compatible with full or no quenching KaiQM . In any case, the $`D`$–meson spectra will be an important experimental challenge for future RHIC undertakings. Some other valuable information on the right “initial distribution” will be learned by $`d+Au`$ experiments VG02 , eg. the size of possible shadowing and the Cronin effect. However, like the medium–induced gluon radiation the here discussed final state interactions with the late stage hadronic matter is a pure heavy ion effect. Here the $`p_{}`$–dependence of the ‘suppression’ pattern might reveal possible differences, as for higher $`p_{}`$ the time of formation for the hadronic wavefunction of the leading hadrons of a jet and the colour transparency effect should suppress these final state interactions. But this is speculation at present.

### Acknowledgments

Stimulating discussions with R. Baier, W. Cassing, T. Falter, P. Jacobs, B. Kämpfer, U. Mosel and J. Nagle are gratefully acknowledged. This work has been supported by BMBF.


text:
# References

Analytical representations for relaxation functions of glasses
R. Hilfe$`\text{r}^{1,2}`$

$`^1`$ICA-1, Universität Stuttgart, Pfaffenwaldring 27, 70569 Stuttgart

$`^2`$Institut für Physik, Universität Mainz, 55099 Mainz, Germany

## Abstract

Analytical representations in the time and frequency domains are derived for the most frequently used phenomenological fit functions for non-Debye relaxation processes. In the time domain the relaxation functions corresponding to the complex frequency dependent Cole-Cole, Cole-Davidson and Havriliak-Negami susceptibilities are also represented in terms of $`H`$-functions. In the frequency domain the complex frequency dependent susceptibility function corresponding to the time dependent stretched exponential relaxation function is given in terms of $`H`$-functions. The new representations are useful for fitting to experiment.

PACS: 77.22.Gm,61.20.Lc,71.55.Jv, 67.40.Pf, 78.30.Ly
published in:
Journal of Non-Crystalline Solids, vol. 305 (2002), page 122

Analytical representaions of relaxation functions in the time domain and susceptibilities in the frequency domain are important to fit experimental data in a broad variety of experiments on glasslike systems. Dielectric spectroscopy, viscoelastic modulus measurements, quasielastic light scattering, shear modulus and shear compliance as well as specific heat measurements all show strong deviations from the normalized exponential Debye relaxation function

$$f(t)=\mathrm{exp}(t/\tau )$$
(1)
where $`\tau `$ is the relaxation time . All relaxation functions in this paper are normalized to $`f(0)=1`$. Relaxation in the frequency domain is described in terms of the normalized complex susceptibility

$$\widehat{\chi }(u)=\frac{\chi (\omega )\chi _{\mathrm{}}}{\chi _0\chi _{\mathrm{}}}=1u\left\{f(t)\right\}(u)$$
(2)
where $`u=\text{i}\omega `$, $`\omega `$ is the frequency, $`\chi (\omega )`$ is a dynamic susceptibility normalized by the corresponding isothermal susceptibility, $`\chi _0=lim_{\omega 0}\mathrm{Re}\chi (\omega )`$ is the static susceptibility, $`\chi _{\mathrm{}}=lim_\omega \mathrm{}\mathrm{Re}\chi (\omega )`$ gives the “instantaneous” response, and $`\left\{f(t)\right\}(u)`$ is the Laplace transform of the relaxation function $`f(t)`$. For the exponential relaxation function this leads to

$$\widehat{\chi }(\omega )=\frac{1}{1+i\omega \tau },$$
(3)
i.e. the well known Debye susceptibility.

Most generalizations of equations (1) and (3) for glasses and other complex materials are obtained by the method of introducing a fractional “stretching” exponent. In the time domain this method leads to the “stretched exponential” or Kohlrausch relaxation function given as

$$f(t)=\mathrm{exp}[(t/\tau _\beta )^\beta ]$$
(4)
with fractional exponent $`\beta `$ . Of course all formulae obtained by the method of stretching exponents are constructed such that they reduce to the exponential Debye expression when the stretching exponent becomes unity. Extending the method of stretching exponents to the frequency domain one obtains the Cole-Cole susceptibility

$$\widehat{\chi }(\omega )=\frac{1}{1+(i\omega \tau _\alpha )^\alpha },$$
(5)
the Davidson-Cole expression

$$\widehat{\chi }(\omega )=\frac{1}{(1+i\omega \tau _\gamma )^\gamma }$$
(6)
or the combined Havrialiak-Negami form

$$\widehat{\chi }(u)=\frac{1}{(1+(u\tau _H)^\alpha )^\gamma }$$
(7)
as empirical expressions for the experimentally observed broadened relaxation peaks. Most surprisingly, the analytical transformations between the time and frequency domain for general values of the parameters in these simple analytical expressions seem to be largely unknown , and authors working in the time domain usually employ the stretched exponential function while authors working in the frequency domain use the stretched susceptibilities. An exception are the results in  where the real and imaginary part of the elastic modulus were obtained for Kohlrausch relaxation. Note however that there is a sign error in the real part in the results of . It is therefore the purpose of this short communication to rederive expressions for the Kohlrausch susceptibility in the frequency domain. Secondly the same methods are used to obtain for the first time the relaxation function corresponding to the Havriliak-Negami susceptibility (and hence also the Cole-Davidson and Cole-Cole susceptibilities) in the time domain. It is hoped that these expressions will be useful for facilitating the fitting of experimental data.

The objective of this paper is achieved by employing a method based on so called $`H`$-functions . The $`H`$-function of order $`(m,n,p,q)^4`$ and with parameters $`A_i{}_{+}{}^{}(i=1,\mathrm{},p)`$, $`B_i{}_{+}{}^{}(i=1,\mathrm{},q)`$, $`a_i(i=1,\mathrm{},p)`$, and $`b_i(i=1,\mathrm{},q)`$ is defined for $`z,z0`$ by a contour integral in the complex plane

$$H_{p,q}^{m,n}\left(z|\begin{array}{c}(a_1,A_1),\mathrm{},(a_p,A_p)\hfill \\ (b_1,B_1),\mathrm{},(b_q,B_q)\hfill \end{array}\right)=\frac{1}{2\pi i}_{}\eta (s)z^s\text{d}s$$
(8)
where the integrand is

$$\eta (s)=\frac{{\displaystyle \underset{i=1}{\overset{m}{}}}\mathrm{\Gamma }(b_i+B_is){\displaystyle \underset{i=1}{\overset{n}{}}}\mathrm{\Gamma }(1a_iA_is)}{{\displaystyle \underset{i=n+1}{\overset{p}{}}}\mathrm{\Gamma }(a_i+A_is){\displaystyle \underset{i=m+1}{\overset{q}{}}}\mathrm{\Gamma }(1b_iB_is)}.$$
(9)
In (8) $`z^s=\mathrm{exp}\{s\mathrm{log}|z|i\mathrm{arg}z\}`$ and $`\mathrm{arg}z`$ is not necessarily the principal value. The integers $`m,n,p,q`$ must satisfy

$$0mq,0np,$$
(10)
and empty products are interpreted as being unity. For the conditions on the other parameters and the path of integration the reader is referred to the literature  (see \[10, p.120ff\] for a brief summary). The importance of these functions for glassy relaxation arises from the facts that (i) they contain most special functions of mathematical physics as special cases and (ii) their Laplace transform is again an $`H`$-function. Moreover they possess series expansions that are generalizations of hypergeometric series.

Based on the convenient properties of $`H`$-functions it is possible to derive time and frequency domain expressions for all non-Debye relaxation functions and suceptibilities. The results of the calculations are summarized in the tables below. Table 1 gives $`H`$-function representations for all relaxation functions in the time domain. Table 2 gives the corresponding power series for all relaxation functions for small and large times. Table 3 summarizes $`H`$-function representations for the susceptibilities in the frequency domain, and Table 4 gives their power series expansions. Note that those power series where the domain of validity is given by a limit are asymptotic series. In these tables the notation

$$\mathrm{\Gamma }(a,x)=_a^{\mathrm{}}y^{a1}e^y\text{d}y$$
(11)
denotes the complementary incomplete Gamma function, and the abbreviation

$$E_a(x)=\underset{k=0}{\overset{\mathrm{}}{}}\frac{x^k}{\mathrm{\Gamma }(ak+1)}$$
(12)
is the Mittag-Leffler function. In addition the short hand notation

$$H_a(x)=H_{11}^{11}\left(x|\begin{array}{c}(1,1)\hfill \\ (1,a)\hfill \end{array}\right)$$
(13)
was introduced for writing the Kohlrausch susceptibility.

In summary the present paper has given unified representations of nonexponential relaxation and non-Debye susceptibilities in terms of $`H`$-functions. These representations lead to computable expressions that can be used to investigate the relations between the Kohlrausch susceptibility and other fit functions . It is hoped that the $`H`$-function representations given here will help to facilitate the computational transformation between the frequency and time domain in theoretical considerations and experiment.