format: text

text:
# References

RCNP–Th02021

FUT–02–01

December 2002

Polarized light-flavor antiquarks from Drell-Yan processes of $`h+\stackrel{}{N}\stackrel{}{\mathrm{}^\pm }+\mathrm{}^{}+X`$
H. Kitagawa222E-mail: khisashi@post.kek.jp

Institute of Particle and Nuclear Studies, High Energy Accelerator Research Organization (KEK),
Tsukuba, Ibaragi 305–0801, Japan

Y. Sakemi999E-mail: sakemi@rcnp.osaka-u.ac.jp

Research Center for Nuclear Physics,
Osaka University, Ibaraki, Osaka 567–0047, Japan
and

T. Yamanishi333E-mail: yamanisi@ccmails.fukui-ut.ac.jp

Department of Management Science,
Fukui University of Technology, Gakuen, Fukui 910–8505, Japan

Abstract

We propose a formula to determine the first moment of difference between the polarized $`\overline{u}`$\- and $`\overline{d}`$-quarks in the nucleon, i.e. $`\mathrm{\Delta }\overline{u}\mathrm{\Delta }\overline{d}`$ from the Drell-Yan processes in collisions of unpolarized hadrons with longitudinally polarized nucleons by measuring outgoing lepton helicities. As coefficients in the differential cross section depend on the $`u`$\- and $`d`$-quark numbers in the unpolarized hadron beam, the difference $`\mathrm{\Delta }\overline{u}\mathrm{\Delta }\overline{d}`$ can be independently tested by changing the hadron beam. Moreover, a formula for estimating the $`K`$-factor in Drell-Yan processes is also suggested.

PACS number(s): 13.88.+e, 13.85.Ni, 14.65.Bt

Several experiments have been devoted to investigate the nucleon structures. In the large momentum transfer region, the nucleon (generally hadron) can be regarded as composed of almost free point-like constituents with spin $`1/2`$, $`i.e.`$ partons. It seems quite natural to identify partons as (almost) massless quarks. Then, how do quarks construct the nucleon? Inclusive deep inelastic scatterings (inclusive DIS) off lepton beams on nucleon targets have revealed some combinations of valence and sea quark distributions in the nucleon. Traditionally, the light-antiquark distributions, $`\overline{u}(x)`$ and $`\overline{d}(x)`$, have been taken to be flavor symmetric in phenomenological analyses for structure functions of nucleons from a point of view, which the strong interaction is independent of the quark flavor for light quark-pair creations from gluons. However, the asymmetry of $`\overline{u}(x)`$ and $`\overline{d}(x)`$ has been predicted , and in 1991 the NM Collaboration (NMC) at CERN has reported on the violation of the Gottfried sum rule from the unpolarized structure functions of the proton and the neutron. This implied a significant flavor asymmetry for the unpolarized light-antiquark distributions, $`\overline{u}(x)\overline{d}(x)`$. Furthermore, this asymmetry was confirmed independently by the NA51 at CERN and the E866 at FNAL through the Drell-Yan (DY) processes with proton beam and proton and deuteron targets at rapidity $`y=0`$ and for a large $`x_F`$ ($`i.e.`$ small $`x`$) region. At present, from these experimental results several approaches such as chiral quark model, Skyrme model, Pauli blocking effects, etc., are proposed to understand light-flavor antiquark asymmetry. However, the discussions are still under going. Since the strong interaction does not depend on the quark flavor for perturbative QCD, the mechanism for its flavor asymmetry may be addressed to nonperturbative QCD effect. Therefore, the studies on its mechanism are related to disclose how quarks build up the hadron, and are a challenging subject in particle and nuclear physics.

We also come to mind a question whether the longitudinally polarized light-antiquark distributions have this asymmetry or not. Unfortunately, we have no idea to estimate the first moment of $`\mathrm{\Delta }\overline{u}\mathrm{\Delta }\overline{d}`$ such as the Gottfried sum rule for the unpolarized light-antiquark flavor asymmetry. It makes the analysis on the light-flavor antiquark polarization difficult. The experimental data for the longitudinally polarized semi–inclusive DIS - have been reanalyzed recently for $`\mathrm{\Delta }\overline{u}(x)`$ and $`\mathrm{\Delta }\overline{d}(x)`$ . The result implied the asymmetry between $`\mathrm{\Delta }\overline{u}(x)`$ and $`\mathrm{\Delta }\overline{d}(x)`$, but had ambiguities because of insufficient information on the fragmentation functions and of statistical error. At present we cannot derive any conclusion on $`\mathrm{\Delta }\overline{u}(x)`$ and $`\mathrm{\Delta }\overline{d}(x)`$ from data, though several models with this asymmetry were proposed. Thus, it is important to find a formula for the first moment of $`\mathrm{\Delta }\overline{u}\mathrm{\Delta }\overline{d}`$, and to measure it directly without ambiguites.

So far, several approches for polarized light-flavor antiquark distributions in the DY process and the weak boson production have been proposed, but these processes can only probe $`\mathrm{\Delta }\overline{u}`$ and $`\mathrm{\Delta }\overline{d}`$ in a limited kinematic region.

Here we propose a formula for the first moment of $`\mathrm{\Delta }\overline{u}\mathrm{\Delta }\overline{d}`$ from the DY process with longitudinally polarized nucleon targets and unpolarized hadron beams by measuring the helicity of one of the produced pair lepton.

Let us consider a process of $`h+\stackrel{}{N}\stackrel{}{\mathrm{}^\pm }+\mathrm{}^{}+X`$. The spin-dependent cross section is defined by

$$\frac{d^3\mathrm{\Delta }\sigma }{dx_1dx_2d\mathrm{cos}\theta _\mu }\frac{d^3\sigma _{++}}{dx_1dx_2d\mathrm{cos}\theta _\mu }\frac{d^3\sigma _+}{dx_1dx_2d\mathrm{cos}\theta _\mu },$$
(1)
where the subscript ($`+`$) means the target nucleon helicity and that of one of the outgoing leptons are parallel and antiparallel, respectively. This cross section is expressed as

$`{\displaystyle \frac{d^3\mathrm{\Delta }\sigma }{dx_1dx_2d\mathrm{cos}\theta _\mu }}`$ $`=`$ $`K_{pol}{\displaystyle \frac{d\mathrm{\Delta }\widehat{\sigma }}{d\mathrm{cos}\theta _\mu }}{\displaystyle \underset{i}{}}e_i^2\left[\overline{q}_i^h(x_1,Q^2)\mathrm{\Delta }q_i^N(x_2,Q^2)+q_i^h(x_1,Q^2)\mathrm{\Delta }\overline{q}_i^N(x_2,Q^2)\right]`$ (2)
$``$ $`K_{pol}{\displaystyle \frac{d\mathrm{\Delta }\widehat{\sigma }}{d\mathrm{cos}\theta _\mu }}\mathrm{\Delta }P^{hN}(x_1,x_2,Q^2),`$
with $`x_1x_{beam}`$ and $`x_2x_{target}`$ in the leading order (LO) of QCD. $`K_{pol}`$, $`\mathrm{\Delta }\stackrel{()}{q_i^N}`$ and $`\stackrel{()}{q_i^h}`$ in eq.(2) are the $`K`$-factor of this DY process, the polarized quark (antiquark) and the unpolarized quark (antiquark) distributions with flavor $`i`$, respectively. Also $`\theta _\mu `$ is the lepton production angle in the center-of-mass frame of $`hN`$ collisions. In the LO QCD, the differential cross sections of the subprocess in eq.(2) is given by

