Magnetic field modification of ultracold molecule-molecule collisions

We present an accurate quantum mechanical study of molecule-molecule collisions in the presence of a magnetic field. The work focusses on the analysis of elastic scattering and spin relaxation in collisions of O2(3Sigma_g) molecules at cold (~0.1 K) and ultracold (~10^{-6} K) temperatures. Our calculations show that magnetic spin relaxation in molecule-molecule collisions is extremely efficient except at magnetic fields below 1 mT. The rate constant for spin relaxation at T=0.1 K and a magnetic field of 0.1 T is found to be as large as 6.1 x 10^{-11} cm3/s. The magnetic field dependence of elastic and inelastic scattering cross sections at ultracold temperatures is dominated by a manifold of Feshbach resonances with the density of ~100 resonances per Tesla for collisions of molecules in the absolute ground state. This suggests that the scattering length of ultracold molecules in the absolute ground state can be effectively tuned in a very wide range of magnetic fields. Our calculations demonstrate that the number and properties of the magnetic Feshbach resonances are dramatically different for molecules in the absolute ground and excited spin states. The density of Feshbach resonances for molecule-molecule scattering in the low-field-seeking Zeeman state is reduced by a factor of 10.


I. INTRODUCTION
The experimental work with cold and ultracold molecules is predicted to lead to many fundamental applications in quantum computation [1], condensed-matter physics [2], precision spectroscopy [3], and physical chemistry [4]. Recent theoretical work has shown that ultracold ensembles of molecules trapped in optical lattices can be used to develop novel schemes for quantum information processing [6], design experiments for quantum simulations of condensed matter physics [2,7], engineer novel phases with topological order [8], and explore many-body dynamics of strongly interacting systems [9]. Polar molecules in external field traps may form chains, which can be used to study rheological phenomena with non-classical behavior [2]. The creation of a Bose-Einstein condensate of molecules may enable the study of Bose-enhanced chemistry [10] and the effects of symmetry breaking on chemical interactions at ultracold temperatures [11]. The realization of these proposals as well as the creation of dense ensembles of ultracold molecules depends critically on the possibility of controlling binary molecule -molecule interactions in molecular gases by external electromagnetic fields. The discovery of magnetic Feshbach resonances in atomic collisions [12,13] opened the door to many groundbreaking experiments with ultracold gases such as the realization of the BEC-BCS crossover [14], the observation of quantum phase transitions [7] and the creation of ultracold molecules using time-varying magnetic fields [15]. Extension of this work to molecular collisions may similarly lead to the development of new research directions such as cold controlled chemistry [4], quantum coherent control of bimolecular reactions [16], and quantum condensed-matter physics with molecular condensates [5].
Fueled by the promise of new discoveries, the experimental research of cold and ultracold molecular collisions is expanding rapidly. Some recent landmark experiments include the observation of threshold collision laws in Xe-OH scattering using slow molecular beams [17], measurements of cross sections for D 2 -OH collisions in a magnetic trap [18], and the detection of magnetic Feshbach resonances in collisions of Cs 2 molecules [19]. Further progress in the experimental study of complex molecule-molecule collisions and chemical reactions at cold and ultracold temperatures requires rigorous theoretical calculations elucidating the mechanisms of energy transfer in ultracold molecular collisions. Theoretical studies of molecule-molecule collisions at low temperatures are also necessary to understand the prospects for evaporative cooling of molecules in a magnetic trap. The evaporative cooling relies on elastic molecule -molecule collisions as the trap depth is gradually reduced [13]. In order to sustain efficient evaporative cooling down to quantum degeneracy, the ratio of the probabilities for elastic scattering and spin relaxation (γ) in molecule -molecule collisions must exceed 10 4 [13]. The magnitudes of the rate constants for elastic and inelastic collisions and the dependence of the scattering observables in molecule -molecule encounters on the magnetic field are, however, completely unknown.
Previous theoretical work by Krems and Dalgarno [20] and Volpi and Bohn [24] identified the main mechanisms of spin relaxation in atom-molecule collisions. It was found that spin relaxation in collisions of 3 Σ molecules with He atoms occurs through coupling to rotationally excited states mediated by the spin-spin interaction [20]. While the mechanisms of spin relaxation in molecule-molecule collisions should be similar, an accurate computational study of molecule-molecule scattering in a magnetic field is urgently needed to elucidate the rates of energy relaxation in dense ensembles of trapped molecules and understand the prospects for evaporative cooling of molecules. The quantum dynamics of ultracold molecule-molecule collisions in the absence of external fields has been studied by several authors (see, e.g., Refs. [21][22][23]). Avdeenkov and Bohn calculated the rates of spin relaxation in collisions of O 2 ( 3 Σ − g ) molecules at zero magnetic field and observed Feshbach resonances due to molecular rotation [25,26]. Krems and Dalgarno [20] presented a formalism for quantum scattering calculations of cross sections for molecule-molecule collisions in a magnetic field. However, they did not consider the symmetrization procedure required to properly describe collisions between identical molecules. Here, we extend the approach of Krems and Dalgarno to study collisions of identical 3 Σ molecules. The use of exchange symmetry reduces the number of scattering channels by a factor of two, which allowed us to obtain converged cross sections for elastic scattering and spin relaxation of O 2 ( 3 Σ − g ) molecules in an external magnetic field. We consider ultracold collisions of O 2 molecules in their ground electronic state of 3 Σ − g symmetry. O 2 -O 2 collisions play an important role in atmospheric chemistry [29,30]. The O 2 -O 2 dimer is an interesting molecular complex [27][28][29][30][31][32][33][34][35] (for a recent review, see Ref. [34]).
The binding in the complex is affected by the Heisenberg exchange interaction and the formation of incipient chemical bond [31,32]. Molecular beam scattering experiments [31,32] and high-level ab initio calculations [33] have previously been used to elucidate the nature of the chemical bond in the O 2 -O 2 dimer. Here, we calculate the magnetic field dependence of the probabilities for spin relaxation in O 2 -O 2 collisions at both ultracold (down to 10 −6 K) and cold (0.5 K) temperatures relevant for magnetic trapping experiments. Our calculations elucidate the possibility of evaporative cooling of paramagnetic 3 Σ molecules in a magnetic trap and indicate that spin depolarization in molecular collisions can be controlled by an external magnetic field. The results of our work will be useful for the interpretation of collision experiments involving slow beams of O 2 molecules produced by Zeeman deceleration [37] and cryogenic cooling [38].

