Estimating Bounds on Collisional Relaxation Rates of Spin-Polarized 87Rb Atoms at Ultracold Temperatures

We present quantum scattering calculations for the collisional relaxation rate coefficient of spin-polarized 87Rb(f = 2,m = 2) atoms, which determines the loss rate of cold Rb atoms from a magnetic trap. Unlike the lighter alkali atoms, spin-polarized 87Rb atoms can undergo dipolar relaxation due to both the normal spin-spin dipole interaction and a second-order spin-orbit interaction with distant electronic states of the dimer. We present ab initio calculations for the second-order spin-orbit terms for both Rb2 and Cs2. The corrections lead to a reduction in the relaxation rate for 87Rb. Our primary concern is to analyze the sensitivity of the 87Rb trap loss to the uncertainties in the ground state molecular potentials. Since the scattering length for the a3Σ+u state is already known, the major uncertainties are associated with the X1Σ+g potential. After testing the effect of systematically modifying the short-range form of the molecular potentials over a reasonable range, and introducing our best estimate of the second-order spin-orbit interaction, we estimate that in the low temperature limit the rate coefficient for loss of Rb atoms from the f = 2,m = 2 state is between 0.4 × 10−15 cm3/s and 2.4 × 10−15 cm3/s (where this number counts two atoms lost per collision). In a pure condensate the rate coefficient would be reduced by 1/2.


Introduction
The recent observation of Bose-Einstein condensation (BEC) in magnetically trapped alkali atoms [1,2,3] has brought to completion a 15 year attempt to achieve BEC in a weakly interacting atomic system [4]. The success of BEC in both 87 Rb [1] and 23 Na [3] and evidence for BEC in 7 Li [2] were remarkable achievements brought about by the development of laser cooling during the past decade, the design of optical and magnetic traps for holding cold atomic samples, and most recently the development of evaporative cooling techniques to cool atoms below the recoil limit. This experimental success has renewed the interest in collisional loss rates for spin aligned alkali systems, since the binary and ternary collisional rates will limit the lifetimes of the experimental condensates. Experimentally, 87 Rb is trapped in the (f a = 2, m a = 2) state, designated the doubly polarized state or the stretched state, for which binary hyperfine changing collisions could be the dominant loss process.
Here f a and m a designate the quantum numbers for the total angular momentum of the 87 Rb atom and its projection on some convenient space-fixed axis. On the other hand, 23 Na atoms have so far been trapped in the (f a = 1, m a = Ϫ 1) state, which is theoretically expected to be more resistant to binary collisional loss. Hence, the lifetime of the 23 Na condensate will probably be limited by three-body rates, although a condensate of stretched state 23 Na atoms may be affected by binary collisional loss as well.
The purpose of this paper is to provide the most accurate calculations possible of the binary collision rates for all inelastic hyperfine scattering processes which can contribute to the loss of spin-polarized ground state 87 Rb atoms at temperatures (T < 1 K) associated with the recent experimental observation [1] of Bose-Einstein Condensation (BEC). This spin-relaxation is due to the following processes: 87 Rb(5s , f a = 2,m a = 2) + 87 Rb(5s , f b = 2,m b = 2) → 87 Rb(5s , f a' , m a' ) + 87 Rb(5s , f b' , m b' ). (1) Having all 87 Rb atoms in the stretched state with m a = f a is ideal for Bose-Einstein condensation, since inelastic collisions between such stretched states have very small rate coefficients. The entrance channel in Eq. (1) has a spin f = f a + f b of magnitude f = 4. Since the 87 Rb atom has nuclear spin 3/2, this entrance channel can only project onto the triplet a 3 ⌺ + u state of the atom pair which has total electron spin S = 1. Consequently a simple spin-exchange model of the collision [5,6] shows that stretched states do not relax during the collision. Our concern is determining the small but significant rate coefficient for the trap loss processes indicated in Eq. (1) that occur when the degeneracy of the moleculefixed projection |⍀ | = 0 or 1 of the S = 1 triplet potential is broken by relativistic forces. This splitting leads to spin-relaxation. We use standard quantum scattering methods to calculate the spin-relaxation event rate coefficient K event , defined by Stoof et al. [7], summed over all f a' , m a' , f b' , m b' channels that lead to loss of trapped atoms, namely channels for which either f a' and/or f b' 2. Such collisions lead to loss of both atoms from the trap because of the large kinetic energy release (equal to one or two units of ground state hyperfine splitting shared equally between the atoms). The total rate of spin relaxation in a Maxwellian gas (number of atoms lost per unit volume per unit time) is Ϫ2K event n 2 , where n is the density of f a = 2,m a = 2 atoms, since two atoms are lost per event represented by Eq. (1). If a condensate is present K event is multiplied by (2 Ϫ 2 )/2, where is the condensate fraction [8].
Two atomic parameters are required to describe the separated atoms: the isotopic mass m 87 = 158425.8m e (where m e is the electron mass = 9.109 ϫ 10 Ϫ31 kg) and the 6834.683 MHz splitting between the f a = 2 and f a = 1 hyperfine components. Assuming that the molecular hyperfine Hamiltonian can be adequately represented by a unitary frame transformation of the asymptotic atomic hyperfine Hamiltonians, the accuracy of the loss rate is basically limited by the accuracy of the molecular interactions we incorporate in our close-coupled scattering codes. To perform the dynamic calculations, we require four accurate molecular potentials V S,⍀ (R): one defined by the ground X 1 ⌺ + g state with S = 0 and with moleculefixed spin projection ⍀ = 0; and three defined by the lowest a 3 ⌺ + u state with S = 1 and spin projections ⍀ = 0, Ϯ 1. These potentials take the following asymptotic form [9,10,11], where V SS ⍀ = 0 (R) = ␣ 2 R Ϫ3 and V SS ⍀ = Ϯ 1 (R ) = Ϫ 1/2␣ 2 R Ϫ3 (␣ ഠ 1/137 is the fine structure constant) are the familiar spin-spin dipole terms that are primarily responsible for dipolar relaxation in hydrogen and the lighter alkalis. These spin-spin dipole terms have been elegantly treated in a series of papers by Verhaar and collaborators in Eindhoven [7,12]. The second-order spin-orbit terms V SO ⍀ which we have included are less well known in atomic collision physics, although they are well recognized as significant terms in molecular spectroscopy [13,14]. These terms, which are induced by spin-orbit interactions mediated through distant electronic states, mimic the effect of the direct spin-spin terms by introducing a splitting in the ⍀ = 0, Ϯ 1 projections of the S = 1 state and can significantly modify the spin-relaxation rates. These terms will be discussed in detail in Sec. 5, where we will show that they are of opposite sign to the ␣ 2 R Ϫ3 terms and tend to diminish the spin relaxation rate for 87 Rb.
Our goal in this paper is to systematically assess the uncertainties in the spin-relaxation rates that are introduced by various uncertainties in the molecular parameters that enter into Eq. (2), and provide realistic bounds on the possible range of the loss rate coefficient K event that can be expected for Rb 2 . All our rates are calculated using molecular potentials obtained from ab initio MCSCF codes employing the highly accurate ab initio pseudopotentials of Krauss and Stevens [10] for R < 20 a o (1 a o = 0.0529177 nm), and joined on to the well-analyzed long range dispersion potentials for the diatoms [11] for R > 20 a o. Because of the usual limitations of pseudopotentials and the typical convergence properties of ab initio calculations combined with the extraordinary demands we make on the required accuracy of the molecular potentials at ultra-cold collision energies, these potentials can only serve as excellent initial estimates of the short range portions of the V S,⍀ (R ) potentials in Eq. (2).
Collaborative work combining photoassociative spectroscopic data from the Texas group with theoretical analysis by the Eindhoven group [15,16] has done an excellent job of characterizing the a 3 ⌺ u + potential, which controls the entrance channel dynamics for the spin relaxation described by Eq. (1). In Sec. 3, we examine the sensitivity of the scattering length to variations in this potential. We will accept the analysis of the Texas/ Eindhoven group and have insured that our potential reproduces both the scattering length A 1 ഠ + 110 a o for the 87 Rb isotope and A 1 ഠ Ϫ 300 a o for the 85 Rb isotope. In addition, we introduce a useful new way to associate the scattering length with the binding energy of the last bound state in an attractive molecular potential.
In Sec. 4, we demonstrate a strong sensitivity of the relaxation rate to the shape of the X 1 ⌺ g + singlet potential, when the a 3 ⌺ u + scattering length is kept fixed at its known value. In this case the sensitivity is not due to the a 3 ⌺ u + entrance channel, but depends on a final-state close coupling effect in the exit channels. The sensitivity to the potential leads to a factor of six uncertainty in the rate coefficient, which is analyzed using generalized MQDT theory [17] and especially the associated half collision amplitude version [18] of the theory. We also describe the interplay between the spin-spin (SS) and second-order spin-orbit (SO) contributions to the relaxation rate.
We previously had speculated that the SO terms would modify Rb spin-relaxation [19]. In Sec. 5 we present new ab initio calculations of these SO terms for Rb 2 and Cs 2 . We also give our current understanding of the uncertainty range of the spin-relaxation rate of 87 Rb. Finally, a summary of our results is presented in Sec. 6.