$`{\displaystyle \frac{d\mathrm{\Delta }\widehat{\sigma }}{d\mathrm{cos}\theta _\mu }}=\mathrm{\Delta }f(x_1,x_2,Q^2,\mathrm{cos}\theta _\mu ),`$ (3)
$`\mathrm{\Delta }f(x_1,x_2,Q^2,\mathrm{cos}\theta _\mu )={\displaystyle \frac{\pi \alpha ^2}{Q^4}}{\displaystyle \frac{8x_1^2x_2^2p_N^2}{\left\{(x_1+x_2)(x_1x_2)\mathrm{cos}\theta _\mu \right\}^3}}\left\{x_1^2(1\mathrm{cos}\theta _\mu )^2x_2^2(1+\mathrm{cos}\theta _\mu )^2\right\},`$
with $`Q^2=x_1x_2s`$ and $`p_N`$ being the nucleon momentum. Here the initial $``$ and $`+`$ signs refer to $`\mathrm{\Delta }q\overline{q}\stackrel{}{\mathrm{}}^+\mathrm{}^{}`$ and $`\mathrm{\Delta }\overline{q}q\stackrel{}{\mathrm{}}^+\mathrm{}^{}`$ for measuring the positive lepton helicity, respectively. For negative lepton, $``$ and $`+`$ refer to $`\mathrm{\Delta }\overline{q}q\mathrm{}^+\stackrel{}{\mathrm{}}^{}`$ and $`\mathrm{\Delta }q\overline{q}\mathrm{}^+\stackrel{}{\mathrm{}}^{}`$, respectively.

When the target is the longitudinally polarized proton target, $`\mathrm{\Delta }P^{hN}(x_1,x_2,Q^2)`$ in eq.(2) can be written by

$`\mathrm{\Delta }P^{hp}(x_1,x_2,Q^2)`$ $`=`$ $`{\displaystyle \frac{4}{9}}\left[\overline{u}^h(x_1,Q^2)\mathrm{\Delta }u^p(x_2,Q^2)u^h(x_1,Q^2)\mathrm{\Delta }\overline{u}^p(x_2,Q^2)\right]`$ (4)
$`+{\displaystyle \frac{1}{9}}\left[\overline{d}^h(x_1,Q^2)\mathrm{\Delta }d^p(x_2,Q^2)d^h(x_1,Q^2)\mathrm{\Delta }\overline{d}^p(x_2,Q^2)\right]`$
$`+{\displaystyle \frac{1}{9}}\left[\overline{s}^h(x_1,Q^2)\mathrm{\Delta }s^p(x_2,Q^2)s^h(x_1,Q^2)\mathrm{\Delta }\overline{s}^p(x_1,Q^2)\right]`$
$`+(\mathrm{contributions}\mathrm{from}\mathrm{heavy}\mathrm{quark}\mathrm{distributions}).`$
Similary, $`\mathrm{\Delta }P^{hn}(x_1,x_2,Q^2)`$ for the neutron target can be also obtained. Assuming the isospin symmetry for the target nucleon, one has an interesting equation such as,

$`{\displaystyle \frac{d^3\mathrm{\Delta }\sigma ^{hp}}{dx_1dx_2d\mathrm{cos}\theta _\mu }}{\displaystyle \frac{d^3\mathrm{\Delta }\sigma ^{hn}}{dx_1dx_2d\mathrm{cos}\theta _\mu }}=K_{pol}{\displaystyle \frac{d\mathrm{\Delta }\widehat{\sigma }}{d\mathrm{cos}\theta _\mu }}\left\{\mathrm{\Delta }P^{hp}(x_1,x_2,Q^2)\mathrm{\Delta }P^{hn}(x_1,x_2,Q^2)\right\}`$
$`=K_{pol}\mathrm{\Delta }f(x_1,x_2,Q^2,\mathrm{cos}\theta _\mu )`$
$`\times [{\displaystyle \frac{1}{9}}[\mathrm{\Delta }u_v^p(x_2,Q^2)\mathrm{\Delta }d_v^p(x_2,Q^2)+2\{\mathrm{\Delta }\overline{u}^p(x_2,Q^2)\mathrm{\Delta }\overline{d}^p(x_2,Q^2)\}]\{4\overline{u}^h(x_1,Q^2)\overline{d}^h(x_1,Q^2)\}`$
$`{\displaystyle \frac{1}{9}}\{\mathrm{\Delta }\overline{u}^p(x_2,Q^2)\mathrm{\Delta }\overline{d}^p(x_2,Q^2)\}[4u_v^h(x_1,Q^2)d_v^h(x_1,Q^2)+2\{4\overline{u}^h(x_1,Q^2)\overline{d}^h(x_1,Q^2)\}{\displaystyle \frac{}{}}]],`$
(5)
where the initial $``$ and $`+`$ signs in the right-hand side of above equation correspond to the measurement of the positively charged lepton and negatively one, respectively. Integrating over $`x_1`$ and $`x_2`$ on eq.(5), we obtain the following relation

$`{\displaystyle _0^1}{\displaystyle _0^1}{\displaystyle \frac{d^3\mathrm{\Delta }\sigma ^{hp}/dx_1dx_2d\mathrm{cos}\theta _\mu d^3\mathrm{\Delta }\sigma ^{hn}/dx_1dx_2d\mathrm{cos}\theta _\mu }{K_{pol}\mathrm{\Delta }f(x_1,x_2,Q^2,\mathrm{cos}\theta _\mu )}}𝑑x_1𝑑x_2`$
$`={\displaystyle \frac{1}{9}}\left[\mathrm{\Delta }u_v^p(Q^2)\mathrm{\Delta }d_v^p(Q^2)+2\left\{\mathrm{\Delta }\overline{u}^p(Q^2)\mathrm{\Delta }\overline{d}^p(Q^2)\right\}\right]\left\{4\overline{u}^h(Q^2)\overline{d}^h(Q^2)\right\}`$
$`{\displaystyle \frac{1}{9}}\left\{\mathrm{\Delta }\overline{u}^p(Q^2)\mathrm{\Delta }\overline{d}^p(Q^2)\right\}\left[4u_v^hd_v^h+2\left\{4\overline{u}^h(Q^2)\overline{d}^h(Q^2)\right\}\right]`$
$`={\displaystyle \frac{1}{9}}\left|{\displaystyle \frac{g_A}{g_V}}\right|\left\{4\overline{u}^h(Q^2)\overline{d}^h(Q^2)\right\}{\displaystyle \frac{1}{9}}\left\{\mathrm{\Delta }\overline{u}^p(Q^2)\mathrm{\Delta }\overline{d}^p(Q^2)\right\}\left[4u_v^hd_v^h+2\left\{4\overline{u}^h(Q^2)\overline{d}^h(Q^2)\right\}\right]`$
(6)
for the measurement of the $`\mathrm{}^{}`$ helicity. In the eq.(6), $`g_A`$ and $`g_V`$ are the nucleon axial and vector coupling constants, respectively, and $`u_v^h`$ and $`d_v^h`$ are numbers of the valence $`u`$\- and $`d`$-quarks in the beam particle $`h`$, respectivley. Here, we drop the label $`Q^2`$, since the valence quark numbers in the hadron are independent of $`Q^2`$. Therefore, appraising $`4\overline{u}^h(Q^2)\overline{d}^h(Q^2)`$ in the unpolarized beam hadron $`h`$ and the $`K`$-factor of this polarized DY process, we get information on the first moment of $`\mathrm{\Delta }\overline{u}^p\mathrm{\Delta }\overline{d}^p`$ from the cross sections for $`h+\{\stackrel{}{p}\text{ and }\stackrel{}{n}\}\stackrel{}{\mathrm{}^\pm }+\mathrm{}^{}+X`$ with wide ranges of both $`x_1`$ and $`x_2`$.

If the polarized light-flavor antiquark distribution is symmetric, i.e. $`\mathrm{\Delta }\overline{u}^p\mathrm{\Delta }\overline{d}^p=0`$, the right-hand side of eq.(6) reduces to $`1/9|g_A/g_V|\{4\overline{u}^h(Q^2)\overline{d}^h(Q^2)\}`$. For instance, choosing the proton as the unpolarized hadron beam $`h`$ and taking $`x_{1min}=10^5`$ and $`Q^2=4`$GeV2, it becomes $`1/91.2670.14\times 10^2`$$`=19.7`$ with $`0.141\times 10^2`$ or $`0.139\times 10^2`$ for $`4\overline{u}^p(Q^2)\overline{d}^p(Q^2)`$ by the parametrization of GRV98LO or MRST98LO, respectively. Also, the difference between the spin-dependent differential cross sections of $`pp`$ and $`pn`$ collisions as a function of $`x_2`$ for several $`x_1`$ values is shown in fig.1 by taking $`K_{pol}=1.8`$. We use AAC with $`\mathrm{\Delta }\overline{u}=\mathrm{\Delta }\overline{d}`$ and GRV98LO parametrizations as polarized and unpolarized distribution functions, respectively. Therefore, we can conclude that the behavior of $`\mathrm{\Delta }\overline{u}`$ and $`\mathrm{\Delta }\overline{d}`$ in the nucleon is asymmetric if we find a discrepancy between the measured values and above predicted ones.