A. Hamiltonian and symmetrized basis functions
The Hamiltonian for two 3 Σ molecules in an external magnetic field can be written as [20,25] where the vectors r A and r B describe the orientation of the molecules A and B in the spacefixed frame, R is the vector joining the molecules,l is the orbital angular momentum of the collision complex, µ is the reduced mass, andV el is the electrostatic interaction potential (see Sec. IIB). The asymptotic HamiltonianĤ as describes the individual molecules A and B in the presence of a magnetic field [20].
The Hamiltonians of an individual 3 Σ − molecule can be written as [39] where the index ν = A, B labels the molecule,N 2 ν is the rotational angular momentum, B e is the rotational constant, µ B is the Bohr magneton, B is the magnetic field vector andŜ ν is the electron spin. The spin-rotation and spin-spin interactions in Eq. (3) are parametrized by the phenomenological constants γ and λ SS , which do not depend on ν since we consider collisions of identical molecules. The spin-spin interaction given by the last term in Eq. (3) can be represented as a sum of tensor products of two identical spin operators and spherical harmonics, which depend on the orientationr ν of the diatomic molecule in the space-fixed (SF) coordinate frame [39]. We use the rigid rotor approximation (r ν =r ν ) and neglect the weak magnetic dipole-dipole interaction and the hyperfine interaction of 17 O [25].
We solve the scattering problem for two identical bosonic molecules by expanding the total wave function in a space-fixed uncoupled basis introduced in Ref. [20]. The uncoupled basis functions are products of the eigenfunctions of the operatorsN 2 ν andŜ 2 ν , their projections along the space-fixed z-axisN νz andŜ νz , and the spherical harmonics |ℓm ℓ where |τ ν = |N ν M Nν |S ν M Sν . The matrix elements of the Hamiltonians (1), (2) in the basis (4) were derived by Krems and Dalgarno [20] for the case of distinguishable molecules.
The basis functions for identical molecules must be the eigenfunctions of the permutation operatorP : r A → r B ; r B → r A ; R → −R, which commutes with the total Hamiltonian (1). The symmetrized basis functions can be constructed by applying the operator 1 +P to Eq. (4) and normalizing the result [20,22,40,41] where ǫ = (−) ℓ and the symmetry of the basis functions with respect to the interchange of identical molecules is given by index η which is equal to 1 for composite bosons and −1 for composite fermions [42]. In what follows we assume that the first basis function in each product on the right-rand side of Eq. (5) depends on the coordinates of molecule A, and the second basis function depends on the coordinates of molecule B [43]. Equaton (5) defines the well-ordered (combined) molecular states with τ A ≥ τ B . The normalization factor [2(1 + δ τ A τ B )] 1/2 takes into account the fact that the basis function (5) with τ A = τ B need not be symmetrized when ηǫ = +1 (i.e. when the colliding molecules are in the same state).
The basis functions given by Eqs. (4) and (5) are also the eigenfunctions of the inversion operator, with eigenvalues given by (−) N A +N B +ℓ [41]. Since for homonuclear molecules both N A and N B are either even or odd, the quantity ǫ = (−) ℓ is conserved as well. Thus, it is convenient to rewrite Eq. (5) in a factorized form where We note that no such symmetry exists for heteronuclear molecules, where the basis functions with ∆ℓ = ±1 are coupled by the interaction potential, and ǫ is not conserved.
The uncoupled symmetrized basis functions given by Eq.
where the coefficients C ǫη τ A τ B ,γ A γ B form the matrix, which diagonalizes the asymptotic Hamiltonian in the symmetrized basis (5). More explicitly, C T H as C = E, where E is the diagonal matrix of asymptotic energies and ǫ γ A γ B = ǫ γ A + ǫ γ B . The energies of individual molecules ǫ ν can be obtained by diagonalization of individual molecule HamiltoniansĤ ν (2) in the subspace of functions |τ ν . Exploiting the properties of the interchange operatorP † =P andP 2 = 1 and using Eq. (5), we obtain the following expression for the matrix elements ofĤ as This expression provides a convenient method of constructing the matrix elements of the asymptotic Hamiltonian (2) in the symmetrized basis. Alternatively, one can obtain the transformation coefficients (8) directly from the eigenvectors ofĤ A defined as The transformation can be derived by multiplying Eq. (10) by a similar expression for |γ B and rearranging the terms. The result is We have verified by numerical tests that Eqs. (8) and (11) give identical results up to an unimportant overall phase. In practice, however, it is more convenient to use Eq. (11) because it requires the evaluation and the diagonalization of smaller matrices.