Scattering Theory of Ground State Alkali Atoms
In a field-free collision both the total angular momentum F = f a + f b + ᐍ = f + ᐍ and the total parity p = Ϯ 1 are constants of the motion. The parity is the symmetry associated with the inversion of all electron and nuclear space-fixed coordinates through the center-of-mass of the dimer. The collision loss rate coefficient K event (E ) for total collision energy E is expressed in terms of sums over F and p involving the scattering matrices S (F ,p ,E ) [8]. Each S -matrix is derived from a multichannel wavefunction calculated from standard close-coupled codes for a given F ,p , and E, using an expansion in a channel state basis |F ,M ,p ;ᐉ , f , f a , f b ͘, which describes the asymptotic properties of the separated atoms, where ᐉ is the angular momentum (partial wave) quantum number of the interatomic coordinate R , f represents the magnitude of the channel angular momentum f = f a + f b , and ␥ gives the spin channel in which the collision starts. The + indicates normal scattering boundary conditions for an incoming state in channel ␥ and outgoing spherical waves in channels ␥' . The channel states are symmetrized with respect to interchange of the identical nuclei. One consequence of this symmetrization is that odd partial waves are missing from Eq. (3) for a collision of two f a = 2, m a = 2 87 Rb atoms.
In the absence of spin-spin and second-order spin-orbit interactions, the molecular Hamiltonian of two colliding ground state alkali atoms in (ns) orbitals, which can be expected to have zero total electronic orbital angular momentum L = (ᐍ a + ᐍ b ) = 0, possesses two additional almost good quantum numbers. These are ᐉ and f . Although this is not true when one or more of the atoms possesses electronic orbital angular momentum, such as an alkali in its first excited (np ) orbit, it is an excellent approximation for the two colliding Rb(5s ) atoms in Eq. (1). The physical reason is that for L = 0 there are no electrostatic interactions that cause locking of the electron spin angular momentum of the system to the internuclear axis. Ultimately, weak spin-spin dipole (SS) and second-order spin-orbit (SO) interactions cause the total electronic spin S = (s a + s b ) to couple to the axis and lead to the small energetic splittings between the molecule-fixed spin projections ⍀ represented in Eq. (2). For the moment, if we ignore these latter interactions the ⍀ projections are perfectly degenerate and we can easily transform the asymptotic channel states |F ,M ,p ;ᐉ , f , f a , f b ͘ into a basis defined by the total electron spin angular momentum S (S = s a + s b ) and the total nuclear spin I (I = i a + i b ), where s a and s b are the atomic electron spin angular momenta and i a and i b are the nuclear spin angular momenta (for 87 Rb, s a = s b = 1/2 and i a = i b = 3/2). This transformation is