For other unpolarized hadron beams, for example charged pions, kaons and so on, it is not difficult to estimate $`4\overline{u}^h(Q^2)\overline{d}^h(Q^2)`$ in eq.(6) from experiments. The value of $`4\overline{u}^h(Q^2)\overline{d}^h(Q^2)`$ in the unpolarized beam hadron $`h`$ can be obtained from data in the same way as above. With unpolarized nucleon targets, the same procedure as eq.(6) leads to

$`{\displaystyle _0^1}{\displaystyle _0^1}{\displaystyle \frac{d^3\sigma ^{hp}/dx_1dx_2d\mathrm{cos}\theta _\mu d^3\sigma ^{hn}/dx_1dx_2d\mathrm{cos}\theta _\mu }{K_{unpol}d\widehat{\sigma }/d\mathrm{cos}\theta _\mu }}𝑑x_1𝑑x_2`$
$`={\displaystyle \frac{1}{9}}\left\{\overline{u}^p(Q^2)\overline{d}^p(Q^2)\right\}(4u_v^hd_v^h)+{\displaystyle \frac{1}{9}}\left[u_v^pd_v^p+2\left\{\overline{u}^p(Q^2)\overline{d}^p(Q^2)\right\}\right]\left\{4\overline{u}^h(Q^2)\overline{d}^h(Q^2)\right\},`$
(7)
where $`K_{unpol}`$ is the $`K`$-factor of the unpolarized DY process, and the differential cross sections of the subprocess is written as

$$\frac{d\widehat{\sigma }}{d\mathrm{cos}\theta _\mu }=\frac{\pi \alpha ^2}{Q^4}\frac{8x_1^2x_2^2p_N^2}{\left\{(x_1+x_2)(x_1x_2)\mathrm{cos}\theta _\mu \right\}^3}\left\{x_1^2(1\mathrm{cos}\theta _\mu )^2+x_2^2(1+\mathrm{cos}\theta _\mu )^2\right\}.$$

Since the value of $`\overline{u}^p(Q^2)\overline{d}^p(Q^2)`$ in eq.(7) was obtained by the NMC and other experiments and also has been studied intensively, one can extract $`4\overline{u}^h(Q^2)\overline{d}^h(Q^2)`$ from the combination of the differential cross sections, $`d\sigma ^{hp}`$ and $`d\sigma ^{hn}`$, for the unpolarized DY processes. Using the NMC result, namely $`\overline{u}^p(Q^2)\overline{d}^p(Q^2)=0.147\pm 0.039`$ at $`Q^2=4`$GeV2, and choosing the proton as the hadron $`h`$, we obtain

$`4\overline{u}^h(Q^2)\overline{d}^h(Q^2)`$
$`=\left[{\displaystyle _0^1}{\displaystyle _0^1}{\displaystyle \frac{d^3\sigma ^{hp}/dx_1dx_2d\mathrm{cos}\theta _\mu d^3\sigma ^{hn}/dx_1dx_2d\mathrm{cos}\theta _\mu }{K_{unpol}d\widehat{\sigma }/d\mathrm{cos}\theta _\mu }}𝑑x_1𝑑x_2+0.114\right]/0.0784.`$ (8)
For $`K_{pol}`$ and $`K_{unpol}`$ appeared in eqs.(6)-(8), in general, it is important to find an equation where the $`K`$-factor cancels out in order to extract physical quantities from DY processes. However, we do not have enough information derived from such equation. In order to get absolute values of the physical quantities induced from the DY processes, the $`K`$-factor must be estimated exactly. Using unpolarized antihadron beams together with unpolarized hadron beams, the mean value of $`K_{pol}`$ in terms of $`x`$ for this polarized DY process is given as follows

$`{\displaystyle _0^1}{\displaystyle _0^1}{\displaystyle \frac{\left\{\frac{d^3\mathrm{\Delta }\sigma ^{\overline{h}p}}{dx_1dx_2d\mathrm{cos}\theta _\mu }\frac{d^3\mathrm{\Delta }\sigma ^{\overline{h}n}}{dx_1dx_2d\mathrm{cos}\theta _\mu }\right\}\left\{\frac{d^3\mathrm{\Delta }\sigma ^{hp}}{dx_1dx_2d\mathrm{cos}\theta _\mu }\frac{d^3\mathrm{\Delta }\sigma ^{hn}}{dx_1dx_2d\mathrm{cos}\theta _\mu }\right\}}{\mathrm{\Delta }f(x_1,x_2,Q^2,\mathrm{cos}\theta _\mu )}}𝑑x_1𝑑x_2`$
$`=K_{pol}{\displaystyle \frac{1}{9}}\left|{\displaystyle \frac{g_A}{g_V}}\right|\left(4u_v^hd_v^h\right),`$ (9)
with $``$ being similar to eq.(5), where we assume the reflection symmetry along the V-spin axis, the isospin symmetry and the charge conjugation invariance for the unpolarized beam hadron with the polarized nucleon target. For $`K_{unpol}`$ of the unpolarized DY process, we also have

$`{\displaystyle _0^1}{\displaystyle _0^1}{\displaystyle \frac{\left\{\frac{d^3\sigma ^{\overline{h}p}}{dx_1dx_2d\mathrm{cos}\theta _\mu }\frac{d^3\sigma ^{\overline{h}n}}{dx_1dx_2d\mathrm{cos}\theta _\mu }\right\}\left\{\frac{d^3\sigma ^{hp}}{dx_1dx_2d\mathrm{cos}\theta _\mu }\frac{d^3\sigma ^{hn}}{dx_1dx_2d\mathrm{cos}\theta _\mu }\right\}}{d\widehat{\sigma }/d\mathrm{cos}\theta _\mu }}𝑑x_1𝑑x_2`$
$`=K_{unpol}{\displaystyle \frac{1}{9}}\left(u_v^pd_v^p\right)\left(4u_v^hd_v^h\right).`$ (10)
Since each term of the right-hand side in eqs.(9) and (10) is constant, $`K_{pol}`$ and $`K_{unpol}`$ can be evaluated from the relevant differential cross sections.

Thus, measuring the $`K`$-factors by using unpolarized hadron and antihadron beams, and evaluating $`4\overline{u}^h(Q^2)\overline{d}^h(Q^2)`$ from the unpolarized DY experiment, the polarized DY process for an unpolarized hadron beam and a polarized nucleon target allows to provide the first moment of $`\mathrm{\Delta }\overline{u}^p(Q^2)\mathrm{\Delta }\overline{d}^p(Q^2)`$ by measuring the helicity of one of the produced pair lepton. Note that it does not require measurements with high precision for large $`x_1`$ and $`x_2`$ region in order to determine the value of $`\mathrm{\Delta }\overline{u}^p\mathrm{\Delta }\overline{d}^p`$, though differential cross sections for the proton and the neutron targets are combined. Because each differential cross section for large $`x_1`$ and $`x_2`$ is quite small as shown in fig.1, a contribution from this region to the integral is small.