B. Interaction Potential
The spin-dependent interaction potential between two 3 Σ molecules can be represented in the formV whereŜ =Ŝ A +Ŝ B is the total spin and M S = M S A + M S B . The representation (12) is often used to describe the interactions of ultracold alkali metal atoms [12]. Here, we choose an alternative approach proposed by van der Avoird and coworkers [27,28], in which the full interaction potential is separated into the spin-independent and spin-dependent partŝ The spin-dependent part, also known as the Heisenberg exchange interaction [27,28], can be parametrized for two 3 Σ molecules aŝ where J(R,r A ,r B ) is a scalar function similar to the interaction potential but decreasing exponentially with R [28]. We expand the interaction potential in the angular basis [20,27,44] with the space-fixed angular basis functions defined as where (:::) denotes a 3-j symbol. An efficient procedure for the evaluation of the radial expansion coefficients in Eq. (15) is described in Appendix A.
The matrix elements of the spin-indepenent interaction potential between the symmetrized basis functions (5) have the form The first (direct) term on the right-hand side is a matrix element in the unsymmetrized basis. Krems and Dalgarno [20] evaluated this matrix element using the expansion of V si in terms of the angular basis functions (15). Using the conservation of the total angular The second (exchange) matrix element appears as a result of the symmetrization procedure (5). It can be obtained from Eq. (18) by interchanging the indexes N ′ The matrix elements of the Heisenberg exchange interaction are obtained using the same procedure: The matrix elements of the prefactor −2J(R, θ A , θ B , ϕ) are given by Eq.
. Equation (19) shows that unlike the spin-independent matrix elements (17) and the exhange matrix element can be obtained from Eq. (20) by interchanging the indexes C. Close-coupling equations, scattering amplitudes, and cross sections The symmetrized wave function of the collision complex can be expanded as Substituting this expansion into the Schrödinger equation, Hψ η = Eψ η , where E is the total energy, we obtain the close-coupling equations for the radial expansion coefficients The asymptotic form of the solutions F η The scattering amplitude for indistinguishable molecules can be written as [23,52] q η where the scattering amplitudes on the right-hand side can be written in terms of unsymmetrized T -matrix elements Substituting this expression into Eq. (24) and rearranging the terms, we find where These T -matrix elements are exactly the same as those obtained from the solution of the close-coupling equations in the symmetrized basis set (as described above). They are related to the S-matrix elements via [48] T η The integral cross section for indistinguishable molecules in fully polarized nuclear spin states can be obtained from the scattering amplitude as described in Appendix B. The final result is where k 2 is the wave vector for the initial collision channel. This equation shows that the elastic cross section of indistinguishable bosons (or fermions) is twice as large as that of distinguishable particles. Our Eq. (29) agrees with the results presented by Avdeenkov and Bohn [25] and Burke [47].