S I f
In this new representation two adiabatic Born-Oppenheimer potentials, uniquely identified by the quantum numbers S = 0 (i.e., the X 1 ⌺ + g state) and S = 1 (i.e., the a 3 ⌺ + u state), appear on the diagonal of the Hamiltonian matrix. In the absence of SS and SO interactions the block of S channel states and the block of I channel states are diagonal and only couplings involving a simultaneous change in I and S are introduced by the hyperfine interactions. These coupling are constrained to subblocks which insure that f = I + S is conserved, and we find both f and ᐉ remain perfectly good quantum numbers at all distances. Of course, at small distances the exchange splitting between the molecular potentials is large compared to the hyperfine splittings and hyperfine coupling is negligible. However, as we shall see, at distances of the order R ഠ 20 a o to 40 a o , the hyperfine interaction becomes important and the I £ S coupling drives the system back into the asymptotically diagonal basis of channel states |F ,M ,p ;ᐉ , f , f a , f b ͘.
As seen in Eq. (2), both SS and SO interactions produce an energy splitting of the ⍀ states which implies a locking of S to the internuclear axis. This effect is highlighted by applying a frame transformation of the |F ,M ,p ;ᐉ , f ,S ,I ͘ channel basis such that ⍀ becomes represented as a ''good'' quantum number. However, in this new basis f and consequently ᐉ are no longer conserved. Fortunately the splitting of the degeneracy of the a 3 ⌺ + u state is small, and hence f and ᐉ still remain good approximate quantum numbers. This allows us to block the Hamiltonian for a given F and p into subspaces of common f and ᐉ values. The atomic stretched state angular momenta, f a = f b = 2 and m a = m b = 2, can only couple to an f = 4 state, which then couples with the ᐉ = 0 s-wave to give a total angular momentum F = 4. This is the only F value for which there is a stretched state s -wave. The spin-spin dipole interaction can only change ᐉ by two units, and thus upon examining the F = 4, p = + 1 Hamiltonian, we find that the (f = 4, ᐉ = 0) subspace can only couple directly to the (f = 3, ᐉ = 2) and (f = 2, ᐉ = 2) subspaces.
Our calculations are carried out with the full compliment of accessible channels associated with a given F ,M ,p and our rates are obtained with a sufficient summation over F ,M ,p to insure convergence. However, we find that at the temperatures relevant to BEC (T < 100 K) only the single set of F = 4, p = + 1 solutions, which involve the close coupling of 20 channels in Eq.
(3), contribute significantly to the stretched state spin relaxation. This is because only the incident s -wave contributes to spin relaxation, since the contributions from incident channels with ᐉ Ն 2 are strongly suppressed at these temperatures due to quantum threshold effects. Furthermore, of the 20 channels contributing to Eq. (3) for F = 4, p = + 1, and coupled by the 20 ϫ 20 interaction matrix U ␥,␥' (R ), only the five diabatic channels labelled ␥ = 1-5 in Table 1 and consisting of the three subspaces described above, play any significant role. Figure 1a shows the diagonal interaction potentials U ␥,␥ (R ) for these five channels. Because of complicated curve crossings and strong interactions in the ␥ basis, more insight comes when we examine the five ␣ = 1,5 ''adiabatic'' potentials V ␣ (R ) shown in Figs. 1a and 1b. These are obtained by diagonalizing the 5 ϫ 5 interaction potential U ␥,␥' (R ) at each R . The five potentials never cross and are labeled in order of increasing energy. Each adiabatic channel ␣ correlates asymptotically with the corresponding ␥ state in Table 1. At short R each ␣ potential corresponds to a very good approximation with either the X 1 ⌺ + g or a 3 ⌺ + u potentials. Thus, at short R , and out to where the exchange interaction term C exc in Eq. (2) remains dominant, the V ␣=1 (R ) potential essentially mimics the pure X 1 ⌺ + g potential, and the other four V ␣ (R ) are basically pure a 3 ⌺ + u potentials.   Table 1. On the scale of this figure the two ␥ = 3,4 and the four ␣ = 2,3,4,5 channels essentially track the a 3 ⌺ + u potential. The ␥ = 1,2,5 potentials are admixtures of V0,0 and V1,⍀ and are strongly coupled by off-diagonal terms proportional to the exchange term in Eq. (2). The adiabatic ␣ = 1 channel tracks the X 1 ⌺ + g potential.  Table 1.
Solving the set of five close-coupled equations and obtaining the 5 ϫ 5 scattering matrix S (␥ ,␥' ) = S(␣,␣') for this abbreviated set of channels is sufficient to quantitatively reproduce the elaborate multichannel calculations to within a few percent at all energies below 100 K. What has been defined [8] as the ''event'' rate constant K event is simply the sum over all inelastic events experienced by the incident stretched state channel ␥ = ␣ = 4, where is the reduced mass of the Rb 2 dimer, and k = ͙(2⑀ /ប 2 ) is the wave number of the ␣ = 4 channel incident with relative kinetic energy ⑀ . At threshold the first three terms |S ␣,4 | 2 ␣ = 1,2,3 vary as k and their contribution to K event approaches a constant at low energies. These inelastic elements measure the coupling to the exothermic channels and invariably cause loss of two Rb atoms from the trap. The fourth term produces disorientation of the stretched-state atoms which ultimately leads to loss from the trap as well. However, as this element vanishes as |S 5,4 | 2 ϰk 2 , its contribution is negligible as k → 0 and can be neglected. The close coupled equations can be numerically solved using either the ␣ or ␥ basis. Since these basis sets are asymptotically equivalent , they lead to exactly the same S -matrix. In actual practice the ␥ basis is vastly more convenient in solving the close-coupling, while the ␣ basis is much more useful in gaining physical insight from the results. For example, examining the non-adiabatic coupling [17] we find that only channels ␣ = 1 and ␣ = 2 are strongly coupled, and this occurs over a very limited region R ഠ 20 a o to 25 a o (see below). All the remaining couplings are weak and perturbative. This latter feature will play an important role in assessing the sensitivity of the rates to the singlet potential in Sec. 4.