In summary, we have proposed a formula for the difference between the polarized light-flavor antiquark density, $`\mathrm{\Delta }\overline{u}^p\mathrm{\Delta }\overline{d}^p`$, from DY processes. It is given by a combination of the cross sections with an unpolarized hadron beam and a longitudinally polarized proton target by measuring one of the produced lepton helicity and that with a neutron target. Then, the formula is described in terms of the neutron $`\beta `$-decay constant and the difference between $`\mathrm{\Delta }\overline{u}^p`$ and $`\mathrm{\Delta }\overline{d}^p`$. As coefficients of these terms depend on $`u`$\- and $`d`$-quark numbers in the unpolarized beam hadron, we can independently get information on the behavior of $`\mathrm{\Delta }\overline{u}^p`$ and $`\mathrm{\Delta }\overline{d}^p`$ by changing the beam hadron.

Recently the DY process in the next-to-leading order (NLO) of QCD has been discussed and compared with that in the LO, though it is the DY process for the longitudinally polarized nucleon and the longitudinally polarized nucleon collisions. The relevant differential cross section has an additional term including contributions from gluons. However, in the difference $`d\mathrm{\Delta }\sigma ^{hp}d\mathrm{\Delta }\sigma ^{hn}`$ the contibutions from gluons will be cancel. Accordingly, eq.(6) is expected to be also kept in the NLO.

Now experiments to solve the nucleon spin problems start at Relativestic Heavy Ion Collider (RHIC) in BNL with colliding polarized protons at high energy, $`\sqrt{s}=200`$GeV. Also, Japan Hadron Facility (JHF) is under construction. It provides 50GeV high intensity proton and antiproton beams, and also can produce high intensity charged pi/K meson beams. One expects that future experiments for processes proposed here will be carried at RHIC and/or JHF with the detector measuring the helicity of high energy leptons over wide ranges of $`x_1`$ and $`x_2`$. It is certainly difficult to measure the high energy lepton helicity. However, the measurement of the helicity of the $`\mu ^+`$ produced in charged current interactions has already been carried out at high-energy antineutrino experiments at CERN-SPS. The polarimeter used there composed of the marbles for stopping the muon, the scintillators for detecting the positron from muon decay and the proportional drift tubes for observing the muon track. This makes the use of the fact that high energy positron is preferentially emitted in the direction of the muon spin because of the V-A intereaction in muon decay. Accordingly, it seems possible to do experiments on the processes discussed here. We wish to get new informations on $`\mathrm{\Delta }\overline{u}^p\mathrm{\Delta }\overline{d}^p`$ at RHIC and/or JHF.

One of us (T. Y.) thanks T. Morii for enlightening comments, and K. Kobayakawa for careful reading of this manuscript.


text:
# Ultracold atomic collisions in tight harmonic traps: Quantum defect model and application to metastable helium atoms

## I Introduction

An understanding of ultracold collision processes between neutral atoms is crucial to the design and operation of atom traps, and to the development of novel quantum processes using trapped atoms Julienne . The elastic collision rate must be high enough to produce efficient thermalization during the evaporative cooling phase of magnetostatic trapping, whereas the inelastic collision rate must be small, because such collisions can generate energetic atoms and change the atomic state, hence destroying the trapping conditions and producing trap loss. Elastic collisions are also important in studies of Bose-Einstein condensates where they determine the mean field of the condensate BEC . Ultracold collisions are usually studied under weak trapping conditions in which the confining inhomogeneous magnetic field is either ignored or is assumed to have a parabolic or harmonic spatial variation of sufficiently low frequency (typically $`10^2`$ Hz) that it can be treated as uniform during the collision. However, recent interest in phenomena such as quantum phase transitions of 87Rb atoms confined in three-dimensional optical lattices Greiner , and far off-resonance three-dimensional optical lattices to store metastable argon atoms Muller  or to implement quantum logic gates and create highly entangled quantum states qcomputing , involve conditions where the trapping frequency is typically $`10^5`$ to $`10^6`$ Hz and the tight trapping environment is expected to significantly modify the properties of the colliding system.

In several existing calculations for tightly confined neutral atoms the exact interatomic potential is replaced by the regularised $`\delta `$\- function pseudopotential

$$V_\delta (𝒓)=\frac{2\pi \mathrm{}^2}{M}a\delta (𝒓)\frac{}{r}r,$$
(1)
where $`a`$ is the scattering length and $`M`$ is the reduced mass of the system. This potential reproduces the s-wave phase shifts in the Wigner threshold regime and also the correct asymptotic behavior of the wavefunction at large $`r`$. This enables an analytical solution for the case of a spherically symmetric harmonic trap to be obtained Busch1998 , with the energy eigenvalues determined from the condition

$$\frac{a}{\xi }=f(E)\frac{1}{2}\mathrm{tan}\left(\frac{\pi E}{2\mathrm{}\omega }+\frac{\pi }{4}\right)\frac{\mathrm{\Gamma }(\frac{E}{2\mathrm{}\omega }+\frac{1}{4})}{\mathrm{\Gamma }(\frac{E}{2\mathrm{}\omega }+\frac{3}{4})},$$
(2)
where $`\omega `$ is the trap frequency and $`\xi =\sqrt{\mathrm{}/M\omega }`$ is the effective range of the ground state wavefunction. The validity of this approach has been investigated by Tiesinga et al. Tiesinga2000  for the Na and Cs systems by comparing the energy eigenvalues with those computed numerically using the best available full interatomic potentials. They find that this approximation is limited to sufficiently weak traps where $`\xi a`$. For the case of two atoms interacting via a hard-sphere potential of range $`a`$, Block and Holthaus Block2002  have shown that the pseudopotential has the form (1) with $`a`$ given by (2). Recently, two groups Bolda2002 ; Blume2002  have advocated a model in which an energy-dependent effective scattering length

$$a_{\text{eff}}(E)=\frac{\mathrm{tan}\delta _0(k)}{k}$$
(3)
is introduced, where $`\mathrm{}k=\sqrt{2ME}`$ and $`\delta _0(k)`$ is the s-wave phase shift for the untrapped atoms scattering in the interatomic potential. The energy eigenvalue condition

$$\frac{a_{\text{eff}}(E)}{\xi }=f(E)$$
(4)
is then solved self-consistently. The procedure given in (3) and (4) reproduces the asymptotic wavefunction and the s-wave phase shifts even at energies above the Wigner threshold regime, and it is found that this self-consistent (SC) method works even when $`a_{\text{eff}}/\xi 1`$.

In this paper we analyze a system of two interacting atoms under strong harmonic confinement from the viewpoint of quantum defect theory (QDT) and derive a generalized SC model which validates (4) and extends it to non s-wave collisions. We then present two highly efficient computational methods for determining the bound state energies and eigenfunctions for the trapped atoms. The first method, based on QDT, is an adaptation of a technique developed by one of the authors (GP) for highly excited states in a modified Coulomb potential. The second method is a direct numerical solution of the radial Schrödinger equation using a discrete variable representation (DVR) of the kinetic energy operator and a scaled radial coordinate grid. Both methods are applicable to a system of two arbitrary ultracold neutral atoms interacting in a harmonic trap and, as an example, we consider the case of spin-polarized metastable helium tightly confined in a spherically symmetric harmonic trap as this system has not been studied before under these conditions.

Interest in collision processes in metastable 23S helium (denoted by He) has been generated by the quest to attain Bose-Einstein condensation in this system Santos2001 ; Robert2001  and, subsequently, to understand these condensates Sirjean2002 . Such condensates are novel in that they are the first excited-state condensates and open up new fields for investigation such as those of atomic correlations and the growth kinematics of the condensate. This experimental success has depended upon the correctness of the theoretical prediction Shlyap1994 ; Fedichev1996 ; Venturi1999 ; Venturi2000 ; Leo2001  that the inelastic Penning ionization processes can be strongly suppressed through spin polarization of the He system in the magnetostatic trap.

This paper is organized as follows. In Sec. II the formalism is developed for collisions of two neutral atoms in an external three-dimensional isotropic harmonic trap, and the general nature of the energy eigenvalues and eigenfunctions of the resultant radial Schrödinger equation discussed. A quantum defect is introduced and shown to be analytic in energy. In Sec. III we formulate a generalization of the SC method, and the QDT and DVR computational methods are presented in Sec. IV. In Sec. V the QDT, DVR and SC methods are applied to ultracold metastable helium atoms under tight harmonic confinement. Results are obtained for s-wave collisions over a range of trapping frequencies and also for d-wave collisions in a 10 MHz trap. Finally, in Sec. VI, we summarize and discuss our results.