D. Computational details
To parametrize the asymptotic Hamiltonian (2), we used the accurate spectroscopic constants of 17 O 2 (in units of cm −1 ): B e = 1.353, γ = −0.00396, λ SS = 1.985 [53]. In order to obtain converged cross sections at collision energies below 0.5 K, we used a basis set comprising 3 rotational states (N A , N B = 0 − 4) and 4 partial waves (ℓ = 0 − 6). These parameters resulted in 2526 coupled differential equations for M = 0. At the lowest collision energy of 10 −6 K, four partial waves were sufficient to achieve convergence. The coupled equations (22) were solved on a grid of R from 2 a 0 to 200 a 0 with a step size of 0.04 a 0 yielding the cross sections converged to within 10%. The spin-independent and spin-dependent parts of the interaction potential for the oxygen dimer were constructed to reproduce the integral cross sections measured in molecular beam scattering experiments [31,32]. The lowestorder isotropic part V 000 (R) was determined with high accuracy from an analysis of the glory pattern in high-velocity collisions of rotationally "hot" O 2 molecules. The anisotropic coefficients V 202 (R), V 220 (R), and V 222 (R) were inferred from the scattering experiments with rotationally cold, supersonic beams of aligned O 2 [31,32]. Recent ab initio studies [35,36] have shown that the O 2 -O 2 interaction may be more anisotropic than suggested in Ref. [31]. The inelastic transition probabilities in molecule-molecule collisions increase with the anisotropy of the interaction potential [20], so the cross sections for spin relaxation presented in this work should be viewed as lower bounds to the actual magnitudes.