Sensitivity Analysis of the Triplet Scattering Length
All our rates are calculated using molecular potentials obtained from ab initio MCSCF codes employing the highly accurate ab initio pseudopotentials of Krauss and Stevens [10] for R < 20 a o , and joined on to the well-analyzed long range dispersion potentials for the diatoms [11] for R > 20 a o. Since the calculation of ultracold collision rates requires extraordinary accuracy of the molecular potentials, these potentials can only serve as excellent initial estimates of the short range portions of the V S,⍀ (R ) in Eq. (2). In particular, they are not accurate enough to confidently calculate threshold properties such as the scattering length. These parameters are determined by an integration of V S,⍀ (R ) over the entire range of R .
Since we have confidence in the long range parameters in Eq. (2), we only vary the short range potential in order to assess the sensitivity of the scattering calculations to these potentials. We have added an adjustable harmonic-like short range term to the a 3 ⌺ u + and X 1 ⌺ g + potentials as follows, where Eq. (5) is only applied at distances smaller than the indicated equilibrium internuclear distances for the attractive singlet and triplet potentials respectively. Note that the functional form in Eq. (4) is chosen arbitrarily and has no physical significance; other short-range forms could be used with equal effect. In particular we want to show that by choosing an appropriate value for C 1 in Eq. (4b) we insure that the associated V 1,⍀ (R ) potentials give a scattering length of 110 a 0 as suggested by the joint theoretical-experimental analysis of Boesten et al . [16].
Our initial fit to the ab initio calculations [10] for the a 3 ⌺ + u potential happened to support 38 vibrational levels. The last level v = 37 has a binding energy of Ϫ 1.45618 ϫ 10 Ϫ8 au (1 au = e 2 /a o = 4.359748 ϫ 10 Ϫ18 J, where e is the electron charge), and a positive scattering length of 21.6 a o . This corresponds to the point at C 1 = 0 in Fig. 2 where we show the variation in scattering length as we systematically vary C 1 in Eq. (5b). The figure shows that the scattering length A 1 is extremely sensitive to the short-range portion of the potential. The scattering length is defined by the s -wave threshold behaviour of the elastic scattering phase shift 1 → Ϫ kA 1 in the limit that the asymptotic wavenumber The variation of 87 Rb2 scattering length A1 for the entrance s -wave channel incident on the a 3 ⌺ + u potential, as a function of the short range C1 parameter, extracted from the wavefunction at ⑀ /kB = 0.1 nK (full circles). The solid line through the points gives the fit of Eq. (6). We also plot the threshold values v (0) of the bound state phase as we vary C1 (full squares). The modular values of v (⑀n ) = n identify the bound state eigenvalues. As the last bound state ⑀n eigenvalue approaches zero the predicted scattering length passes from a large positive value to a large negative value as the level ''pops'' out of the potential. This occurs at about C1 = Ϫ 2.6 ϫ 10 Ϫ5 au/ao 2 for 87 Rb2 where v (0) → 38. k goes to zero (ប 2 k 2 /2 = ⑀ with ⑀ the collision energy). It is well known [20] that the actual value and the sign of A 1 is critically dependent on the position of the last bound state that can be supported by a given potential. This in turn is related to what we like to call [17,21] the bound state phase v ( ), which is defined for negative energies ⑀ = Ϫប 2 2 /2 which lie below the threshold at ⑀ = 0, and where in turn is defined as a continuous positive real variable. The modular-value of the bound state phase v ( n ) = n locates the position of the n th vibrational eigenvalue [21,25], such that ⑀ n = Ϫប 2 n 2 / 2 . Actually the deviation of the threshold value of v (0) from modular-is a useful measure of the last bound state position. This quantity plays a prominent role in many descriptions of threshold behaviour (see Stwalley [22], LeRoy and Bernstein [23] and the Eindhoven group [16]), where it has been denoted as v D and is sometimes called the effective vibrational quantum number at the dissociation limit. For our initial fit with C 1 = 0 we found v (0) = 37.699541. This quantity increases or deceases monotonically as we systematically make the potential more or less attractive by varying C 1 in Eq. (4b). At 9.75 a o , the zero energy turning point of the a 3 ⌺ + u potential for C 1 = 0, a value of C 1 = Ϯ 5 ϫ 10 Ϫ5 au/a o 2 produces a Ϯ 50 cm Ϫ1 change in the potential. For comparison, the triplet state potential is 204 cm Ϫ1 deep at its potential minimum R e1 . Figure 2 shows that the scattering length as a function of the C 1 shift parameter passes from plus to minus infinity as the last bound state is pushed out of the potential. The position of the singularity is easily located by examining the threshold behaviour of the bound state phase v (0). We see that this quantity approaches 38 just as C 1 approaches Ϫ 2.6 ϫ 10 Ϫ5 au/a o 2 and A 1 passes through infinity. In fact a nice analytic relationship exists between the scattering length A and v (0), The parameter s in cot(/ (s Ϫ 2)) is defined by the leading asymptotic power law Ϫ C s /R s for the potential. Both the a 3 ⌺ + u and the X 1 ⌺ + g potential have as the leading term C 6 R Ϫ6 and we have cot(/(s Ϫ 2)) = 1. (Instead of the pure a 3 ⌺ + u potential we actually prefer to use the adiabatic potential designated as ␣ = 4 in Table 1, in which case the lower order ␣ 2 R Ϫ3 terms in Eq. (2b) and (2c) are rigorously removed by the diagonalization of the interaction matrix and do not contribute to the threshold behaviour of this s -wave channel). If we evaluate v (⑀ ) at an eigenvalue ⑀ = ⑀ n then v (0) ഠ n Ϫ Ѩv Ѩ n = n + ␦ n . In the special case, as ⑀ n → Ϫ 0 and the last bound state lies just below the dissociation limit, such that ␦ n ഠ Ϫ Ѩv Ѩ n << 1 and tan␦ n → ␦ n , we obtain the usual perturbative expression (Ref. [20] p. 48) for the scattering length Although the conventional derivation of this expression is limited to a potential with a single bound state, we find perfect agreement with this limiting behavior even for wells supporting many bound states. Actually, if we have a situation where v (0) is not quite ready to support the n th bound state, i.e., v (0) = n Ϫ ␦ n such that ␦ n ≡ Ϫ Ѩv Ѩ ( n ) << 1, then we can visualize a ''pseudobound'' state lying just above the dissociation limit with an ''eigenvalue'' ⑀ n = + ប 2 2 n /2 . In the limit where tan v (0) → Ϫ ␦ n we obtain an expression which complements Eq. (7), A → Ϫ 1 n , and predicts a negative scattering length whenever a pseudobound state lies just above the dissociation limit. This behavior is well substantiated, and quantitatively confirmed, by the results in Fig. 2 Fig. 2 and find perfect agreement over the entire (modular ) range of v (0) with the calculated points. Note that this expression also predicts the exact locations of the zeros in the scattering length, which are always located at v (0) = n + 0.75. It is interesting to note that over the modular range the scattering length is predicted to be positive for 3/4 of the range, and negative for only 1/4. This means that, if we know nothing about the short range potential, we can at least predict there is a 3:1 probability that the scattering length will be positive! The functional form of Eq. (6) was also confirmed for an asymptotic R Ϫ8 by setting C 6 = 0 in Eq. (2), such that cot(/(s Ϫ 2)) = cot(/6) = ͙3.
If we choose C 1 = 3.128 ϫ 10 Ϫ5 au/a o 2 we obtain a scattering length A 1 (87) = + 109.1 a o and a threshold value of v 87 (0) → 37.3761, which predicts 38 bound vibrational levels with a last bound state at E 37 = Ϫ 2.947925 ϫ 10 Ϫ9 au. This choice is made to conform to the scattering length obtained by the Texas/Eindhoven group [15,16]. Further confidence is obtained from Fig.  3 where we compare the scattering lengths for the 87 Rb and the 85 Rb isotopes. In calculating the two scattering lengths, the only difference is the mass of the two isotopes. We see that at the same value of C 1 the calculated scattering length A 1 (85) = Ϫ 309.1 a o for 85 Rb is in good agreement with [15,16]. In addition, the threshold value v 85 (0) = 36.9367 predicts 37 vibrational levels with a last bound state at E 36 = Ϫ 3.305245 ϫ 10 Ϫ8 au. One final confirmation of the validity of the a 3 ⌺ u + potential we are using is obtained by examining the d -wave shape resonance structure defined by the adiabatic channel ␣ = 5 in Table 1. This corresponds to an incident channel which correlates with the triplet state at short distance, and correlates with the f a = f b = 2 atomic states at large distance entering with an ᐉ = 2 partial wave. This gives rise to the centrifugal barrier shown in Fig. 4, with a barrier height of 420 K. The radial wavefunctions f ␣=4 (R ) and f ␣=5 (R ) are shown in Fig. 5 for four incident kinetic energies: 100 K, 200 K, 350 K, and 500 K. These f ␣ (R ) are the single-channel, energy normalized elastic scattering wavefunction associated with the V ␣ (R ) adiabatic potential. Boesten et al.  [16] have concluded that there is a shape resonance enhancement of the photoassociation from channel ␣ = 5 relative to the ␣ = 4 channel with incident ᐉ = 0. The solid curve shows the ␣ = 4 wavefunction for 350 K. The ␣ = 5 channel amplitude increases with energy up to ⑀ ഠ 350 K, after which it decreases. One would infer from these plots that the d -wave shape resonance is very broad with a ''width'' extending from 200 K to 500 K. This behaviour is consistent with the analysis given in reference [16].