## II Two-atom collisions in a harmonic trap

Consider two atoms $`j=1,2`$ of mass $`M_j`$ and position $`𝒓_j`$ relative to the centre of the trap. For a central interatomic potential $`V(r)`$, where $`r=|𝒓|=|𝒓_1𝒓_2|`$, and an isotropic harmonic trap, the two-atom Hamiltonian is separable into centre-of-mass and relative motions. The equation for the relative motion of angular momentum $`l`$ is

$$\left[\frac{\mathrm{}^2}{2M}_r^2+\frac{1}{2}M\omega ^2r^2+V(r)\right]\psi (𝒓)=E\psi (𝒓),$$
(5)
where $`E`$ is the energy eigenvalue and the reduced mass $`M=M_1M_2/(M_1+M_2)`$. The trap potential $`V_{\text{trap}}(r)=M\omega ^2r^2/2`$ has been assumed to be independent of the atomic state, which is generally valid in far-detuned optical lattices Grimm . As the interaction is spherically symmetric, $`\psi (𝒓)`$ has the form

$$\psi (𝒓)=\frac{1}{r}F_{kl}(r)Y_{lm}(\theta ,\varphi ),$$
(6)
where $`Y_{lm}(\theta ,\varphi )`$ are spherical harmonics and $`F_{kl}(r)`$ satisfies

$$\left[\frac{\mathrm{}^2}{2M}\frac{d^2}{dr^2}+\frac{l(l+1)\mathrm{}^2}{2Mr^2}+\frac{1}{2}M\omega ^2r^2+V(r)\right]F_{kl}(r)=EF_{kl}(r).$$
(7)
It is convenient to introduce the dimensionless variables $`\rho =r/\xi `$ and $`\kappa =2E/\mathrm{}\omega =\xi ^2k^2`$ and rewrite (7) as

$$\left[\frac{d^2}{d\rho ^2}\frac{l(l+1)}{\rho ^2}+\kappa \rho ^2\frac{2V(\rho )}{\mathrm{}\omega }\right]F(\rho )=0.$$
(8)
For large values of $`\rho `$, $`V(\rho )`$ is negligible and $`F(\rho )`$ has the asymptotic form

$$F_a(\rho )=z^{l/2+1/2}e^{z/2}w(z);z=\rho ^2$$
(9)
where $`w(z)`$ is a linearly combination of the two independent solutions given by

$$w_1(z)={}_{1}{}^{}F_{1}^{}(a;c;z);w_2(z)=z^{1c}{}_{1}{}^{}F_{1}^{}(1+ac;2c;z)$$
(10)
on using the notation adopted by Luke Luke . In (10)

$$a=\frac{l}{2}+\frac{3}{4}\frac{\kappa }{4};c=l+\frac{3}{2}$$
(11)
and $`{}_{1}{}^{}F_{1}^{}(a;c;z)`$ is the confluent hypergeometric function

$${}_{1}{}^{}F_{1}^{}(a;c;z)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{\mathrm{\Gamma }(a+n)}{\mathrm{\Gamma }(a)}\frac{\mathrm{\Gamma }(c)}{\mathrm{\Gamma }(c+n)}\frac{z^n}{n!}.$$
(12)
For the case of the unperturbed oscillator the solution $`w_1(z)`$, which is regular at the origin, is bounded as $`z\mathrm{}`$ provided that $`a=n_r`$, where $`n_r`$ is a non-negative integer. The energy eigenvalues are then given by

$$E_0=\mathrm{}\omega (2n_r+l+\frac{3}{2});n_r=0,1,2,\mathrm{},$$
(13)
where $`n_r`$ denotes the number of nodes in the corresponding radial wavefunction $`F_{kl}^0(r)`$.

In the presence of collisions, the energy eigenvalues $`E`$ are no longer equal to $`E_0`$ but can be written in the form

$$E=\mathrm{}\omega (2n_r^{}+l+\frac{3}{2}).$$
(14)
Thus the effect of the collisions is to replace $`n_r`$ by

$$n_r^{}=n_r\mu $$
(15)
where $`\mu `$ will be called a quantum defect and is in general not an integer; this quantum defect has previously been introduced by Blume and Greene Blume2002 .

The present problem, characterized by a long asymptotic trap region beyond the effective range of the interatomic potential, shows similar features to the modified Coulomb problem in which a long-range attractive Coulomb potential is supplemented by a short-range interaction. This is an extremely well-known situation in connection with the analysis of atomic energy levels and spectra where it is common to specify the energy of a level with quantum numbers $`nl`$ in terms of an effective principal quantum number $`n^{}`$. The quantum defect is then defined to be the difference between $`n`$ and $`n^{}`$. In this context there has been extensive development of a quantum defect theory, see Seaton Seaton , and in the following, we demonstrate that equation (8) for the perturbed harmonic oscillator can be mathematically transformed into an equation that is identical in form to that for the modified Coulomb problem. This enables quantum defect theory to be directly applied here and valuable theoretical insight into the underlying behaviour of the energy level structure obtained.

Asymptotically, the exponentially decaying eigenfunction (9) is given by

$$w_3(z)=\psi (a;c;z)=\frac{\mathrm{\Gamma }(1c)}{\mathrm{\Gamma }(1+ac)}w_1(z)+\frac{\mathrm{\Gamma }(c1)}{\mathrm{\Gamma }(a)}w_2(z)$$
(16)
where $`a=n_r^{}`$, and as $`z\mathrm{}`$,

$$\psi (a;c;z)z^a{}_{2}{}^{}F_{0}^{}(a,1+ac;\frac{1}{z}),$$
(17)
where

$${}_{2}{}^{}F_{0}^{}(a;c;z)=\underset{n=0}{\overset{\mathrm{}}{}}\frac{\mathrm{\Gamma }(a+n)}{\mathrm{\Gamma }(a)}\frac{\mathrm{\Gamma }(c+n)}{\mathrm{\Gamma }(c)}\frac{z^n}{n!}.$$
(18)
The quantum defect $`\mu `$ is an analytic function of energy for $`E>0`$. To prove this we make the transformation

$$F_a(\rho )=z^{1/4}Y(y)$$
(19)
where $`y=\kappa z/8`$. The function $`Y(y)`$ then satisfies

$$\left[\frac{d^2}{dy^2}\frac{\lambda (\lambda +1)}{y^2}+\frac{2}{y}+ϵ\right]Y(y)=0,$$
(20)
where $`\lambda =l/21/4`$ and

$$ϵ=\frac{1}{n^2};n^{}\frac{\kappa }{4}=n_r^{}+\lambda +1=\nu \mu $$
(21)
with $`\nu =n_r+\lambda +1`$. Equation (20) has exactly the form of the Coulomb equation for a bound state $`nl`$, in which $`n`$ and $`l`$ have been replaced by $`\nu `$ and $`\lambda `$. The present problem differs from that studied by Seaton in that, here, $`\lambda `$ cannot be an integer. Consequently we introduce the two linearly independent solutions of (20)

$$Y_1(n^{},\lambda ;y)=\frac{(n^{}z)^{\lambda +1}e^{z/2}}{\mathrm{\Gamma }(2\lambda +2)}{}_{1}{}^{}F_{1}^{}(\lambda +1n^{},2\lambda +2;z)=\frac{(n^{}z)^{\lambda +1}e^{z/2}}{\mathrm{\Gamma }(2\lambda +2)}w_1$$
(22)
and

$$Y_2(n^{},\lambda ;y)=\frac{(n^{}z)^\lambda e^{z/2}}{\mathrm{\Gamma }(2\lambda )}{}_{1}{}^{}F_{1}^{}(\lambda n^{},2\lambda ;z)=\frac{(n^{})^\lambda z^{\lambda +1}e^{z/2}}{\mathrm{\Gamma }(2\lambda )}w_2,$$
(23)
where $`a=\lambda +1n^{}`$ and $`c=2\lambda +2`$ in (10). The functions $`Y_1`$ and $`Y_2`$ are identical to the functions $`y_1`$ and $`y_2`$ defined by Seaton and are analytic functions of $`ϵ`$. We then introduce the general solution