III. RESULTS AND DISCUSSION
We consider collisions of two identical oxygen molecules in the ground electronic state  Fig. 2 shows that spin-changing collisions are suppressed at low magnetic fields (< 1 mT) leading to a large value of the elastic to inelastic ratio γ ∼ 10 4 which is favorable for evaporative cooling. Since M is conserved, the spin relaxation transitions must be accompanied by the transition ℓ = 0 → ℓ ′ = 2, which leads to a centrifugal barrier in the outgoing reaction channel suppressing the s-wave inelastic scattering. Because the ℓ = 2 barrier has a fixed height of 13 mK [25], this suppression occurs only at low magnetic fields where the energy difference between the initial and final Zeeman states does not exceed the barrier height [24,25]. At larger magnetic fields, the spin-changing transitions are much more efficient, as shown in Fig. 2(b). This suggests that evaporative cooling of 3 Σ molecules in a deep magnetic trap will be challenging. However, our results indicate that evaporative cooling may be possible in shallow magnetic traps (< 1 mT deep). The maximum temperature of the molecules that can be held in such a trap is determined by the parameter where T is the temperature, k B is the Boltzmann's constant, and µ = 2µ 0 for 3 Σ molecules.
Spin relaxation in collisions of Σ-state molecules in fully spin-stretched states can only occur through coupling to the rotationally excited states assisted by the spin-spin interaction [20]. This is a two-step mechanism illustrated schematically in Fig. 1(b). The spin-spin interaction (denoted by V SS in Fig. 1 Figure   1(b) illustrates that at certain magnetic fields, the energy in the incoming collision channel becomes degenerate with that of a quasibound state supported by the uppermost curve.
This results in the formation of a long-lived complex in which one of the molecules is in the N = 2 rotationally excited state [25]. The lifetime of the complex is determined by the strength of the couplings induced by the spin-rotation interaction and the intermolecular potential as well as the presence of inelastic loss processes. In particular, the peaks shown in Fig. 3 are relatively broad because inelastic spin relaxation leads to suppression of the S-matrix poles [54].
The magnetic field dependence of the scattering cross sections for collisions of molecules in the absolute ground high-field-seeking state |M S A = −1, M S B = −1 is dramatically different (see Figure 4). The s-wave elastic scattering cross section displays a manifold of resonance peaks which we attribute to the combined action of the interaction potential and the spinspin interaction (1). As a result of an interplay between these interactions, avoided crossings occur between the incoming scattering state |M S A = −1, M S B = −1 and the quasibound states of the O 2 -O 2 complex shown in Fig. 1(b) leading to the resonant variation of the elastic cross section. This suggests that the resonances depicted in Fig. 4 are similar to the magnetic Feshbach resonances in collisions of the alkali metal atoms [12] and that the s-wave scattering length in an ultracold gas of 3 Σ molecules can be efficiently manipulated with magnetic fields.
In order to elucidate the possibility for evaporative cooling of 3 Σ molecules, it is useful to analyze the dependence of the spin relaxation rates on the rotational constant [20,55].
In ∆α where ∆α = α || − α ⊥ is the polarizability anisotropy, χ is the angle between the molecular axis and the laser polarization vector (which we assume to be parallel to B), and ǫ 0 is the amplitude of the laser field. Figure 6 shows the energy levels of 17 O 2 as functions of the laser field strength (we transform away the χ-independent terms by properly choosing the zero of energy). The effective rotational constant defined as the splitting between the closest field-dressed Zeeman levels in the N = 0 and N = 2 manifolds, increases by a factor of 1.6 as the laser field intensity is varied from zero to 5 × 10 11 W/cm 2 . For polar molecules, similar effects can be induced by dc and microwave laser fields [56].
As shown in Fig. 6, any substantial modification of molecular energy levels requires laser field intensities of order 10 11 W/cm 2 . The maximum cw laser field intensity currently available in the laboratory is about 2.5 × 10 5 W/cm 2 for a 0.5 mm size sample. However, if the trapped cloud is compressed to one-tenth its original size, the maximum attainable laser intensity increases by a factor of 100, so the required field strengths of 10 11 W/cm 2 are well within reach for microscopic clouds (10 −3 mm). Samples of such size can be produced and stored in magnetic microtraps recently demonstrated for cold polar molecules [60].
An alternative solution is to use pulsed lasers, whose intensity is typically much stronger than that of cw lasers. However, the time-independent description used in this work can only be applied if the duration of the laser pulse is long compared to the collision time.
This condition is necessary to ensure the validity of the quasistatic approximation for the molecule-field interaction [57]. These resonances can be used to create tetra-atomic molecules via magnetoassociation [19], manipulate molecule-molecule interactions in optical lattices [7], control chemical reactions of polyatomic molecules [4], or facilitate evaporative cooling of molecules in an optical dipole trap [61]. Our results show that the number of Feshbach resonances for molecule-molecule scattering in the absolute ground state (∼ 100 T −1 ) is larger by a factor of 10 than the density of the resonances for the low-field-seeking states.

APPENDIX A: MATRIX ELEMENTS OF THE INTERACTION POTENTIAL
In principle, the radial coefficients V λ A λ B λ (R) can be obtained by inverting Eq. (15) using the orthonormality properties of the spherical harmonics. This procedure requires the evaluation of six-dimensional integrals over the angles (R,r A ,r B ), which poses a difficult practical problem. A more convenient method of evaluating the expansion coefficients in Eq. (15) involves a transformation to the body-fixed (BF) frame. The z-axis of the BF frame coincides with the vector R [40,41]. The tensor product in Eq. (16) is invariant under rotations of the coordinate system [45,46], so we can substituteR =0 into Eq. (15) and use Eq. (16) to obtain where the primes indicate the BF angles. In deriving Eq. (A1) we have used the fact that Y λm λ (0) = [(2λ + 1)/4π] 1/2 . We define the BF angular basis functions [31,40,41] to expand the interaction potential as follows We note that the radial coefficients V λ A λ B λ (R) in Eqs. (A3) and (15) are exactly the same.
This is not the case for the basis functions expressed in different coordinate systems.
The BF basis functions span a coordinate subspace defined by 4 angles,r ′ A = (θ A , ϕ A ) and r ′ B = (θ B , ϕ B ). However, the interaction potential depends on the dihedral angle ϕ = ϕ A −ϕ B and not on the azimuthal angles separately. By separating the sum in Eq. (A2) into three contributions from m > 0, m = 0, and m < 0, respectively, and changing the sign of the summation indexes with m < 0, we obtain where Θ λ B ,−m (θ B ) are the normalized associated Legendre polynomials. Here, we have used the relation (−) λ A +λ B +λ = +1 that holds for two identical homonuclear molecules.
In practice, it is more convenient to use the basis functions that depend on the angles θ A , θ B , and ϕ. Even though the right-hand side of Eq. (A4) depends on the three angles, the functions on the left-hand side are defined in a four-angle space. Therefore, we need to renormalize Eq. (A4) to exclude the integration over one extra angle. We define a reduced orthonormal basis set of functions, which span the three-dimensional subspace of angles We emphasize that unlike Eq. (A4), this expression defines an orthonormal basis in the subspace of three independent angular variables. This suggests the following expansion We note that the coefficientsṼ λ A λ B λ (R) are not the same as in Eq. (15). Multiplying Eq.
and integrating over all angles, we find Comparing the expansion (A6) with Eq. (A3) and using Eq. (A5), we obtain the relation between the original coefficients defined by Eqs. (15), (A1), (A3) and those given by Eq. (A7) To summarize, the expansion coefficients V λ A λ B λ (R) can be obtained by (i) evaluating the reduced basis functions (A5); (ii) integrating the interaction potential with these functions following Eq. (A7), and (iii) using Eq. (A8) to obtain the original expansion coefficients (15).

TICAL MOLECULES
The differential cross section (DCS) for collisions of indistinguishable molecules can be expressed in terms of the symmetrized scattering amplitudeq η The integral cross section can be obtained by integrating the DCS over the coordinates of the final flux and averaging over all possible directions of the initial flux.
Since the scattered molecules in the same internal state are not distinguishable after the collision, the integration overR in Eq. (B1) has to be restricted over half-space [25,49] for the final states satisfying γ ′ A = γ ′ B . Thus, Eq. (B1) can be written as whereR i = (θ i , φ i ),R = (θ, φ), and we have defined θ max to be π/2 if γ ′ A = γ ′ B and π otherwise. The integration in Eq. (B2) can be performed trivially for θ max = π to yield [40,41] All that remains is to consider the special case γ ′ A = γ ′ B . Expanding the square modulus in Eq. (B2) yields The first integral is the usual normalization integral of two spherical harmonics. The second integral can be separated in two parts [45] 2π 0 dφ π/2 The one-dimensional integral over φ is equal to , so that the intergal (B5) becomes Since for identical molecules ℓ ′ 1 and ℓ ′ 2 are of the same parity (either even or odd) [62], and the index m ′ ℓ 1 is fixed, the associated Legendre polynomials in Eq. (B6) are both even or odd functions of cos θ [45,46]. Therefore, their product is always an even function of cos θ.
Since in this case the substitution θ → π − θ leaves the integrand in Eq. (B6) unchanged, we can halve the θ integration range and obtain