Sensitivity of Loss Rates to the Singlet Potential
As we systematically varied the short-range X 1 ⌺ + g potential by varying C 0 in Eq. (5) we find a surprisingly large, and seemingly erratic variation in the stretched state loss rate coefficient. The spin-relaxation rates are determined by the very small S (␣ ,4) matrix elements, which typically have a magnitude of 10 Ϫ3 . We found that the variation was intimately associated with the variation of the large |S (␣ = 1,␣ = 2)| 2 element shown in Fig.  6. As seen from Fig. 1 and Table 1, this S -matrix element measures the probability for the f a , f b = 1,2 → f a , f b = 1,1 transition. The coupling which determines this transition probability is the short-range exchange potential, which causes the strongest and only non-perturbative inelastic event associated with the channels in Table 1. Since S (1,2) is evaluated at a total energy determined by the ␣ = 4 entrance channel, namely 0.1 K above the ␣ = 4 channel threshold, the asymptotic kinetic energies in channels ␣ = 1 and ␣ = 2 are 656 mK and 328 mK, respectively (see Table 1). Before we can understand the strong influence of the X 1 ⌺ + g potential on the perturbative S (␣ ,4) elements we must first examine the profound effect of this potential on the behaviour of S (1,2). The change to the short range potential is sufficient to make small displacements in the nodes of f ␣=1 (R ) in the coupling region which is important for determining S (1,2), but it is not immedi-short distances by increasing C 0 , with nodes being systematically pushed to larger distances, Fig. 9 shows that K event tracks the associated values of |S (1,2)| 2 with K event varying by roughly a factor of 4. The figure also shows that several of the individual components contributing to the spin-relaxation show similar qualitative variation with differences in detail.  Figure 10 shows the behavior of K event with respect to the variation in |S (1,2)| 2 . We show the dependence of the stretched state rate coefficients on the X 1 ⌺ + g potential for three different situations. The SS-only curve gives K event when we only include the usual spin-spin splitting which varies as 3␣ 2 /2R 3 . The SO-only curve shows the very small rate coefficient, almost independent of S (1,2), that results if the SS splitting is removed and only the short-range second-order spin-orbit splitting is included. In Sec. 5 we will show that the SO-only contribution primarily occurs at short distance, R < 20 a o , where significant molecular interactions can occur. Consequently, the dependence of the SO-only rate coefficient on |S (1,2)| 2 is minuscule. Finally the SS-plus-SO curve shows the rate coefficient when both SS and SO interactions are included in the close-coupling calculation. This curve demonstrates up to a factor of two reduction in the rate coefficient when compared to the SS-curve The rate coefficients in Fig. 10 are double valued as a function of |S (1,2)| 2 except at the extremes of the range. This property can be traced to Fig. 9, where we have followed |S (1,2)| 2 over a range of C 0 that changes v (0)/ by unity and ''pushes'' one bound state out of the X 1 ⌺ + g potential well. We find two identical values of |S (1,2)| 2 in this interval. The second value occurs when the f 1 function is shifted by approximately half a deBroglie wavelength in the peak region of Q 12 , thereby maintaining the same value of the distorted wave integral for |S (1,2)| 2 (see Fig. 8). This changes the sign of S (1,2) but results in the same |S (1,2)| 2 . However, this change in sign results in slightly different interference effects in the evaluation of S (␣ ,4) and thus K event at distances beyond R ഠ 26 a o , as we will discuss below. Actually, in Fig. 10 we varied C 0 enough to ''eject'' two bound states from the potential, and, for a given resultant |S (1,2)| 2 the rate coefficients in Fig. 10 can not be distinguished.
At the extrema of K event in Figs. 9 or 10, it is possible to associate a single A 0 scattering length for the X 1 ⌺ + g potential with the particular |S (1,2)| 2 value. This is not possible away from the extrema, since there are two values that correspond to the same K event . The minimum in K event corresponds to a scattering length for the X 1 ⌺ + g potential of A 0 = + 95 a o and the maximum corresponds to a value of A 0 = + 54 a o . If the scattering length were measured to be near one of these values, then the relaxation rate will be near one of its extreme values.
We will now present a qualitative argument why varying the magnitude of S (1,2) affects the magnitude of the spin-relaxation rate involving the S (␣ ,4), ␣ 4, matrix elements. Since the coupling is weak it is an excellent approximation to represent these elements as follows: ately obvious why this leads to such large variations in the above-threshold S (1,2). Fig. 6 shows |S (1,2)| 2 varies from a minimum of 0 at C 0 ഠ Ϫ 9.4 ϫ 10 Ϫ5 au/a o 2 to a maximum of 0.834 at C 0 ഠ Ϫ 16.4 ϫ 10 Ϫ5 au/a o 2 . This variation is best understood by considering the matrix element of the radial coupling operator Ѩ ѨR [17] between the adiabatic states ␣ = 1,2 in Table 1 and Fig. 7. In the first order distorted wave approximation (Ref. [20] p. 349, and Refs. [26 and 27]) the S (1,2) matrix element is proportional to the integral, where Q 12 (R ) is determined by the R -variation of the orthogonal 5 ϫ 5 matrix M ij (R ) which diagonalizes the diabatic interaction matrix, Not surprisingly we find the coupling is highly localized in the vicinity of R ഠ 22 (see Fig. 7) where the spin-exchange splitting between the singlet and triplet potential in Eq. (2) becomes comparable to the hyperfine splitting of the atoms. Although Q 12 causes strong nonadiabatic mixing between channels 1 and 2, the distorted wave approximation above is suitable for the qualitative argument we make below, even though it is not suitable for quantitative calculations. Fig. 7. The adiabatic potentials for channels ␣ = 1 and ␣ = 2 together with the non-adiabatic coupling operator (in arbitrary units) between these channels. The latter is well localized in the region of 18 a0 < R < 26 a0, where the exchange splitting between the X 1 ⌺ + g and a 3 ⌺ + u potentials becomes comparable to the hyperfine splitting.
Changing the singlet potential has a negligible effect on the adiabatic f 2 (R ), which is shown by the solid curve in Fig. 8 and primarily portrays a pure triplet state, at least up to the vicinity of strong coupling near 22 a o . The two different short range singlet potentials, associated with the indicated minimum and maximum |S (1,2)| 2 in Fig. 6, are used to obtain the two different f 1 functions which we designate as f 1,min and f 1,max in Fig. 8. Note the structure of these functions in the vicinity of the peak of the non-adiabatic coupling operator Q 12 . The f 2,min function has an almost perfect overlap with the f 1 function. Since the complete non-adiabatic operator is equal to Q 12 times the radial derivative Ѩ/ѨR this perfect overlap implies a very poor overlap between f 1 and Ѩf 2,min (R )/ѨR , and therefore implies that the distorted wave integral should be quite small. In fact, for this case our exact close-coupling results yield |S (1,2)| 2 = 0.000017 (presumably by varying C 0 slightly we could have found a perfect cancellation with |S (1,2)| 2 ≡ 0). The second function f 2,max in Fig. 8 is phase shifted with respect to f 1 , implying improved overlap between f 1 and Ѩf 2,max (R )/ѨR and a larger distorted wave integral. For this case our close-coupling predicts the maximum |S (1,2)| 2 = 0.83. Fig. 8. The adiabatic f2 wavefunction (solid curve) and the adiabatic wavefunctions f1,max (dotted curve) and f1,min (dashed curve). The latter were calculated using the singlet shift parameters C0 = Ϫ 16.4 ϫ 10 Ϫ5 au/ao 2 and C0 = Ϫ 9.80 ϫ 10 Ϫ5 au/ao 2 , which produce the maximum and the minimum inelastic |S12| 2 elements in Fig. 7 respectively. The Q12 operator (in arbitrary units) is also shown, and locates the region of strong non-adiabatic coupling between channels 1 and 2.
Our initial fit to the ab initio calculations [4] for the 1 ⌺ + g potential, combined with the 3 ⌺ + u potential with a scattering length A 1 = 109 a 0 prescribed by Ref. [16], just happened to yield a value of |S (1,2)| 2 = 0.0731. When the 1 ⌺ + g potential was systematically modified at where F Ϫ ␣'␣ are defined in the same manner as in Eq. (3), but now in the adiabatic channel basis and with an outgoing state in channel ␣ . Each column vector in the matrix F Ϫ (R ) defines the radial components of a fivechannel close-coupled outgoing wavefunction for each of the adiabatic channels ␣ = 1,5. The function f + 4 represents the incoming wavefunction of the ␣ = 4 adiabatic potential of Fig. 1. This can be viewed as a generalized, multichannel version of the Distorted Wave Approximation [20,26,27]. The Q ␣',4 matrix elements introduce the weak spin-dependent coupling between adiabatic channels 4 and ␣'.
At short distances, to the left of the strong coupling region in Fig. 8 the solutions F Ϫ (R ) are simply proportional to the 5 ϫ 5 diagonal matrix of adiabatic reference functions f ␣ (9a) Applying generalized MQDT theory [17] and its associated half collision amplitude [18] to the outgoing multichannel functions at distances beyond about R = 26 a o , the exact close-coupled wavefunctions can be represented rigorously as follows where g (R ) = ␦ ␣,␣' g ␣ (R) and g ␣ (R) is an irregular solution for the adiabatic V ␣ (R ) potential. Alternatively, this may be written using running wave reference functions: where the real symmetric Y matrix is related to the close-coupled scattering matrix S = e i ⌺e i = e i (1 + iY ) The only Y matrix element of any magnitude in this ultracold five channel system is Y 1,2 ; therefore, we need only consider the ␣ = 1 and 2 channels and we can reduce F Ϫ to a 2 ϫ 2 matrix for these two channels. For the weak coupling case in Fig. 8 even the Y 1,2 element is negligible and the structure of F Ϫ (R ) remains diagonal as in Eq. (9a) for all R . However, if the chosen X 1 ⌺ + g potential yields a large S (1,2), the structure in Eq. (9b) can strongly influence Eq. (8). It is not a bad approximation [17] to represent a 2 ϫ 2 Y matrix as follows, where, especially if one uses an adiabatic representation, the diagonal elements d are generally negligible. For the strong coupling case, |S (1,2)| = |⌺(1,2)| can have a maximum value approaching unity. For the simple model used here, we use x = 1 for the strong coupling case and x = 0 for the weak coupling case. For these two extreme cases the function F Ϫ in Eq. (9b) takes a simple form. For the two channels ␣ = 1 and ␣ = 2 the relevant matrix elements of F Ϫ are: We designate these two cases ''s'' and ''w'' respectively, for strong and weak inelastic scattering probability measured by S (1,2). Using these limits in Eq. (8) we find The radial functions f 1 and f 2 are just the elastic scattering standing waves shown in Fig. 8 and oscillate strongly against the standing wave f 4 defined by the incident channel. Thus, we expect small values for the matrix elements in Eq. (12w). In the strong coupling regime the matrix elements involve the overlap of f 4 with pure outgoing (or incoming) running waves with amplitudes which do not oscillate with R , and we can easily understand why the matrix elements in Eq. (12s) are much larger than those in Eq. (12w). The enhancement of the spin-relaxation rate coefficient in the presence of strong final state interactions is well demonstrated in Figs. 9 and 10. Further insight into this final state effect can be gleaned by evaluating K event while systematically limiting the range of the spin-spin coupling in Eq. (2). We do this by assuming the fine structure constant ␣ has its usual constant value 1/137 up to some cutoff distance R cut and then vanishes identically for R > R cut . The results of such a study are shown in Fig. 11. The solid curves show the result using the X 1 ⌺ + g potential which yields the maximum |S (1,2)| 2 . In this case we expect Eqs. (10s) and (11s) to prevail at distances larger than the peak of the Q 12 operator, which is shown by the dotted curve in Fig.  11. The dashed curves show the corresponding results for the minimum |S (1,2)| 2 case, for which Eqs. (11w) and (12w) should be valid for all distances. There are two sets of curves. Each set has one line with the SS interaction only and one with our best estimate of the second-order spin-orbit terms included. Since these SO terms are short ranged and their principle contributions occur at distances to the left of the Q 12 operator, we did not apply a cutoff to these terms. We now discuss the interesting effect that the total loss rates are decreased when we include the V SO terms as well as the V SS terms. The distorted wave approximation in Eq. (8) implies that the S (4,␣ ) S -matrix elements can be separated into an SS and an SO term. The loss rate depends on the square of S -matrix elements: |S SS + S SO | 2 = |S SS | 2 + |S SO | 2 + S SS S SO* + S SS* S SO .
The first two terms sum the individual contributions to the rate, whereas the last terms exhibit the interference between them. Figures 10 and 11 show that |S SO | 2 is independent of the strength of the exit channel coupling measured by |S (1,2)| 2 and is always small compared to |S SS | 2 . On the other hand, |S SS | 2 does depend very strongly on the strength of the exit channel coupling. Obviously the interference effect is more significant for the case of weak SS coupling than for the case of strong SS coupling. In the former case, K event is decreased by a factor of 2 when SO coupling is included, whereas in the latter case, the decrease is only 20 %. Including the SO terms causes a decrease in K event , since the SO coupling has an opposite sign from the SS coupling, as discussed in the next section.