$$Y_3(n^{},\lambda ;y)=z^{\lambda +1}e^{z/2}\left[\alpha (n^{})\frac{(n^{})^{\lambda +1}}{\mathrm{\Gamma }(2\lambda +2)}w_1+\beta (n^{})\frac{(n^{})^\lambda }{\mathrm{\Gamma }(2\lambda )}w_2\right],$$
(24)
where $`\alpha (n^{})`$ and $`\beta (n^{})`$ are analytic functions of $`n^{}`$ ($`E`$), see (14) and (21). On expressing $`w_2`$ in terms of $`w_1`$ and the exponentially decaying solution $`w_3`$ in (16), the condition for (24) to decay as $`z\mathrm{}`$ is that the coefficient of $`w_1`$ must vanish. This requires that

$$\alpha (n^{})A(n^{},\lambda )B(n^{},\lambda )\beta (n^{})=0,$$
(25)
where

$$A(n^{},\lambda )=\frac{\mathrm{\Gamma }(\lambda +1+n^{})}{\mathrm{\Gamma }(n^{}\lambda )(n^{})^{2\lambda +1}}1+O(1/n^2)$$
(26)
for large values of $`n^{}`$, and

$$B(n^{},\lambda )=\frac{\mathrm{sin}[\pi (n^{}+\lambda )]}{\mathrm{sin}[\pi (n^{}\lambda )]}=(1)^l\mathrm{cot}\pi \mu .$$
(27)
In practice for the cases considered in this paper for which $`l2`$ , $`A(n^{},\lambda )`$ remains very close to unity except for the lowest trap states. Therefore from (25)–(27)

$$\alpha (n^{})\mathrm{sin}(\pi \mu )\beta (n^{})(1)^l\mathrm{cos}(\pi \mu )[1+O(1/n^2)]=0,$$
(28)
and so $`\mu `$ can be written in the form

$$\mu =a+bn^{}+cn^2+\mathrm{}.$$
(29)
In general, the interatomic potential $`V(r)`$ supports a number of bound states, $`n_b`$ say, for $`E<0`$, so that the lowest trap state $`(E>0)`$ has $`n_r=n_b`$, that is, the lowest trap state eigenfunction has $`n_b`$ nodes. If however a pseudopotential is used such as $`V_\delta (𝒓)`$ in (1), there are no bound states with $`E<0`$ and the wave function for the lowest state with $`E>0`$ has no nodes. In this case the number of nodes, $`n_r^{}`$, and the quantum defect, $`\mu ^{}`$, are defined by

$$n_r^{}=n_rn_b;\mu ^{}=\mu n_b.$$
(30)
The differences between these two types of potential and the effects of their use are discussed in detail by Peach Peach .

## III Generalization and Validation of the SC method

Beyond the effective range of the interatomic potential there is an extended region, typically $`10^2a_0r10^4a_0`$, where the effect of the trapping potential is extremely small and the asymptotic form of the wavefunction is essentially that of a free wave. Therefore in this region, to an excellent approximation $`F_{kl}(r)`$ can be written as

$$F_{kl}(r)=kr[\mathrm{cos}\delta _lj_l(kr)\mathrm{sin}\delta _ln_l(kr)],$$
(31)
where $`j_l(kr)`$ and $`n_l(kr)`$ are spherical Bessel functions, see Abram , and $`\delta _l(k)`$ is the $`l`$-wave phase shift. However in the previous section it has been shown that the correct asymptotic form of the wave function is obtained from (19) and (24), so that (31) must correspond to

$$F_{kl}(r)=CF_a(\rho )=Cz^{1/4}\left[\alpha (n^{})Y_1+\beta (n^{})Y_2\right],$$
(32)
where $`C`$ is a constant. In equation (32), the coefficients $`\alpha (n^{})`$ and $`\beta (n^{})`$ satisfy equations (25)-(27) and so this form of the solution does decay exponentially at very large separations as it must for a bound state. We can ensure that the two forms (31) and (32) are identical throughout the region by examining their behavior for small values of $`r`$. The Bessel functions $`j_l(kr)`$ and $`n_l(kr)`$ in (31) take the form

$$krj_l(kr)\frac{2^l\mathrm{\Gamma }(l+1)}{\mathrm{\Gamma }(2l+2)}(kr)^{l+1};krn_l(kr)\frac{\mathrm{\Gamma }(2l+1)}{2^l\mathrm{\Gamma }(l+1)}(kr)^l,$$
(33)
whereas from (10) and (12), for small values of $`z=r^2/\xi ^2`$,

$$z^{1/4}Y_1\frac{(n^{})^{\lambda +1}z^{\lambda +3/4}}{\mathrm{\Gamma }(2\lambda +2)};z^{1/4}Y_2\frac{(n^{})^\lambda z^{\lambda 1/4}}{\mathrm{\Gamma }(2\lambda )}.$$
(34)
Since $`\lambda =l/21/4`$, the behavior of $`F_{kl}(r)`$ as $`r0`$ in (32) is

$$F_{kl}(r)C\left[\frac{(n^{})^{\lambda +1}\alpha (n^{})}{\xi ^{l+1}\mathrm{\Gamma }(l+\frac{3}{2})}r^{l+1}+\frac{(n^{})^\lambda \beta (n^{})\xi ^l}{\mathrm{\Gamma }(l+\frac{1}{2})}r^l\right],$$
(35)
and then from (31), (33) and (35) we obtain

$$\frac{1}{(n^{})^{2\lambda +1}}\frac{\beta (n^{})}{\alpha (n^{})}=\frac{\mathrm{tan}\delta _l}{(\xi k)^{2l+1}}\frac{\mathrm{\Gamma }(l+\frac{1}{2})}{\mathrm{\Gamma }(l+\frac{3}{2})}\frac{\mathrm{\Gamma }(2l+1)\mathrm{\Gamma }(2l+2)}{[2^l\mathrm{\Gamma }(l+1)]^2}$$
(36)
where $`\xi ^2k^2=4n^{}`$. Combining this with the expression for $`\beta (n^{})/\alpha (n^{})`$ derived from equations (25)-(27), our final result is

$$\frac{\mathrm{tan}\delta _l(k)}{(\xi k)^{2l+1}}=f_l(E),$$
(37)
where $`n^{}=E/2\mathrm{}\omega `$ and

$`f_l(E)={\displaystyle \frac{1}{2^{2l+1}}}\mathrm{tan}[\pi (n^{}{\displaystyle \frac{l}{2}}+{\displaystyle \frac{1}{4}})]{\displaystyle \frac{\mathrm{\Gamma }(n^{}\frac{l}{2}+\frac{1}{4})}{\mathrm{\Gamma }(n^{}+\frac{l}{2}+\frac{3}{4})}}.`$ (38)
For $`l=0`$, (38) reduces to the result given by equations (3) and (4). Thus not only have we obtained a generalization of the SC method to the case of $`l0`$, but we have proved mathematically that relation (4) which was originally introduced empirically, is in fact rigorously true.

The behavior of (37) for $`l>0`$ warrants some discussion. For potentials $`V(r)`$ with the asymptotic form $`r^n`$ at large $`r`$, the threshold behavior at small $`k`$ of the phaseshifts $`\delta _l(k)`$ is Geltman