Evaluation of Second-Order Spin-Orbit Coupling
The doubly spin polarized atomic states on the left hand side of Eq. (1) can only be relaxed by the weak coupling between the two electron spins. Two terms, shown in Eq. (2), contribute to the effective spin-spin interaction Hamiltonian: the direct spin-spin dipole interaction V SS ⍀ (R ), varying at long range as 1/R 3 , and the indirect second-order spin-orbit interaction V SO ⍀ (R ). The latter originates when the atomic charge clouds overlap as a molecule is formed, and the interaction between the ground state spins are modified due to couplings mediated through distant excited electronic states of the molecule. These interactions are well known in molecular spectroscopy and mimic the direct spin-spin coupling for a 3 ⌺ state [13,14]. The reason these interactions mimic spin-spin coupling is that they split the ⍀ = 0 and Ϯ 1 components in a similar manner as the direct spin-spin terms. For heavy species like Rb and Cs we will see that these indirect terms can be much larger than the V SS ⍀ (R ) term at short distances and strongly influence the spin-relaxation rate.
We have calculated ab initio molecular spin-orbit matrix elements to obtain estimates of the second-order spin-orbit correction terms V SO ⍀ (R ) in Eq. (2). For the a 3 ⌺ u + state these terms are mediated through distant electronic states of 1 ⌸ u , 3 ⌸ u , 5 ⌸ u , and 3 ⌺ u Ϫ symmetry. We calculate the R -dependent spin-orbit matrix elements for these states using all-electron wavefunctions for the molecular states generated by standard ab initio methods. The second-order interaction for the a 3 ⌺ u + state is dominated by the matrix elements ͗ a 3 ⌺ + u,⍀ |H SO (R )| 2S+1 ⌸ u,⍀ ͘ involving 1 ⌸ u and 3 ⌸ u states which correlate to the first excited 2 S + 2 P atomic asymptotes. In this case the second-order coupling term due to a specific 2S+1 ⌸ u,⍀ state takes the form [14] V ⍀ and Here H SO (R ) is the electronic spin-orbit coupling operator and ⍀ = ⌳ + ⌺ , where ⌳ is the projection of the electronic orbital angular momentum and ⌺ is the projection of the electron spin angular momentum on the internuclear axis. Since the important aspect of the second-order spin orbit coupling for dipolar relaxation is the splitting it introduces between ⍀ = 0 and Ϯ 1 components of the 3 ⌺ u + state, we represent the second-order spin-orbit couplings as an ''effective'' spin-spin coupling term V ⍀ SO (R ) as follows, Since P S (R ) decays exponentially to large R , and its magnitude at short R is small compared to uncertainties in the short range 3 ⌺ u + potential, we can ignore the mean contribution of these interactions which equals [3P 1 (R ) + P 0 (R)]/2 and assume that our adjustment of the 3 ⌺ u + potential to fit the experimental-theoretical estimate of the scattering length in Sec. 3 actually incorporates this mean spin-orbit contribution. Figure 12 shows our calculated P 0 (R ) and P 1 (R) for Rb 2 and Cs 2 . The b 3 ⌸ u state, which is energetically closest to the 3 ⌺ u + state, has the largest coupling, and P 1 is an order of magnitude larger than P 0 . Figure 12 also shows that P 0 and P 1 have a much shorter range than the 1/R 3 spin-spin term. To a good approximation we find we can fit the numerical ab initio results to the following expression, 10 Ϫ18 J). Note that V SO ⍀ = 0 is negative compared to the positive ␣ 2 /R 3 spin dipolar term in Eq. (2b) and at about R ഠ 10.33 a o these two terms exactly cancel in the case of Rb 2 . The same is true for Eq. (2c). Since it is the difference between the potentials from Eqs. (2b) and (2c) that gives rise to spin-relaxation, and the overall magnitude of these differences is reduced by the addition of the spin-orbit interaction, we might anticipate the spin-dipolar relaxation rates will be reduced accordingly. However, we should emphasis that although this is true for spin aligned Rb, where the two contributions are of similar magnitude, the opposite effect appears likely for Cs, since the magnitude of the second-order spin-orbit term for Cs 2 , although still negative, is much larger than the spin-spin term. The overall effect of the spinorbit term in Cs 2 is to increase the spin-relaxation rate relative to that predicted by spin-spin only.
To test the sensitivity of K event to the magnitude of the spin-orbit coupling we have arbitrarily multiplied Eq. (15) by a factor of two and recalculated the loss-rates. Although it is impossible to place unequivocal error bounds on our calculations of V ⍀ so , it is unlikely that our error in estimating V SO is as much as a factor of two. Such a variation does give a reasonable bound to the possible effect of the coupling. In Fig. 13 we compare the rates for the SS-only calculation to those calculated using V SO and 2 ϫ V SO added to V SS . This figure also summarizes the status of our current confidence in K event for 87 Rb. Since we are not able to determine the 1 ⌺ + g potential with sufficient accuracy to place any constraint on the magnitude of |S 12 | 2 , the rate coefficient is spanned by the range shown by the curve labeled ''SS plus SO'' in Fig. 13. The loss rate coefficient, 2K event , lies between 0.4 ϫ 10 Ϫ15 cm 3 /s and 2.4 ϫ 10 Ϫ15 cm 3 /s. If our calculated V SO should be erroneous, we still expect the loss rate coefficient to lie between 0.1 ϫ 10 Ϫ15 cm 3 /s and 3.0 ϫ 10 Ϫ15 cm 3 /s. The lower bound is determined by the ''SS plus 2 ϫ SO'' curve, and the upper bound by the ''SS only'' curve. It must also be remembered that this assessment is based on having a precise scattering length for the a 3 ⌺ + u potential. If the estimate of A 1 = 110 a o supplied by the Texas/Eindhoven should be in error, then the rate coefficient will be modified accordingly. Fig. 13. Rate coefficient Kevent versus |S12| 2 as C0 is varied for the X 1 ⌺ + g potential. The upper two curves are the same as in Fig. 10. The three curves indicate the uncertainty in Kevent due to the uncertainty in the X 1 ⌺ + g potential and in the SO interaction. The three curves show three different cases of the strength of the second-order spin-orbit (SO) interaction: the upper curve has no SO interaction, the middle curve uses our calculated ab initio values, and the lower curve uses twice the calculated values.
We wish to make one final point before we conclude our discussion of the spin-order effects. Figure 12 shows that V SO is about an order of magnitude larger for Cs 2 than for Rb 2 . In fact, the contribution to spin-relaxation of doubly spin-polarized Cs atoms from the second-order spin-orbit term will dominate that due to the spinspin term. Model calculations for Cs suggest that spinrelaxation rates for Cs are insensitive to the singlet potential. This agrees with our result for Rb shown by the nearly horizontal curve labeled ''SO only'' in Fig.  10. The spin-relaxation rate for Cs atoms is likely to be much larger than that for Rb atoms, but the actual value will still depend on the a 3 ⌺ + u potential.