$$\mathrm{tan}\delta _l(k)\{\begin{array}{cc}k^{2l+1};\hfill & n>2l+3\hfill \\ k^{2l+1}\mathrm{ln}k;\hfill & n=2l+3\hfill \\ k^{n2};\hfill & n<2l+3\hfill \end{array}.$$
(39)
In our case, $`n=6`$ and the behavior as $`k0`$ is given by

$$\frac{\mathrm{tan}\delta _l(k)}{k^{2l+1}}a_l;l=0,1$$
(40)
and

$$\frac{\mathrm{tan}\delta _l(k)}{k^4}b_l;l2,$$
(41)
where $`a_l`$ and $`b_l`$ are constants.

## IV Computational methods for trapped atoms

The accurate numerical determination of bound state energies and eigenfunctions for interacting atoms in a trap requires computational techniques optimized to deal with both the small inner interatomic region and the large outer trap region. Two such methods are now presented in some detail as they will form the basis of subsequent investigations into more complex processes involving trapped atoms.

The first approach is an adaptation of a method (QDT) developed by one of us (GP) for highly excited states in a modified Coulomb potential. We consider the region $`rr_0(20a_0`$ to $`40a_0)`$ where the asymptotic solution specified in (24) is only weakly perturbed by the presence of the interatomic potential $`V(\rho )`$. We first of all generate this asymptotic solution throughout the region $`rr_0`$ and then calculate a multiplicative correction factor in order to obtain an accurate solution of (20). Equation (20) is written as

$$\left[\frac{d^2}{dy^2}g(y)\right]Y(y)=0;g(y)=\frac{\lambda (\lambda +1)}{y^2}\frac{2}{y}ϵ,$$
(42)
and the region is divided up into $`N`$ ranges, $`y=y_n,n=0,1,2,\mathrm{}N`$, say. $`y_N`$ is chosen so that for $`yy_N`$ the decaying solution $`Y_3(y)`$ can be evaluated by using its asymptotic form, see (17) and (18). Within each range $`y_{n1}yy_n`$ the following expansions are made

$$g(y)=\underset{m=0}{\overset{\mathrm{}}{}}\frac{(yy_n)^m}{m!}g^{(m)}(y_n);Y(y)=\underset{m=0}{\overset{M_0}{}}a_m(yy_n)^m$$
(43)
and substituted into (42). The coefficients $`a_m`$ can then be obtained by equating the coefficients of powers of $`(yy_m)`$ to zero, and the values of $`Y(y_{n1})`$ and $`Y^{}(y_{n1})`$ obtained from the range $`y_{n1}yy_n`$ provide the input values for the solution in the next range $`y_{n2}yy_{n1}`$, etc. The choice of $`M_0`$ in each interval clearly depends on its length, but it is found that the accuracy of the solution $`Y(y)`$ is very insensitive to the precise choice of $`N`$ and $`M_0`$. One of the $`y_n`$ is chosen to be at the outer turning point given by $`g(y)=0,y=y_a`$ where

$$y_a=n^2\left[1+\sqrt{1\frac{\lambda (\lambda +1)}{n^2}}\right]$$
(44)
and $`y_0`$ and $`\rho _a`$ are defined by, see (20) and (21),

$$y_0=n^{}(\rho _0)^2/2;y_a=n^{}(\rho _a)^2/2,$$
(45)
where $`\rho _0r_0/\xi `$ and $`r_a\xi \rho _a`$. For $`yy_a`$ (42) is solved to obtain $`Y_3(y)`$ only, but at $`y=y_a`$ the solution $`Y_4(y)`$ is introduced. We consider the Airy functions Ai($`q`$) and Bi($`q`$) Abram , where $`yy_a`$ and

$$q=[g^{}(y_a)]^{1/3}(yy_a).$$
(46)
The functions $`Y_3(y_a)`$ and $`Y_3^{}(y_a)`$ are approximately proportional to Ai(0) and Ai(0), and a second solution $`Y_4(y)`$ of (20), which is exponentially increasing as $`y\mathrm{}`$, is chosen to be proportional to Bi(0) at $`y=y_a`$. Hence

$$Y_4(y_a)=\sqrt{3}Y_3(y_a);Y_4^{}(y_a)=\sqrt{3}Y_3^{}(y_a),$$
(47)
and then the complex function

$$Y_5(y)=Y_3(y)\text{i}Y_4(y)$$
(48)
is propagated inwards over the range $`y_0yy_a`$.

We now consider the evaluation of the correction factor $`R(x)`$ where $`x=1/\rho `$. The solution $`F_a(\rho )`$ of (8) is written as

$$F(\rho )=R(x)F_a(\rho ),$$
(49)
where $`F(\rho )`$ is given by (19) in which

$$Y(y)Y_5(y),\rho _0\rho \rho _a;Y(y)Y_3(y),\rho >\rho _a.$$
(50)
We also define the functions $`\varphi (\rho )`$ and $`\varphi _a(\rho )`$ by

$$\varphi (\rho )=\frac{1}{F}\frac{dF}{d\rho };\varphi _a(\rho )=\frac{1}{F_a(\rho )}\frac{dF_a(\rho )}{d\rho }.$$
(51)
The function $`R(x)`$ then satisfies the equation

$$x^4\frac{d^2R}{dx^2}+2x^3\frac{dR}{dx}2\varphi _a(r)x^2\frac{dR}{dx}\frac{2V(x)}{\mathrm{}\omega }R(x)=0,$$
(52)
which is solved using the boundary conditions

$$R(0)=1;\frac{dR}{dx}|_0=0.$$
(53)
By making the choice for $`F_a(\rho )`$ given by (50) we ensure that $`\varphi _0(\rho )`$ is a slowly varying function for all $`\rho >\rho _0`$ and therefore that $`R(x)`$ is also a slowly varying function of $`x`$ over the range $`0xx_0`$, where $`x_0=1/\rho _0`$. Since $`r_a`$ may be of the order $`10^3a_0`$ to $`10^4a_0`$, the range $`r_0rr_a`$ can be quite large. The choice of the complex function here is crucial because it behaves like $`\mathrm{exp}(\text{i}\vartheta )`$ where $`\vartheta `$ varies approximately linearly with $`r`$ thus guaranteeing that $`\varphi _0(\rho )`$ changes slowly. Therefore the differential equation (52) for $`R(x)`$ can always be solved very accurately throughout the outer region $`0xx_0`$, using a grid method containing a maximum of 66 points and finally the solution $`F_{kl}(r)`$ required is obtained by taking the real part of the function in (49). In the inner region, $`0<rr_0`$, equation (7) is integrated numerically outwards using the Numerov algorithm. By using an iterative process to pinpoint the precise value of the energy, the values of $`\varphi (\rho )`$ in (51) obtained from the two regions are matched at $`\rho =\rho _0`$. The value of $`r_0`$ is varied within the range $`20a_0r_040a_0`$ and the results are shown to be extremely insensitive to its precise choice.

The second approach to the solution of the energy eigenvalue equation (5) is to use a scaled discrete variable representation (DVR) Tiesinga1998  (usually a Fourier grid) of the kinetic energy operator, and convert the bound state problem into one involving the diagonalization of a $`N\times N`$ matrix where $`N`$, for this single-channel problem, is equal to the number of grid points. Under a general real invertible transformation of the radial coordinate $`\rho `$ given by

$$t=u(\rho );\rho =u^1(t)U(t)$$
(54)
equation (8) becomes

$$\left[f^2\frac{\text{d}^2}{\text{d}t^2}f^2+P(t)\right]\stackrel{~}{F}(t)=\kappa \stackrel{~}{F}(t)$$
(55)
where

$$f(t)\left[\frac{\text{d}U}{\text{d}t}\right]^{\frac{1}{2}};\stackrel{~}{F}(t)\frac{F[U(t)]}{f(t)}$$
(56)
and

$$P(t)\frac{l(l+1)}{\rho ^2}+\rho ^2+\frac{2V(\rho )}{\mathrm{}\omega }+f^3(t)\frac{\text{d}^2f}{\text{d}t^2}.$$
(57)
Equation (55) is solved with the boundary conditions $`\stackrel{~}{F}(t_1)=0`$ and $`\stackrel{~}{F}(t_N)=0`$ by constructing a DVR using a set of basis functions $`\{\varphi _m(t)\}`$ and a finite set of coordinate points $`\{t_m\}`$ over the interval $`[t_1,t_N]`$. Then (55) is converted into $`N`$ linear equations specified by

$$\underset{j=1}{\overset{N}{}}\left[f^2(t_i)T_{ij}f^2(t_j)+P(t_j)\delta _{ij}\right]\stackrel{~}{F}(t_j)=\kappa \stackrel{~}{F}(t_i);i=1,2,\mathrm{}N.$$
(58)
For the Fourier basis defined by

$$\varphi _m(t)=\sqrt{\frac{2}{N}}\mathrm{sin}\left[\frac{m\pi (tt_1)}{t_Nt_1}\right];m=1,2,\mathrm{}N$$
(59)
the DVR of the kinetic energy operator is

$`T_{ij}={\displaystyle \frac{\pi ^2}{2(t_Nt_1)^2}}\{\begin{array}{cc}(1)^{ij}\left[\mathrm{csc}^2(\frac{\pi }{2N}(ij))\mathrm{csc}^2(\frac{\pi }{2N}(i+j))\right];\hfill & ij\hfill \\ & \\ \frac{1}{3}(2N^2+1)\mathrm{csc}^2\left(\frac{\pi i}{N}\right);\hfill & i=j\hfill \end{array}.`$ (63)
The scaling transformation used

$$t=u(\rho )=\left(\frac{\rho }{\zeta }\right)^{\frac{1}{p}};\rho =U(t)=\zeta t^p$$
(64)
involves the parameters $`\zeta `$ and $`\rho `$. Equation (64) gives $`t_1`$ and $`t_N`$ in terms of $`\rho _1`$ and $`\rho _N`$ where theoretically the boundaries are at $`\rho _1=0`$ and $`\rho _N=\mathrm{}`$. In practice, because of the extremely high potential barrier near the origin and the nature of the scaling, $`\rho _1`$ is chosen so that $`\rho _1\xi =2a_0`$. The choice of $`\rho _N`$ depends on the number of states being investigated as the associated eigenfunctions must have exponentially decayed sufficiently on the outer boundary. Taking $`\rho _N`$ = 15 allows us to calculate the positions of about the first 50 trap states accurately. The scaling parameters $`\zeta `$ and $`p`$ are chosen to ensure that a good proportion of the mesh points are inside the potential well around the point $`r=7a_0`$. On taking $`\zeta =20a_0`$ and $`p=10`$ about $`17\%`$ of the scaled mesh points lie between $`r_1`$ and $`\zeta `$, and with $`N=500`$, eight to ten digit convergence is obtained for the eigenvalues.

## V Application: Ultracold metastable helium atoms

As an application of the computational methods adopted in this paper, we consider the case of spin-polarized metastable helium tightly confined in harmonic traps of various frequencies. The colliding atoms are in the $`{}_{}{}^{5}\mathrm{\Sigma }_{g}^{+}`$ molecular state for which we use the potential of Stärck and Meyer Starck94 . This potential has a scattering length of $`156.777a_0`$ and supports 15 bound states. Calculated quantum defects for the 31 lowest trap states with $`l=0`$ are shown in Tables 1 and 2 for trapping frequencies ranging from 1 to 100 MHz. Also shown are the results obtained using the SC solution of (4). Very recently Gadéa et al. Gadea  have constructed a new $`{}_{}{}^{5}\mathrm{\Sigma }_{g}^{}`$ potential that supports the same number of states but has a scattering length of $`291a_0`$. Experiments Santos2001 , Robert2001  suggest that these two scattering lengths represent lower and upper limits on the exact value.

The QDT results are in excellent agreement with those obtained from the DVR method, with the absolute differences being $`O(10^6)`$. Also the agreement between the results from the QDT and SC approaches is very good, the absolute differences being only $`O(10^7)`$ at 1 MHz, increasing to $`O(10^5)`$ for 10, 20 and 50 MHz, and $`O(10^4)`$ at 100 MHz. For the trapping frequencies considered, the bound states in the $`{}_{}{}^{5}\mathrm{\Sigma }_{g}^{+}`$ potential are relatively unaffected by the presence of the harmonic trap. The number of bound states with $`E<0`$ is still $`n_b=15`$, but the most loosely bound state does show some unusual behavior, being shifted downwards at 1 MHz and upwards at 100 MHz. The trapping frequencies of 1, 5, 10, 50 and 100 MHz correspond to the ratios $`a/\xi =`$ 0.1167, 0.3691, 0.5221, 0.8255 and 1.167 respectively. The pseudopotential result (2), based upon the use of an energy-independent scattering length $`a`$, breaks down at the higher trapping frequencies. Allowance for the variation of $`a_{\text{eff}}(E)`$ with energy $`E`$ is essential in order to obtain the correct eigenvalues; $`a_{\text{eff}}(E)`$ changes sign as it passes through a divergence at $`k=0.01356a_0^1`$, producing a rapid variation in $`a_{\text{eff}}(E)/\xi `$ around $`E/\mathrm{}\omega 165.6/\nu (\text{MHz})`$.

As a test of the generalized self-consistent (GSC) method (37) we have considered the case of spin-polarized metastable helium atoms in a 10 MHz trap. We have chosen $`l=2`$ since for identical atoms, parity considerations exclude odd values of $`l`$, and in this case the $`{}_{}{}^{5}\mathrm{\Sigma }_{g}^{}`$ potential supports 14 bound states and $`b_2=2.383887\times 10^5a_0^4`$ in (41). Scaled energy eigenvalues $`n^{}`$ calculated for the 31 lowest trap states using the QDT and GSC methods are given in Table 3. The two sets of results agree to six significant figures. The associated quantum defects $`\mu ^{}`$ for $`l=2`$ are much smaller than the $`l=0`$ quantum defects, increasing from $`O(10^4)`$ for the lowest states to $`O(10^2)`$ for the highest states considered.

## VI Summary and discussion

A QDT analysis of colliding ultracold atoms tightly confined in harmonic potentials has been undertaken and a generalized SC model obtained. Two highly efficient computational methods are presented for calculating the bound state energies and eigenfunctions, one based upon QDT and the other on a scaled DVR method. The perturbed harmonic oscillator problem is characterized by a long asymptotic region beyond the effective range of the interatomic potential and the QDT method is very efficient for integrating inwards through this outer region. The radial Schrödinger equation for the relative motion of the harmonically confined atoms is transformed into that for a modified Coulomb potential and a quantum defect introduced that is an analytic function of energy. At large separations $`r`$, each eigenfunction for the confined atoms is expressed as a product of the unperturbed harmonic oscillator eigenfunction and a residual function $`R(x)`$ that is slowly varying in $`x=1/r`$. The unperturbed function is calculated by dividing the region into a number of ranges and expanding the function in a power series within each range. The differential equation for $`R(x)`$ is then solved accurately using a grid method requiring relatively few mesh points. The Schrödinger equation is integrated directly through the inner region and the two solutions matched at an intermediate separation in the range 20 to 40 $`a_0`$. The eigenvalues are determined very efficiently by iteration on an initial estimate obtained by extrapolating the quantum defects downwards in energy from highly excited trap states.

The two computational methods have been applied to the case of trapped spin-polarized metastable helium atoms. Energy eigenvalues calculated with the QDT method agree closely with those computed directly from the radial Schrödinger equation for the trapped atoms using a DVR method, and with those obtained using a self-consistent method (SC) involving an energy-dependent effective scattering length.

The range of trapping frequencies and the number of trap states considered in the present investigation was motivated by the desire to test the validity and robustness of our computational methods and not by current experimental feasibility. The higher trap frequencies are beyond those presently used and the higher trap states may not be physically realistic. The harmonic approximation to the confining potential of the optical lattice is only valid for the lower lying states of the atoms confined near the nodes or antinodes of the lattice. Also, the rate of quantum tunnelling from these higher states to neighboring wells in the lattice may be significant, thus requiring a calculation of the band structure type for the entire lattice.

The present calculations need to be extended to include collisional loss processes and to study the loss rate as a function of trapping frequency. In particular, for trapped metastable helium atoms, it will be important to include the magnetic dipole-dipole interactions that couple the spin-polarized $`{}_{}{}^{5}\mathrm{\Sigma }_{g}^{+}`$ state to the $`{}_{}{}^{1}\mathrm{\Sigma }_{g}^{+}`$ state from which there is a high probability of loss through Penning and associative ionization at small interatomic separations.