Conclusions: Assessment of Current Accuracy of Calculated Rb2 Loss Rates
In conclusion, we have determined the uncertainty in the spin-relaxation rate coefficient for stretched state 87 Rb atoms associated with uncertainties in the molecular parameters which control the magnitude of the relaxation cross section. The stretched state relaxation affects the lifetime of the experimentally observed 87 Rb condensate [1,24]. The condensate described in [1] remains dilute enough that its lifetime, ignoring collision with background thermal atoms, is determined by the binary stretched state loss rates [24]. Ternary rates will become dominant in more dense condensates.
Our calculations use the a 3 ⌺ + u or stretched state scattering length A 1 = 110 a 0 as required by the experimental photoassociation data for spin-polarized 87 Rb atoms [16]. Our calculations show that the lack of experimental knowledge of the X 1 ⌺ + g potential provides the largest source of uncertainty in determining the spin-relaxation rate coefficient. The s -wave entrance channel for collision of the f a = 2 + f b = 2 stretched state atoms only couples very weakly to the two possible exit channels of the spin-relaxation process, which produce f a = 1 + f b = 2 or f a = 1 + f b = 1 separated atoms with increased kinetic energy. The uncertainty in spin-relaxation rate is associated with the strength of mixing between the two components in these exit channels. We varied the inner wall of the X 1 ⌺ + g potential in order to determine the range of uncertainty. Our calculations show that the loss rate coefficient due to spin-relaxation, 2K event , is uncertain to about a factor of 6, lying in the range 0.4 ϫ 10 Ϫ15 cm 3 /s to 2.4 ϫ 10 Ϫ15 cm 3 /s. In a pure condensate the rate coefficient is simply K event [18]. Doubly polarized Rb may have the smallest collisional loss rate coefficient of any the alkali species if the rate coefficient lies near the lower end of its estimated range.
We have provided ab initio estimates for the secondorder spin-orbit terms V SO ⍀ (R ) which contribute to the effective spin-spin interaction. For the heavier alkali atoms these terms have an important effect on the spin relaxation rate. In 87 Rb, the V SO ⍀ (R ) terms causes a reduction in the rate coefficient for the collisional relaxation of stretched state atoms throughout the possible parameter space for varying the X 1 ⌺ + g potential. The decrease ranges between a factor of 2 % and 20 %, depending on whether the rate coefficient lies near the lower or higher end of its range of uncertainty. Our calculations suggest that the contribution of the V SO ⍀ (R ) terms to the spin-relaxation rate coefficient of stretched state Cs atoms will be much larger than that from the spin-spin dipole term.
Our analysis points out a critical need for more precise determinations of the X 1 ⌺ + g and a 3 ⌺ + u potentials for Rb 2 and Cs 2 and of the spin-coupling parameters. Experimental determination of V SO ⍀ (R ) could be approached by precision spectroscopy on the fine and hyperfine structure of the a 3 ⌺ + u state. The quantum chemistry community could provide more complete and accurate calculations of the second-order spin orbit interactions as a function of R .