text:
# 1 Introduction

## 1 Introduction

In a previous paper , we studied the $`U(1)`$ and $`U(2)`$ instanton solutions on general noncommutative $`𝐑^\mathrm{𝟒}`$ , where for the $`U(1)`$ 1-instanton, 2-instanton and $`U(2)`$ 1-instanton cases we obtain explicitly the gauge field by ADHM construction .

In (noncommutative) ADHM construction, there are two steps to obtain the instanton configuration (the gauge field) explicitly. The first step is to solve the ADHM equation to get the matrix $`\mathrm{\Delta }`$ (we follow closely the notations of ). The second step is to construct the matrix $`U`$ from $`\mathrm{\Delta }`$, or in other words, to find all the zero modes of $`\mathrm{\Delta }`$. The gauge field can then be constructed from $`U`$. The difficulty of the first step has been well-known for a long time. But it was only recently realized that the second step is also a difficult problem in the noncommutative case.  discussed the $`U(1)`$ $`\theta _2>0`$ case, and our previous paper  resolved the problem in some of the $`U(1)`$ and $`U(2)`$ cases (see also ), but for the general cases there were no systematic methods to obtain the matrix $`U`$ yet. It would be a formidable task to find the matrix $`U`$ for higher $`N`$ and $`k`$.

Keeping this problem in mind, we will try to develop a method to accomplish the second step systematically. In fact there is a general method in the commutative case . What we will do is to extend this method to the noncommutative case. By using this method we have rechecked all the results obtained in our previous paper . As a new application of this method, we will demonstrate how to deal with some of the $`U(2)`$ 2-instanton cases. As far as we know there is no general solution of noncommutative $`U(2)`$ 2-instanton ADHM equation. So we will use this method to study a special case in the coincident instanton limit.

This paper is organized as follows: in section 2 we recall briefly the noncommutative $`𝐑^\mathrm{𝟒}`$ and set our notations and then the ADHM construction. In section 3 we present the modified matrix partition method and in section 4 we solve a special case of $`U(2)`$ 2-instanton. Finally we mention that our method can be used to give a more general class of $`U(2)`$ multi-instantons.

## 2 $`𝐑_{\mathrm{NC}}^\mathrm{𝟒}`$, the (anti-)self-dual equations and ADHM construction

### 2.1 $`𝐑_{\mathrm{NC}}^\mathrm{𝟒}`$ and the (anti-)self-dual equations

First let us recall briefly the noncommutative $`𝐑^\mathrm{𝟒}`$ and set our notations111For general reviews on noncommutative geometry and field theory, see, for example, .. For a general noncommutative $`𝐑^\mathrm{𝟒}`$ we mean a space with (operator) coordinates $`x^m`$, $`m=1,\mathrm{},4`$, which satisfy the following relations:

$$[x^m,x^n]=i\theta ^{mn},$$
(1)
where $`\theta ^{mn}`$ are real constants. If we assume the standard (Euclidean) metric for the noncommutative $`𝐑^\mathrm{𝟒}`$, we can use the orthogonal transformation with positive determinant to change $`\theta ^{mn}`$ into the following standard form:

$$(\theta ^{mn})=\left(\begin{array}{cccc}0& \theta ^{12}& 0& 0\\ \theta ^{12}& 0& 0& 0\\ 0& 0& 0& \theta ^{34}\\ 0& 0& \theta ^{34}& 0\end{array}\right),$$
(2)
with $`\theta ^{12}>0`$ and $`\theta ^{12}+\theta ^{34}0`$. If the original $`\theta ^{mn}`$ is non-degenerate, $`\theta ^{34}`$ in this standard form is non-vanishing. Otherwise the considered noncommutative $`𝐑^\mathrm{𝟒}`$ can be identified with a direct product $`𝐑_{\mathrm{NC}}^\mathrm{𝟐}\times 𝐑^\mathrm{𝟐}`$. One can refer  for this degenerate case, which we will not dwell on in this paper. By using this form of $`\theta ^{mn}`$, the only non-vanishing commutators are as follows:

$$[x^1,x^2]=i\theta ^{12},[x^3,x^4]=i\theta ^{34},$$
(3)
and the other twos obtained by using the anti-symmetric property of commutators. Introducing complex coordinates:

$$\begin{array}{cc}\hfill z_1=x^2+ix^1,& \overline{z}_1=x^2ix^1,\hfill \\ \hfill z_2=x^4+ix^3,& \overline{z}_2=x^4ix^3,\hfill \end{array}$$
(4)
the non-vanishing commutation relations are

$$[\overline{z}_1,z_1]=2\theta ^{12}\theta _1,[\overline{z}_2,z_2]=2\theta ^{34}\theta _2.$$
(5)
By a noncommutative gauge field $`A_m`$ we mean an operator valued field. The (anti-hermitian) field strength $`F_{mn}`$ is defined similarly as in the commutative case:

$$F_{mn}=\widehat{}_{[m}A_{n]}+A_{[m}A_{n]}\widehat{}_mA_n\widehat{}_nA_m+[A_m,A_n],$$
(6)
where the derivative operator $`\widehat{}_m`$ is defined as follows:

$$\widehat{}_mfi\theta _{mn}[x^n,f],$$
(7)
where $`\theta _{mn}`$ is the inverse of $`\theta ^{mn}`$. For our standard form (2) of $`\theta ^{mn}`$ we have

$$\widehat{}_1f=\frac{i}{\theta ^{12}}[x^2,f],\widehat{}_2f=\frac{i}{\theta ^{12}}[x^1,f],$$
(8)
which can be expressed by the complex coordinates (4) as follows:

$$_1f\widehat{}_{z_1}f=\frac{1}{\theta _1}[\overline{z}_1,f],\overline{}_1f\widehat{}_{\overline{z}_1}f=\frac{1}{\theta _1}[z_1,f],$$
(9)
and similar relations for $`x^{3,4}`$ and $`z_2,\overline{z}_2`$.

For a general metric $`g_{mn}`$ the instanton equations are

$$F_{mn}=\pm \frac{ϵ^{pqrs}}{2\sqrt{g}}g_{mp}g_{nq}F_{rs},$$
(10)
and the solutions are known as self-dual (SD, for “+” sign) and anti-self-dual (ASD, for “$``$” sign) instantons. Here $`ϵ^{pqrs}`$ is the totally anti-symmetric tensor ($`ϵ^{1234}=1`$ etc.) and $`g`$ is the metric. We will take the standard metric $`g_{mn}=\delta _{mn}`$ and take the noncommutative parameters $`\theta _{1,2}`$ as free parameters. We also note that the notions of self-dual and anti-self-dual are interchanged by a parity transformation. A parity transformation also changes the sign of $`\theta ^{mn}`$. In the following discussion we will consider only the ASD instantons. So we should not restrict $`\theta _2`$ to be positive.

### 2.2 ADHM construction for ordinary gauge theory

For ordinary gauge theory all the (ASD) instanton solutions are obtained by ADHM (Atiyah-Drinfeld-Hitchin-Manin) construction . In this construction we introduce the following ingredients (for $`U(N)`$ gauge theory with instanton number $`k`$):

* complex vector spaces $`V`$ and $`W`$ of dimensions $`k`$ and $`N`$,
* $`k\times k`$ matrix $`B_{1,2}`$, $`k\times N`$ matrix $`I`$ and $`N\times k`$ matrix $`J`$,
* the following quantities:

$`\mu _r`$ $`=`$ $`[B_1,B_1^{}]+[B_2,B_2^{}]+II^{}J^{}J,`$ (11)
$`\mu _c`$ $`=`$ $`[B_1,B_2]+IJ.`$ (12)

The claim of ADHM is as follows:

* Given $`B_{1,2}`$, $`I`$ and $`J`$ such that $`\mu _r=\mu _c=0`$, an ASD gauge field can be constructed;
* All ASD gauge fields can be obtained in this way.

It is convenient to introduce a quaternionic notation for the 4-dimensional Euclidean space-time indices:

$$xx^n\sigma _n,\overline{x}x^n\overline{\sigma }_n,$$
(13)
where $`\sigma _n=(i\stackrel{}{\tau },1)`$ and $`\tau ^c`$, $`c=1,2,3`$ are the three Pauli matrices, and the conjugate matrices $`\overline{\sigma }_n=\sigma _n^{}=(i\stackrel{}{\tau },1)`$. In terms of the complex coordinates (4) we have

$$(x_{\alpha \dot{\alpha }})=\left(\begin{array}{cc}z_2& z_1\\ \overline{z}_1& \overline{z}_2\end{array}\right),(\overline{x}^{\dot{\alpha }\alpha })=\left(\begin{array}{cc}\overline{z}_2& z_1\\ \overline{z}_1& z_2\end{array}\right).$$
(14)
Then the basic object in the ADHM construction is the $`(N+2k)\times 2k`$ matrix $`\mathrm{\Delta }`$ which is linear in the space-time coordinates:

$$\mathrm{\Delta }=a+b\overline{x},$$
(15)
where the constant matrices

$$a=\left(\begin{array}{cc}I^{}& J\\ B_2^{}& B_1\\ B_1^{}& B_2\end{array}\right),b=\left(\begin{array}{cc}0& 0\\ 1& 0\\ 0& 1\end{array}\right).$$
(16)
Consider the conjugate operator of $`\mathrm{\Delta }`$:

$$\mathrm{\Delta }^{}=a^{}+xb^{}=\left(\begin{array}{ccc}I& B_2+z_2& B_1+z_1\\ J^{}& B_1^{}\overline{z}_1& B_2^{}+\overline{z}_2\end{array}\right).$$
(17)
It is easy to check that the ADHM equations (11) and (12) are equivalent to the so-called factorization condition:

$$\mathrm{\Delta }^{}\mathrm{\Delta }=\left(\begin{array}{cc}f^1& 0\\ 0& f^1\end{array}\right),$$
(18)
where $`f(x)`$ is a $`k\times k`$ hermitian matrix. From the above condition we can construct a hermitian projection operator $`P`$ as follows:222We use the following abbreviation for expressions with $`f`$:
$$\mathrm{\Delta }f\mathrm{\Delta }^{}\mathrm{\Delta }\left(\begin{array}{cc}f& 0\\ 0& f\end{array}\right)\mathrm{\Delta }^{}=\mathrm{\Delta }(f1_2)\mathrm{\Delta }^{}.$$
(19)

$`P`$ $`=`$ $`\mathrm{\Delta }f\mathrm{\Delta }^{},`$ (20)
$`P^2`$ $`=`$ $`\mathrm{\Delta }ff^1f\mathrm{\Delta }^{}=P.`$ (21)
Obviously, the null-space of $`\mathrm{\Delta }^{}(x)`$ is of $`N`$ dimension for generic $`x`$. The basis vectors for this null-space can be assembled into an $`(N+2k)\times N`$ matrix $`U(x)`$:

$$\mathrm{\Delta }^{}U=0,$$
(22)
which can be chosen to satisfy the following orthonormal condition:

$$U^{}U=1.$$
(23)
The above orthonormal condition guarantees that $`UU^{}`$ is also a hermitian projection operator. Now it can be proved that the completeness relation

$$P+UU^{}=1$$
(24)
holds if $`U`$ contains the whole null-space of $`\mathrm{\Delta }^{}`$. In other words, this completeness relation requires that $`U`$ consists of all the zero modes of $`\mathrm{\Delta }^{}`$.333The proof is sketched as follows: the two projection operators $`P`$ and $`UU^{}`$ are orthogonal to each other, and so $`1PUU^{}`$ is also a hermitian projection operator. Now this can always be written as the form $`VV^{}`$. Here $`\mathrm{\Delta }`$ and $`f`$ are both of maximum rank and $`PVV^{}=0`$, then $`V`$ must consist of some zero modes of $`\mathrm{\Delta }^{}`$ other than those in $`U`$. This conclusion is in conflict with the assumption that $`U`$ contains all the zero modes of $`\mathrm{\Delta }^{}`$. Although the notion of maximum rank is ambiguous in the infinite-dimensional case as we will encounter in noncommutative ADHM construction, it can be defined as follows so that the above proof is also true in the infinite-dimensional case. We say that an $`\mathrm{}\times \mathrm{}`$ matrix $`X`$ is of maximum rank if it has no zero modes, i.e., $`X\mathrm{\Phi }0`$ if $`\mathrm{\Phi }0`$.

The (anti-hermitian) gauge potential is constructed from $`U`$ by the following formula:

$$A_m=U^{}_mU.$$
(25)
Substituting this expression into (6), we get the following field strength:

$`F_{mn}`$ $`=`$ $`_{[m}(U^{}_{n]}U)+(U^{}_{[m}U)(U^{}_{n]}U)=_{[m}U^{}(1UU^{})_{n]}U`$ (26)
$`=`$ $`_{[m}U^{}\mathrm{\Delta }f\mathrm{\Delta }^{}_{n]}U=U^{}_{[m}\mathrm{\Delta }f_{n]}\mathrm{\Delta }^{}U=U^{}b\overline{\sigma }_{[m}\sigma _{n]}fb^{}U`$
$`=`$ $`2i\overline{\eta }_{mn}^cU^{}b(\tau ^cf)b^{}U.`$
Here $`\overline{\eta }_{ij}^a`$ is the standard ’t Hooft $`\eta `$-symbol, which is anti-self-dual:

$$\frac{1}{2}ϵ_{ijkl}\overline{\eta }_{kl}^a=\overline{\eta }_{ij}^a.$$
(27)
### 2.3 Noncommutative ADHM construction

The above construction has been extended to noncommutative gauge theory . We recall this construction briefly here. By introducing the same data as above but considering the $`z_i`$’s as noncommutative we see that the factorization condition (18) still gives $`\mu _c=0`$, but $`\mu _r`$ no longer vanishes. It is easy to check that the following relation holds:

$$\mu _r=\zeta \theta _1+\theta _2.$$
(28)
In this case the two ADHM equations (11) and (12) can be combined into one :

$$\tau ^{c\dot{\alpha }}{}_{\dot{\beta }}{}^{}(\overline{a}^{\dot{\beta }}a_{\dot{\alpha }})_{ij}^{}=\delta _{ij}\delta ^{c3}\zeta .$$
(29)
As studied mathematically by various people (see, for example, the lectures by H. Nakajima ), the moduli space of the noncommutative instantons is better behaved than their commutative counterpart. In the noncommutative case the operator $`\mathrm{\Delta }^{}\mathrm{\Delta }`$ always has maximum rank (see ).

Though there is no much difference between the noncommutative ADHM construction and the commutative one, we should study the noncommutative case in more detail. In order to study the instanton solutions precisely, we use a Fock space representation as follows ($`n_1,n_20`$):

$`z_1|n_1,n_2`$ $`=`$ $`\sqrt{\theta _1}\sqrt{n_1+1}|n_1+1,n_2,`$ (30)
$`\overline{z}_1|n_1,n_2`$ $`=`$ $`\sqrt{\theta _1}\sqrt{n_1}|n_11,n_2,`$ (31)
by using the commutation relation (5). Similar expressions for $`z_2`$ and $`\overline{z}_2`$ also apply (but paying a little attention to the sign of $`\theta _2`$ which is not restricted to be positive). In this representation the $`z_i`$’s are infinite-dimensional matrices, and so are the operator $`\mathrm{\Delta }`$, $`\mathrm{\Delta }^{}`$ etc. Because of infinite dimensions are involved we can not determine the dimension of null-space of $`\mathrm{\Delta }^{}`$ straightforwardly from the difference of the numbers of its rows and columns. But it turns out that $`\mathrm{\Delta }^{}`$ also has infinite number of zero modes, and they can be arranged into an $`(N+2k)\times N`$ matrix with entries from the (noncommutative) algebra generated by the coordinates, which resembles the commutative case.

## 3 The modified matrix partition method

In the commutative case, there is a standard method , which we call the ‘matrix partition method’ in ADHM construction, to obtain the matrix $`U`$ from $`\mathrm{\Delta }`$ as follows. First we introduce the following decomposition of the matrices $`U`$ and $`\mathrm{\Delta }`$:

$$U=\left(\begin{array}{c}V_{N\times N}\\ U_{2k\times N}^{}\end{array}\right),\mathrm{\Delta }=\left(\begin{array}{c}K_{N\times 2k}\\ \mathrm{\Delta }_{2k\times 2k}^{}\end{array}\right).$$
(32)
So we have

$$K=\left(\begin{array}{cc}I^{}& J\end{array}\right),\mathrm{\Delta }^{}=\left(\begin{array}{cc}B_2^{}+\overline{z}_2& B_1z_1\\ B_1^{}+\overline{z}_1& B_2+z_2\end{array}\right).$$
(33)
Substituting (32) into the completeness relation (24), we obtain

$$\left(\begin{array}{cc}VV^{}& VU^{}\\ U^{}V^{}& U^{}U^{}\end{array}\right)=\left(\begin{array}{cc}1KfK^{}& Kf\mathrm{\Delta }^{}\\ \mathrm{\Delta }^{}fK^{}& 1\mathrm{\Delta }^{}f\mathrm{\Delta }^{}\end{array}\right).$$
(34)
Now we can choose

$$V=V^{}=(1KfK^{})^{1/2}$$
(35)
as a solution of the (1,1) element of the matrix equation (34). Then by the (2,1) element of (34) the matrix $`U^{}`$ can be expressed as

$$U^{}=\mathrm{\Delta }^{}fK^{}V^1.$$
(36)
The choice (35) of $`V`$ is known as the ‘singular gauge’, and any other choices of $`V`$ which solve (34) are related to (35) by a gauge transformation.

Corresponding to the singularity arising in the commutative case, the choice (35) also brings problems if we extend the matrix partition method to the noncommutative case. In the following subsections we will show that slight modifications of this method can remedy these problems.

In the noncommutative case, the solutions of the (1,1) element of (34) are not always related to each other by gauge transformations. This can be seen as follows. If $`V`$ is a solution of the (1,1) element of (34) and $`S`$ is an operator such that $`SS^{}=1`$, then $`VS`$ is also a solution:

$$(VS)(VS)^{}=VSS^{}V^{}=VV^{},$$
(37)
while we can have $`S^{}S1`$. This $`S`$ induces a new solution $`US`$ of the completeness relation (24). However, $`US`$ does not satisfy the orthonormal condition (23):

$$(US)^{}(US)=S^{}U^{}US1,$$
(38)
which is not a problem in the commutative case: any solution of (24) must satisfy (23). The above freedom of $`V`$ can be used to resolve the problems of the matrix partition method in noncommutative ADHM construction, while we must consider the completeness relation (24) and the orthonormal condition (23) as two independent equations.

The two cases, $`\theta _2>0`$ and $`\theta _2<0`$, have an essential difference: for $`\theta _2>0`$ the matrix $`\mathrm{\Delta }^{}`$ has zero modes and $`\mathrm{\Delta }^{}`$ is of maximum rank, while for $`\theta _2<0`$ the matrix $`\mathrm{\Delta }^{}`$ is of maximum rank and $`\mathrm{\Delta }^{}`$ has zero modes. As we can see in the following subsections, the two cases must be considered separately just because of this difference.

### 3.1 $`\theta _2>0`$ case

The immediate consequence of that $`\mathrm{\Delta }^{}`$ has zero modes in this case is that the matrix $`V`$ in (35) is not invertible and we we can not use (36) to determine $`U^{}`$. Here we will try to directly project out the zero modes of $`VV^{}`$ first.

Let $`\mathrm{\Pi }_0`$ be the projector whose range is the null-space of $`(1KfK^{})`$. We can still use (35) and let

$$U^{}=\mathrm{\Delta }^{}fK^{}(1KfK^{})_\mathrm{\Pi }^{1/2},$$
(39)
where the inversion of $`(1KfK^{})`$ is restricted on $`\mathrm{\Pi }1\mathrm{\Pi }_0`$. Noting that

$$\begin{array}{cc}\hfill \mathrm{\Delta }^{}\mathrm{\Delta }^{}fK^{}\mathrm{\Pi }_0& =\mathrm{\Delta }^{}\mathrm{\Delta }^{}(K^{}K+\mathrm{\Delta }^{}\mathrm{\Delta }^{})^1K^{}\mathrm{\Pi }_0\hfill \\ & =[1K^{}K(K^{}K+\mathrm{\Delta }^{}\mathrm{\Delta }^{})^1]K^{}\mathrm{\Pi }_0\hfill \\ & =K^{}(1KfK^{})\mathrm{\Pi }_0=0\hfill \end{array}$$
(40)
and $`\mathrm{\Delta }^{}`$ is of maximum rank by (33), we have

$$\mathrm{\Delta }^{}fK^{}\mathrm{\Pi }_0=0.$$
(41)
The (2,1) element of the matrix equation (34) is now satisfied:

$$\begin{array}{cc}\hfill U^{}V^{}& =\mathrm{\Delta }^{}fK^{}(1KfK^{})_\mathrm{\Pi }^{1/2}(1KfK^{})^{1/2}\hfill \\ & =\mathrm{\Delta }^{}fK^{}(1\mathrm{\Pi }_0)=\mathrm{\Delta }^{}fK^{}.\hfill \end{array}$$
(42)
Moreover we can easily check that the (2,2) element of (34) holds:

$$\begin{array}{cc}\hfill \mathrm{\Delta }^{}U^{}U^{}& =\mathrm{\Delta }^{}\mathrm{\Delta }^{}fK^{}(1KfK^{})_\mathrm{\Pi }^1Kf\mathrm{\Delta }^{}\hfill \\ & =K^{}(1KfK^{})(1KfK^{})_\mathrm{\Pi }^1Kf\mathrm{\Delta }^{}\hfill \\ & =K^{}(1\mathrm{\Pi }_0)Kf\mathrm{\Delta }^{}=K^{}Kf\mathrm{\Delta }^{}\hfill \end{array}$$
(43)
and

$$\begin{array}{cc}\hfill \mathrm{\Delta }^{}(1\mathrm{\Delta }^{}f\mathrm{\Delta }^{})& =[1\mathrm{\Delta }^{}\mathrm{\Delta }^{}(K^{}K+\mathrm{\Delta }^{}\mathrm{\Delta }^{})^1]\mathrm{\Delta }^{}\hfill \\ & =K^{}Kf\mathrm{\Delta }^{},\hfill \end{array}$$
(44)
so

$$U^{}U^{}=1\mathrm{\Delta }^{}f\mathrm{\Delta }^{}.$$
(45)
Even though (35) together with (39) satisfy equation (34), they do not satisfy the orthonormal condition (23):

$$\begin{array}{cc}\hfill U^{}U& =V^{}V+U^{}U^{}\hfill \\ & =1KfK^{}+(1KfK^{})_\mathrm{\Pi }^{1/2}Kf\mathrm{\Delta }^{}\mathrm{\Delta }^{}fK^{}(1KfK^{})_\mathrm{\Pi }^{1/2}\hfill \\ & =1KfK^{}\hfill \\ & +(1KfK^{})_\mathrm{\Pi }^{1/2}Kf(f^1K^{}K)fK^{}(1KfK^{})_\mathrm{\Pi }^{1/2}\hfill \\ & =1KfK^{}\hfill \\ & +(1KfK^{})_\mathrm{\Pi }^{1/2}(KfK^{}KfK^{}KfK^{})(1KfK^{})_\mathrm{\Pi }^{1/2}\hfill \\ & =1KfK^{}+\mathrm{\Pi }KfK^{}\mathrm{\Pi }\hfill \\ & =1\mathrm{\Pi }_0KfK^{}KfK^{}\mathrm{\Pi }_0+\mathrm{\Pi }_0=\mathrm{\Pi }1.\hfill \end{array}$$
(46)
The solution to this problem is to introduce a shift operator $`s`$ as follows:

$$ss^{}=1,s^{}s=\mathrm{\Pi },$$
(47)
and we set

$`V`$ $`=`$ $`(1KfK^{})^{1/2}s^{},`$ (48)
$`U^{}`$ $`=`$ $`\mathrm{\Delta }^{}fK^{}(1KfK^{})^{1/2}s^{}.`$ (49)
The operator $`s^{}`$ removes the zero modes of the operator $`(1KfK^{})`$ and makes its inversion well-defined. Then we can easily see that the completeness relation (34) still holds and the orthonormal condition is now satisfied:

$$\begin{array}{cc}\hfill U^{}U& =s(1KfK^{})s^{}\hfill \\ & +s(1KfK^{})^{1/2}Kf\mathrm{\Delta }^{}\mathrm{\Delta }^{}fK^{}(1KfK^{})^{1/2}s^{}\hfill \\ & =s\mathrm{\Pi }s^{}=1.\hfill \end{array}$$
(50)
So (48) and (49) constitute the required matrix $`U`$.

Now we determine the null-space of $`(1KfK^{})`$. As

$$\begin{array}{cc}\hfill (1KfK^{})K& =K[1(K^{}K+\mathrm{\Delta }^{}\mathrm{\Delta }^{})^1K^{}K]\hfill \\ & =Kf\mathrm{\Delta }^{}\mathrm{\Delta }^{},\hfill \end{array}$$
(51)
we can first determine the null-space of $`\mathrm{\Delta }^{}`$. It is easy to see from (33) that the zero modes of $`\mathrm{\Delta }^{}`$ must take the form:

$$\mathrm{\Phi }=\left(\begin{array}{c}\mathrm{{\rm Y}}\\ 0\end{array}\right)$$
(52)
because the second column of $`\mathrm{\Delta }^{}`$ is of maximum rank. By (51), we introduce an operator

$$\mathrm{\Psi }K\mathrm{\Phi }=I^{}\mathrm{{\rm Y}},$$
(53)
which belongs to the null-space of $`(1KfK^{})`$:

$$(1KfK^{})\mathrm{\Psi }=Kf\mathrm{\Delta }^{}\mathrm{\Delta }^{}\mathrm{\Phi }=0.$$
(54)
Then a hermitian operator $`G`$ can be constructed as follows:

$$G=\mathrm{\Psi }^{}\mathrm{\Psi }=\mathrm{{\rm Y}}^{}II^{}\mathrm{{\rm Y}}=\mathrm{{\rm Y}}^{}[II^{}+\underset{\alpha =1,2}{}(B_\alpha +z_\alpha )(B_\alpha ^{}+\overline{z}_\alpha )]\mathrm{{\rm Y}}=\mathrm{{\rm Y}}^{}f^1\mathrm{{\rm Y}},$$
(55)
which is positive definite because $`f^1`$ has no zero mode. Here we have used the fact

$$(B_\alpha ^{}+\overline{z}_\alpha )\mathrm{{\rm Y}}=0,\alpha =1,2.$$
(56)
If $`\mathrm{\Psi }`$ contains all the zero modes of $`(1KfK^{})`$, it is easy to show that

$$\mathrm{\Pi }_0=\mathrm{\Psi }G^1\mathrm{\Psi }^{}.$$
(57)
We will prove this completeness of $`\mathrm{\Psi }`$ as follows. From (41) we have

$`fI\mathrm{\Pi }_0`$ $``$ $`\mathrm{{\rm Y}},`$ (58)
$`J^{}\mathrm{\Pi }_0`$ $`=`$ $`0,`$ (59)
where the notation ‘$``$’ now means ‘included in the space spanned by the column vectors of’, and (59) gives

$$(1KfK^{})\mathrm{\Pi }_0=(1I^{}fI)\mathrm{\Pi }_0=0,$$
(60)
which leads to

$$\mathrm{\Pi }_0=I^{}fI\mathrm{\Pi }_0I^{}\mathrm{{\rm Y}}.$$
(61)
So (51) does span the whole null-space of $`(1KfK^{})`$.

For $`[B_1,B_2]0`$, we can not give any more explicit results. But when $`[B_1,B_2]=0`$, we have :

$$\mathrm{{\rm Y}}=e^{_\alpha \theta _\alpha ^1B_\alpha ^{}z_\alpha }|0,0.$$
(62)
In general, the condition $`[B_1,B_2]=0`$ can not be satisfied. When $`N=1`$ it can be proved  that we can always take $`J=0`$, so the above condition holds. Moreover, it also holds in the $`U(2)`$ 1-instanton case, which we had considered in . A special $`U(2)`$ 2-instanton case also satisfies this condition. We will give the explicit solution in the next section.

### 3.2 $`\theta _2<0`$ case

In this case the problem is that $`\mathrm{\Delta }^{}`$ has zero modes. Here $`\mathrm{\Delta }^{}`$ is of maximum rank as we can see from (33), so is $`(1KfK^{})`$ for the reason stated in the previous subsection. But if we use (35) and (36) as the solution of (34), equations (43) (where $`\mathrm{\Pi }_0=0`$) and (44) do not give (45).

Suppose that the null-space of $`\mathrm{\Delta }^{}`$ is of $`k^{}`$ dimension and $`\mathrm{\Gamma }`$ is the matrix composed of all the orthonormal zero modes of $`\mathrm{\Delta }^{}`$:

$$\mathrm{\Gamma }^{}\mathrm{\Gamma }=1_{k^{}\times k^{}},\mathrm{\Delta }^{}\mathrm{\Gamma }=0,$$
(63)
we have

$$1\mathrm{\Delta }^{}f\mathrm{\Delta }^{}U^{}U^{}\mathrm{\Gamma }.$$
(64)
Noting that a matrix $`\stackrel{~}{U}`$ composed of (35) and (36) satisfies

$$\stackrel{~}{U}^{}\stackrel{~}{U}=1,\mathrm{\Delta }^{}\stackrel{~}{U}=0,$$
(65)
the only thing we must do is to add some extra zero modes to $`\stackrel{~}{U}`$ so as to make it satisfy the completeness relation. In what follows, we will prove that these extra zero modes can be assembled into a matrix $`X`$ of the form

$$X=\left(\begin{array}{c}0\\ \mathrm{\Gamma }\end{array}\right).$$
(66)
Firstly, $`X`$ belongs to the null-space of $`\mathrm{\Delta }^{}`$:

$$\mathrm{\Delta }^{}X=\mathrm{\Delta }^{}\mathrm{\Gamma }=0.$$
(67)
Secondly, it is not difficult to see that a non-vanishing (1,1) element of $`X`$ will break down the (1,1) element of equation (34). Finally, any vectors belong to the space spanned by $`\mathrm{\Gamma }`$ appearing in the (2,1) element of $`X`$ will violate (64). So we conclude that the combination of $`\stackrel{~}{U}`$ and $`X`$ gives all the zero modes of $`\mathrm{\Delta }^{}`$.

To make the above conclusion more explicit, we introduce the following shift operator:

$$SS^{}=1,S^{}S=1\underset{n=0}{\overset{k^{}1}{}}|n,0n,0|$$
(68)
and a $`k^{}\times N`$ matrix

$$Y=\left(\begin{array}{cccc}0,0|& 0& \mathrm{}& 0\\ 1,0|& 0& \mathrm{}& 0\\ \mathrm{}& \mathrm{}& \mathrm{}& \mathrm{}\\ k^{}1,0|& 0& \mathrm{}& 0\end{array}\right).$$
(69)
An $`N\times N`$ matrix operator $`𝒮`$ is defined as

$$𝒮=\mathrm{diag}(S,1,\mathrm{},1).$$
(70)
Then the required matrix $`U`$ can be expressed as follows:

$`V`$ $`=`$ $`(1KfK^{})^{1/2}𝒮,`$ (71)
$`U^{}`$ $`=`$ $`\mathrm{\Gamma }Y\mathrm{\Delta }^{}fK^{}(1KfK^{})^{1/2}𝒮,`$ (72)
which can be straightforwardly checked to satisfy both conditions (23) and (24).

Again for $`[B_1,B_2]=0`$ we can explicitly give

$$\mathrm{\Gamma }=\left(\begin{array}{c}e^{\theta _2^1B_2\overline{z}_2\theta _1^1B_1^{}z_1}|0,0\\ 0\end{array}\right)W^{1/2},$$
(73)
where the normalization factor $`W`$ is:

$$W=0,0|e^{\theta _2^1B_2^{}z_2\theta _1^1B_1\overline{z}_1}e^{\theta _2^1B_2\overline{z}_2\theta _1^1B_1^{}z_1}|0,0,$$
(74)
and so we have $`k^{}=k`$.

## 4 A special $`U(2)`$ 2-instanton solution

To explicitly give the moduli space of a $`U(N)`$ multi-instanton is a very difficult problem . In the commutative case,  solved this problem when the instanton number $`k=2`$. For the noncommutative case a similar solution has not appeared yet. However, it is much easier to find some special solutions of $`U(N)`$ multi-instanton. To illustrate our method, we construct one of noncommutative $`U(2)`$ 2-instanton in the limit of coincident instantons.

In a special case, these ADHM data (15) can be written as follows:

$$\mathrm{\Delta }=\left(\begin{array}{cccc}0& \sqrt{2\zeta +\rho ^2}& 0& 0\\ 0& 0& \rho & 0\\ \overline{z}_2& 0& z_1& \sqrt{\zeta +\rho ^2}\\ 0& \overline{z}_2& 0& z_1\\ \overline{z}_1& 0& z_2& 0\\ \sqrt{\zeta +\rho ^2}& \overline{z}_1& 0& z_2\end{array}\right),$$
(75)
$$\mathrm{\Delta }^{}=\left(\begin{array}{cccccc}0& 0& z_2& 0& z_1& \sqrt{\zeta +\rho ^2}\\ \sqrt{2\zeta +\rho ^2}& 0& 0& z_2& 0& z_1\\ 0& \rho & \overline{z}_1& 0& \overline{z}_2& 0\\ 0& 0& \sqrt{\zeta +\rho ^2}& \overline{z}_1& 0& \overline{z}_2\end{array}\right),$$
(76)
where $`\rho `$ is the scale of the instanton size. In other words, we have

$$I=\left(\begin{array}{cc}0& 0\\ \sqrt{2\zeta +\rho ^2}& 0\end{array}\right),J=\left(\begin{array}{cc}0& 0\\ \rho & 0\end{array}\right),$$
(77)
$$B_1=\left(\begin{array}{cc}0& \sqrt{\zeta +\rho ^2}\\ 0& 0\end{array}\right),B_2=0.$$
(78)
These ADHM data lead to

$$f=\left(\begin{array}{cc}(Z_1^{\prime \prime }+Z_2^{}+\zeta +\rho ^2)(Z_\rho ^{\prime \prime }{}_{}{}^{})^1& \sqrt{\zeta +\rho ^2}(Z_\rho ^{\prime \prime }{}_{}{}^{})^1\overline{z}_1\\ \sqrt{\zeta +\rho ^2}z_1(Z_\rho ^{\prime \prime }{}_{}{}^{})^1& (Z_1^{}+Z_2^{}\theta _1+\rho ^2)(Z_\rho ^{}{}_{}{}^{})^1\end{array}\right).$$
(79)
In order to make the above expressions more similar to that of the $`U(1)`$ 2-instanton discussed in , we have carefully chosen our notations as follows:

$$Z_1z_1\overline{z}_1,Z_2z_2\overline{z}_2,$$
(80)
$`Z_1^{}Z_1+\theta _1,Z_1^{\prime \prime }Z_1+2\theta _1,\mathrm{etc}.`$ (81)
$`Z_2^{}Z_2+\theta _2,Z_2^{\prime \prime }Z_2+2\theta _2,\mathrm{etc}.`$ (82)
$`Z`$ $``$ $`(Z_1+Z_2)^2\theta _1Z_1+(\theta _2+\rho ^2)Z_2,`$ (83)
$`Z^{}^{}`$ $``$ $`(Z_1^{}+Z_2^{})^2\theta _1Z_1^{}+(\theta _2+\rho ^2)Z_2^{},\mathrm{etc}.`$ (84)
$`Z_\rho `$ $``$ $`Z+\rho ^2(Z_1+Z_2+\zeta +\rho ^2),`$ (85)
$`Z_\rho ^{}^{}`$ $``$ $`Z^{}{}_{}{}^{}+\rho ^2(Z_1^{}+Z_2^{}+\zeta +\rho ^2),\mathrm{etc}.`$ (86)
It is not difficult to obtain

$$1KfK^{}=1\left(\begin{array}{cc}(2\zeta +\rho ^2)f_{22}& 0\\ 0& \rho ^2f_{11}\end{array}\right)=\left(\begin{array}{cc}Z/Z_\rho ^{}^{}& 0\\ 0& Z^{\prime \prime }{}_{}{}^{}/Z_\rho ^{\prime \prime }^{}\end{array}\right),$$
(87)
$`(\mathrm{\Delta }^{}fK^{})_1`$ $`=`$ $`\sqrt{2\zeta +\rho ^2}\left(\begin{array}{c}\sqrt{\zeta +\rho ^2}\overline{z}_2\overline{z}_1\\ \overline{z}_2(Z_1+Z_2+\theta _2+\rho ^2)\\ \sqrt{\zeta +\rho ^2}\overline{z}_1\overline{z}_1\\ \overline{z}_1(Z_1+Z_2\theta _1)\end{array}\right)/Z_\rho ^{}{}_{}{}^{},`$ (92)
$`(\mathrm{\Delta }^{}fK^{})_2`$ $`=`$ $`\rho \left(\begin{array}{c}z_1(Z_1^{\prime \prime }+Z_2^{})\\ \sqrt{\zeta +\rho ^2}z_1z_1\\ z_2(Z_1^{\prime \prime }+Z_2^{}+\zeta +\rho ^2)\\ \sqrt{\zeta +\rho ^2}z_2z_1\end{array}\right)/Z_\rho ^{\prime \prime }{}_{}{}^{},`$ (97)
where the subscripts ‘\*1’ and ‘\*2’ mean the first column and the second column of the matrix $`\mathrm{\Delta }^{}fK^{}`$ respectively. Now we solve the two cases $`\theta _2>0`$ and $`\theta _2<0`$ in sequence.

### 4.1 The $`\theta _2>0`$ case

Because $`[B_1,B_2]=0`$ for this special solution, we have

$$\mathrm{\Psi }=I^{}e^{_\alpha \theta _\alpha ^1B_\alpha ^{}z_\alpha }|0,0=\sqrt{2\zeta +\rho ^2}\left(\begin{array}{cc}\sqrt{\frac{\zeta +\rho ^2}{\theta _1}}|1,0& |0,0\\ 0& 0\end{array}\right),$$
(98)
so

$$\mathrm{\Pi }_0=\mathrm{\Psi }(\mathrm{\Psi }^{}\mathrm{\Psi })^1\mathrm{\Psi }^{}=\left(\begin{array}{cc}|0,00,0|+|1,01,0|& 0\\ 0& 0\end{array}\right).$$
(99)
In fact, a simple comparison between (56) in  and the above (83) makes it clear that here $`Z`$ also annihilates two states: $`|0,0`$ and $`|1,0`$, so (99) can be obtained directly.

By (47) we introduce a shift operator

$$s=\left(\begin{array}{cc}t& 0\\ 0& 1\end{array}\right),$$
(100)
$$tt^{}=1,t^{}t=1|0,00,0||1,01,0|.$$
(101)
The required $`U`$ is given as follows:

$$V=\left(\begin{array}{cc}(Z/Z_\rho ^{}{}_{}{}^{})^{1/2}t^{}& 0\\ 0& (Z^{\prime \prime }{}_{}{}^{}/Z_\rho ^{\prime \prime }{}_{}{}^{})^{1/2}\end{array}\right),$$
(102)
$`U_1^{}`$ $`=`$ $`\sqrt{2\zeta +\rho ^2}\left(\begin{array}{c}\sqrt{\zeta +\rho ^2}\overline{z}_2\overline{z}_1\\ \overline{z}_2(Z_1+Z_2+\theta _2+\rho ^2)\\ \sqrt{\zeta +\rho ^2}\overline{z}_1\overline{z}_1\\ \overline{z}_1(Z_1+Z_2\theta _1)\end{array}\right)(ZZ_\rho ^{}{}_{}{}^{})^{1/2}t^{},`$ (107)
$`U_2^{}`$ $`=`$ $`\rho \left(\begin{array}{c}z_1(Z_1^{\prime \prime }+Z_2^{})\\ \sqrt{\zeta +\rho ^2}z_1z_1\\ z_2(Z_1^{\prime \prime }+Z_2^{}+\zeta +\rho ^2)\\ \sqrt{\zeta +\rho ^2}z_2z_1\end{array}\right)(Z^{\prime \prime }{}_{}{}^{}Z_{\rho }^{\prime \prime }{}_{}{}^{})^{1/2}.`$ (112)
### 4.2 The $`\theta _2<0`$ case

Unlike the $`U(1)`$ 2-instanton case in , the (2,2) element of the matrix $`f`$ is well-defined on the vacuum $`|0,0`$:

$$f_{22}|0,0=(2\theta _1+\theta _2+\rho ^2)^1|0,0.$$
(113)
Note that this $`f_{22}`$ tends to that of the $`U(1)`$ 2-instanton case when $`\rho `$ tends to zero, as we have expected.

Now we have

$$e^{\theta _2^1B_2\overline{z}_2\theta _1^1B_1^{}z_1}|0,0=\left(\begin{array}{cc}|0,0& 0\\ \sqrt{\frac{\zeta +\rho ^2}{\theta _1}}|1,0& |0,0\end{array}\right),$$
(114)
and so

$$\mathrm{\Gamma }=\left(\begin{array}{cc}\sqrt{\frac{\theta _1}{2\theta _1+\theta _2+\rho ^2}}|0,0& 0\\ \sqrt{\frac{\zeta +\rho ^2}{2\theta _1+\theta _2+\rho ^2}}|1,0& |0,0\\ 0& 0\\ 0& 0\end{array}\right),$$
(115)
Then the required $`U`$ is given as follows:

$$V=\left(\begin{array}{cc}(Z/Z_\rho ^{}{}_{}{}^{})^{1/2}t& 0\\ 0& (Z^{\prime \prime }{}_{}{}^{}/Z_\rho ^{\prime \prime }{}_{}{}^{})^{1/2}\end{array}\right),$$
(116)
$$\begin{array}{cc}\hfill U_1^{}& =\left(\begin{array}{c}\sqrt{\frac{\theta _1}{2\theta _1+\theta _2+\rho ^2}}|0,0\\ \sqrt{\frac{\zeta +\rho ^2}{2\theta _1+\theta _2+\rho ^2}}|1,0\\ 0\\ 0\end{array}\right)0,0|+\left(\begin{array}{c}0\\ |0,0\\ 0\\ 0\end{array}\right)1,0|\hfill \\ & +\sqrt{2\zeta +\rho ^2}\left(\begin{array}{c}\sqrt{\zeta +\rho ^2}\overline{z}_2\overline{z}_1\\ \overline{z}_2(Z_1+Z_2+\theta _2+\rho ^2)\\ \sqrt{\zeta +\rho ^2}\overline{z}_1\overline{z}_1\\ \overline{z}_1(Z_1+Z_2\theta _1)\end{array}\right)(ZZ_\rho ^{}{}_{}{}^{})^{1/2}t\hfill \end{array}$$
(117)
and $`U_2^{}`$ is the same as given in (112).

## 5 Discussion

To end this paper, we will explain here how our method can be used to deal with an interesting class of $`U(2)`$ multi-instantons. This case is quite similar to the elongated $`U(1)`$ multi-instanton.

We will take the following ADHM data:

$$I=\sqrt{k\zeta +\rho ^2}\left(\begin{array}{cc}e_k& 0\end{array}\right),J=\rho \left(\begin{array}{c}0\\ e_1^{}\end{array}\right),$$
(118)
$$B_1=\underset{i=1}{\overset{k1}{}}\sqrt{i\zeta +\rho ^2}e_ie_{i+1}^{},B_2=0,$$
(119)
where

$$e_i^{}=(\stackrel{1}{0},\mathrm{},\stackrel{i1}{0},\stackrel{i}{1},\stackrel{i+1}{0},\mathrm{},\stackrel{k}{0}).$$
(120)
One can easily check that the above data do form a $`U(2)`$ $`k`$-instanton solution, which includes the $`U(2)`$ 2-instanton solution in the previous section as a special case. It is the $`U(2)`$ analog of the elongated $`U(1)`$ $`k`$-instanton , and we may call it the ‘elongated $`U(2)`$ $`k`$-instanton’.

Because the above data always satisfy the condition $`[B_1,B_2]=0`$, it is easy to obtain

$$\mathrm{\Pi }_0=\left(\begin{array}{cc}_{n=0}^{k1}|n,0n,0|& 0\\ 0& 0\end{array}\right),$$
(121)
and to calculate $`\mathrm{\Gamma }`$ by (73). Then there is no difficulty to work out the instanton configuration.

## Acknowledgments

Chuan-Jie Zhu would like to thank Prof. Zhu-Jun Zheng and the hospitality at the Institute of Mathematics, Henan University.

## Appendix A ADHM Construction for Noncommutative ’t Hooft Instantons

In this appendix we apply our method to the construction of noncommutative ’t Hooft instantons . We will only discuss the ASD ’t Hooft instantons on ASD $`𝐑_{\mathrm{NC}}^\mathrm{𝟒}`$ (i.e., $`\zeta =\theta _1+\theta _2=0`$ while $`\theta _{1,2}0`$).

Take the following ADHM data :

$$a=\left(\begin{array}{c}\rho 1_2\\ a^n\overline{\sigma }_n\end{array}\right),\rho =(\rho _1,\mathrm{},\rho _k),a^n=\mathrm{diag}(a_1^n,\mathrm{},a_k^n),$$
(122)
where $`\rho _i`$ and $`a_i^n`$ are constants parameterizing the scale and the position of the $`i`$th instanton respectively. By direct calculation one finds

$$f=r^2r^2\rho ^{}\varphi ^1\rho r^2,r^2r^nr^n$$
(123)
where

$$r^n=x^na^n\mathrm{diag}(r_1^n,\mathrm{},r_k^n),\varphi =1+\rho _i^2r_i^2.$$
(124)
Summation of repeated indices is assumed here for both $`n`$ and $`i`$. It is straightforward to obtain

$$1KfK^{}=\varphi ^11_2,$$
(125)
$$\mathrm{\Delta }^{}fK^{}=\mathrm{\Delta }^{}r^2\rho ^{}\varphi ^1,$$
(126)
where a tensor product with $`1_2`$ is omitted in the second equation.

From the ADHM data (122) we see that the matrices $`B_1`$ and $`B_2`$ take very simple diagonal form. So we have also $`[B_1,B_2]=0`$ and the extra zero modes from (73) are just $`k`$ coherent states respectively shifted by the complex components of $`a_1^n,\mathrm{},a_k^n`$. Denote these (orthonormal) coherent states as $`|a_i`$, we will have

$$\mathrm{\Gamma }=\left(\begin{array}{c}\mathrm{diag}(|a_1,\mathrm{},|a_k)\\ 0_k\end{array}\right).$$
(127)
Then the required $`U`$ is again obtained by (68-72) where $`k^{}=k`$:

$$U=\left(\begin{array}{c}0\\ \mathrm{\Gamma }\end{array}\right)Y+\left(\begin{array}{c}1_2\\ \mathrm{\Delta }^{}r^2\rho ^{}\end{array}\right)\varphi ^{1/2}𝒮.$$
(128)


format: text

text:
# Inelastic Collision Rates of Trapped Metastable Hydrogen

## I Introduction

The ability to create large clouds of metastable ($`2S`$) H atoms in a magnetic trap makes it possible to study the rich collisional physics of cold metastable H. In an encounter between two metastables, the possible outcomes include Penning ionization, the formation of molecular ions, excitation transfer to short-lived $`2P`$ states, and hyperfine transitions . If the $`2S`$ atoms are generated from a background $`1S`$ gas, then several $`1S`$-$`2S`$ collision processes may also occur. Though these states of H are among the simplest atomic states in nature, an accurate description of their inelastic collisions at low energies remains a theoretical challenge .

At low temperatures, the rate of quenching collisions between metastables is small enough to allow dense ($`>10^{10}`$ cm-3) $`2S`$ clouds to exist for tens of milliseconds, which is an appreciable fraction of the 122 ms natural lifetime of the $`2S`$ state. Such a cloud can serve as an atom source for precision spectroscopy of transitions originating in the $`2S`$ state. In particular, the absolute frequencies of transitions from the $`2S`$ to higher-lying states can be combined with the well-known $`1S`$-$`2S`$ interval to simultaneously determine the $`1S`$ Lamb shift (a sensitive test of QED) and the Rydberg constant (which relates several more fundamental constants) . Taking advantage of the high signal-to-noise ratio possible in a cold, trapped sample, spectroscopy of metastable H may lead to a significant advance in accuracy for the determination of these quantities .

Cold metastable H is also potentially interesting for “quantum atom optics” experiments involving single-atom detection. The large (10 eV) internal energy of a $`2S`$ atom is easily registered by an ionizing collision on an electron-multiplying detector, or by detecting the Lyman-$`\alpha `$ photon emitted at a surface or other localized electric field. As demonstrated already with metastable noble gases, single-atom detection allows experiments involving atom-atom correlations, atom interferometry, or atom holography . A metastable beam ejected from a trapped H sample, potentially much brighter than other cold metastable beams to date, may be a useful source for quantum atom optics research.

After improving our apparatus at MIT, we have been able to generate metastable H clouds which are larger and longer-lived than in our previous work . In this paper, we report results from decay studies of trapped $`2S`$ clouds, conducted at temperatures ranging from a few mK to below 100 $`\mu `$K. These results include the first experimental determinations of the $`2S`$-$`2S`$ two-body loss rate constant $`K_2`$ and upper limits for the $`1S`$-$`2S`$ inelastic rate constant $`K_{12}`$. In the case of $`K_2`$, we have observed a temperature dependence below 230 $`\mu `$K which is not predicted in the current theory for cold metastable H collisions . This work serves as a stepping stone to high-resolution spectroscopy experiments which probe fundamental physics. Knowledge of metastable collision parameters may also be instrumental in the development of single-atom detection experiments involving H.

## II Experimental Methods

The methods we use for trapping and cooling atomic H in a cryogenic apparatus (Fig. 1(a)) are described elsewhere . Molecular hydrogen is dissociated in a discharge chamber that is thermally anchored at one end to a dilution refrigerator and opens to a trapping cell at the other end. Atoms in the $`F=1,m_F=1`$ hyperfine state of the $`1S`$ manifold are captured in a Ioffe-Pritchard magnetic trap , initially 500 mK deep. With typical aspect ratios between 100:1 and 400:1, the cloud is shaped like a highly elongated ellipsoid, cylindrically symmetric about the trap axis. We refer to the direction along the trap axis as “axial” and to the perpendicular direction as “radial.”

If an atom has sufficient energy, it can escape over a magnetic field saddlepoint located at one end of the trap. The sample can be cooled to a temperature of $`200\mu `$K by progressively lowering the saddlepoint field, thus forcing evaporation . If further cooling is desired, rf-induced evaporation is employed to reach temperatures as low as $`20\mu `$K . The number of atoms in the trap ranges from more than $`10^{14}`$ after initial loading to $`10^{10}`$ at the lowest temperatures. For the experiments described in this paper, the samples were not Bose degenerate, and typical peak $`1S`$ densities ranged from $`10^{13}`$ cm-3 to $`10^{14}`$ cm-3.

To generate metastable hydrogen, we excite the two-photon $`1S`$-$`2S`$ transition with a 243 nm laser . To achieve large excitation rates, typically $`15`$ mW of laser power is focused to a radius of $`20`$ $`\mu `$m near the trap minimum. In accordance with the two-photon selection rules, the metastables produced are in the $`F=1,m_F=1`$ state of the $`2S`$ manifold; except for a small relativistic correction $`(10^4)`$, these experience the same magnetic trapping potential as the $`1S`$ atoms. In the presence of a weak electric field of strength $`E`$, the $`2S`$ state quenches at a rate $`\gamma _s=2800E^2(`$cm$`{}_{}{}^{2}/`$V$`{}_{}{}^{2})`$s-1 due to Stark-mixing with the $`2P`$ state. To detect the metastables, a field of 10 V/cm is applied, quenching the $`2S`$ atoms in a few microseconds. A small fraction of the resulting 122 nm Lyman-$`\alpha `$ photons are counted on a microchannel plate (MCP) detector. The metastable decay behavior is observed by cycling through a series of different wait times between excitation and quench pulses (Fig. 1(b)). During the trap lifetime of the ground state sample, typically several minutes, we can make many $`2S`$ decay measurements. Samples of differing initial $`2S`$ density are obtained by stepping the laser frequency across the resonance.

For the measurements discussed here, the detection efficiency was approximately $`2\times 10^6`$. This was calibrated using a 230 $`\mu `$K sample where the initial number of $`1S`$ atoms was known from bolometric density measurements  and the trap geometry. The sample was depleted by resonant $`1S`$-$`2S`$ excitation and subsequent quenching to an untrapped state. With appropriate corrections, the detection efficiency can be inferred from the total number of signal photons counted and the initial number of atoms . The corrections account for other loss mechanisms, for the branching ratio back to the trapped ground state, and for reabsorption (“radiation trapping”) of the Lyman-$`\alpha `$ photons by the ground state cloud. The photons emitted by the quenched metastables have a frequency which is detuned from a 100 MHz-wide $`1S`$-$`2P`$ absorption resonance by only 1.1 GHz, and at the densities achieved in our experiments, the effects of radiation trapping are significant. Using Monte Carlo simulations of Lyman-$`\alpha `$ propagation in experimentally realistic cloud geometries, it was determined that typically $`30`$% of potential signal photons are absorbed and scattered in the atom cloud . Other factors, such as the small detection solid angle ($`1.1\times 10^2`$ sr), the transmission of two MgF2 windows between the atoms and the MCP, and the MCP quantum efficiency, lead to the small overall detection efficiency.

Based on the detection efficiency, we can estimate $`N_{2S}`$, the population of the metastable cloud, and $`n_{2S}`$, the peak metastable density in the trap. In a very cold sample (87 $`\mu `$K), as many as 150 Lyman-$`\alpha `$ photons were observed on the MCP after a single laser shot, implying $`N_{2S}8\times 10^7`$. The effective volume, $`V_{2S}=N_{2S}/n_{2S}`$, of the metastable clouds for these samples is $`1\times 10^3`$ cm3, which means that densities of $`10^{11}`$ cm-3 have been achieved.

Due to scatter from optics inside the trapping cell, each laser pulse induces some background fluorescence which is detected with low efficiency at the MCP. The background, which is typically small compared to the metastable signal, decays in a few milliseconds to a negligible level. To establish a background correction for our decay data, the laser is tuned far off resonance at the end of each trapping sequence. The background fluorescence is then recorded for the same cycle of wait times used when generating metastables.

In earlier experiments, performed in an entirely nonmetallic cell, we found that stray electric fields limited metastable lifetimes to 1 ms or less. In the current version of the apparatus, a copper film has been added to the inside surface of the trapping cell. The film is thin enough to allow rf to penetrate the trapping cell from outside coils, yet thick enough to significantly suppress stray dc fields. As depicted in Fig. 2(a), the copper film has been divided into several contiguous electrodes which are used to apply both the quench pulses and a compensating dc field. The dc voltage which best compensates for stray fields is found by measuring the metastable decay rate in a low density sample for several applied voltages. As shown in Fig. 2(b), we observe the expected quadratic dependence of the decay rate on electric field. With optimal stray field compensation, metastable lifetimes as long as 100 ms have been observed, approaching the natural lifetime of the $`2S`$ state.

Our decay measurements have focused on metastable clouds generated in four magnetic trap configurations, spanning a range of temperatures and ground state densities. These trap configurations are designated Traps W, X, Y, and Z, and each has its own standard evaporation sequence for sample preparation. To increase the signal-to-noise ratio, decay measurements were repeated in each trap over several consecutive trap cycles.

Densities and temperatures of the standard samples and the approximate trajectories by which they were created are depicted in Fig. 3 . To ensure that the colder samples (W, X, and Y) had reached thermal equilibrium, they were held undisturbed for 30 s after the end of the forced evaporation. The sample densities were measured to 20% accuracy by the bolometric method , and the temperatures were determined by a combination of methods. For the coldest sample (W), the temperature was calibrated to within 5% by measuring the width of the Doppler-sensitive $`1S`$-$`2S`$ line . The temperature of Trap X was determined by comparing the widths of the Doppler-free $`1S`$-$`2S`$ resonances of Traps W and X in the low density limit . For Traps Y and Z, the temperature was established using results from a numerical simulation of evaporative cooling . The temperature uncertainty for Traps X, Y, and Z is estimated to be 10%.

## III Evidence for $`\mathrm{𝟐}𝑺`$-$`\mathrm{𝟐}𝑺`$ Two-Body Loss

Our investigation of inelastic $`\mathrm{𝟐}𝑺`$-$`\mathrm{𝟐}𝑺`$ collisions was initiated by observations of how the decay behavior depended on metastable density. In each of the standard traps the $`\mathrm{𝟐}𝑺`$ density was varied by generating metastables at different laser detunings and also by exciting different $`\mathrm{𝟏}𝑺`$ densities as the ground state sample decayed. The decay curves were binned according to the total number of metastable signal counts in each decay curve. Since, to a first approximation, the signal is proportional to the metastable density, this means that decay curves corresponding to approximately the same initial $`\mathrm{𝟐}𝑺`$ density were averaged together. For a preliminary analysis, the average decay curves were fit to a simple exponential function, with only the decay rate and initial amplitude as free parameters. Figure 4 shows the results of this analysis. In the two coldest traps, W and X, the decay rate clearly increases with increasing metastable density, consistent with the presence of two-body inelastic loss processes. Furthermore, among the four traps, the decay rate is largest in Trap W, where the highest metastable densities are achieved.

The peak $`\mathrm{𝟐}𝑺`$ densities achieved in each trap can be estimated from the number of signal counts observed, but the uncertainty is large, reflecting significant uncertainties in both the metastable cloud shape and the detection efficiency. The shape of the $`\mathrm{𝟐}𝑺`$ cloud depends sensitively on the location of the laser beam waists with respect to the magnetic trap. For the density values listed in Fig. 4, a detection efficiency of $`\mathrm{𝟐}\mathbf{\times }\mathrm{𝟏𝟎}^\mathbf{}\mathrm{𝟔}`$ was assumed, and the counterpropagating laser beams were assumed to be perfectly overlapped and well-aligned with the atom cloud; with a plausible misalignment, the actual densities could be a factor of 3 smaller.

As illustrated in Fig. 5, a comparison of the average decay curves at high and low metastable signal provides further evidence for two-body loss at high densities. At low $`\mathrm{𝟐}𝑺`$ densities, the decay curves are well described by an exponential decay. By contrast, the decay curves at high $`\mathrm{𝟐}𝑺`$ densities deviate significantly from a simple exponential. They are better fit by a “one-plus-two” model with a free two-body decay parameter and a one-body decay rate fixed at its low density value. The one-plus-two model is described further in Sec. IV.

In principle, three-body inelastic collisions could also contribute to the observed density dependence of the decay. However, the three-body loss rate constant for metastables would have to be unexpectedly large to play a significant role at our experimental densities ($`\mathbf{}\mathrm{𝟏𝟎}^{\mathrm{𝟏𝟏}}`$ cm-3) . In the following analysis, we assume that two-body effects dominate.

## IV Determination of the Two-Body Loss Rate Constant

In this section, we derive the rate constant $`𝑲_\mathrm{𝟐}`$ for two-body loss of metastables based on decay data of the form $`𝒔\mathbf{(}𝒕\mathbf{)}`$, the number of signal counts observed when quenching at wait time $`𝒕`$. The constant is defined as follows to give the rate of change of the local $`\mathrm{𝟐}𝑺`$ density $`𝒏_{\mathrm{𝟐}𝑺}\mathbf{(}𝐫\mathbf{)}`$ due to two-body processes:

$$\mathbf{.}\dot{𝒏}_{\mathrm{𝟐}𝑺}\mathbf{(}𝐫\mathbf{)}\mathbf{|}_{\mathrm{𝐭𝐰𝐨}\mathbf{}\mathrm{𝐛𝐨𝐝𝐲}}\mathbf{=}\mathbf{}𝑲_\mathrm{𝟐}𝒏_{\mathrm{𝟐}𝑺}^\mathrm{𝟐}\mathbf{(}𝐫\mathbf{)}\mathbf{.}$$
(1)
If the relative distribution of metastables in the trap is not changing in time (“static approximation”), $`𝑲_\mathrm{𝟐}`$ is the product of the metastable detection efficiency, a factor describing the cloud geometry, and a two-body decay parameter extracted from fits to $`𝒔\mathbf{(}𝒕\mathbf{)}`$. After analyzing the decay process in the static approximation, we introduce dynamic corrections to obtain final results for $`𝑲_\mathrm{𝟐}`$.

Given the geometry of our apparatus and the ground state densities achieved, it is a reasonable first approximation to assume that the metastable cloud does not change shape as it decays. The initial extent of the metastable cloud along the trap axis, about 2 cm between the $`\mathrm{𝟏}\mathbf{/}𝒆`$ points, is defined by the depth of focus of the laser. This is shorter than the thermal length of the ground state cloud, which ranges from 4 cm in Trap W to 23 cm in Trap Z. On the time scale of the decay, the $`\mathrm{𝟐}𝑺`$ atoms collide frequently with $`\mathrm{𝟏}𝑺`$ atoms. At the peak ground state density in Trap W, for example, the $`\mathrm{𝟐}𝑺`$ mean free path is 100 $`𝝁`$m, and the mean time between collisions is 70 $`𝝁`$s. (To calculate these values, we use the theoretical value for $`𝒂_{\mathrm{𝟏}𝑺\mathbf{}\mathrm{𝟐}𝑺}`$, the elastic $`\mathrm{𝟏}𝑺`$-$`\mathrm{𝟐}𝑺`$ scattering length .) Thus, the metastables diffuse slowly along the trap axis. Numerical simulations of the time evolution of the $`\mathrm{𝟐}𝑺`$ cloud have shown that over a cloud lifetime of 100 ms, the shape of the initial axial distribution does not change dramatically . In the radial direction, the cloud shape evolves even less, although for a different reason. Due to the geometry of the magnetic trap, the metastable cloud is cylindrically symmetric and confined to the same thermal radius as the ground state cloud. This thermal radius, about 110 $`𝝁`$m for Trap W, is comparable to the $`\mathrm{𝟐}𝑺`$ mean free path. In a time shorter than the 2 ms excitation pulse, the metastables establish an equilibrium spatial distribution across the short radial dimension. The relative spatial distribution of metastables is quasi-static, then, because (1) axial diffusion through the $`\mathrm{𝟏}𝑺`$ gas is slow and (2) radial equilibrium exists over the lifetime of the $`\mathrm{𝟐}𝑺`$ cloud.

In the static approximation, the rate of signal decay can be described by the differential equation

$$\dot{𝒔}\mathbf{=}\mathbf{}𝜶_\mathrm{𝟏}𝒔\mathbf{}𝜶_\mathrm{𝟐}𝒔^\mathrm{𝟐}\mathbf{,}$$
(2)
where $`𝜶_\mathrm{𝟏}`$ and $`𝜶_\mathrm{𝟐}`$ are, respectively, parameters corresponding to one- and two-body decay. The solution to Eq. 2 is

$$𝒔\mathbf{(}𝒕\mathbf{)}\mathbf{=}\frac{𝜶_\mathrm{𝟏}𝑨_𝒐𝒆^{\mathbf{}𝜶_\mathrm{𝟏}𝒕}}{𝜶_\mathrm{𝟏}\mathbf{+}𝜶_\mathrm{𝟐}𝑨_𝒐\mathbf{(}\mathrm{𝟏}\mathbf{}𝒆^{\mathbf{}𝜶_\mathrm{𝟏}𝒕}\mathbf{)}}\mathbf{,}$$
(3)
which we have named the “one-plus-two” model. The parameter $`𝜶_\mathrm{𝟏}`$ is the one-body decay rate, the total rate of decay due to processes which are independent of $`\mathrm{𝟐}𝑺`$ density. Contributions to the one-body rate include the natural decay of the $`\mathrm{𝟐}𝑺`$ state, quenching due to stray electric fields, and inelastic collisions with the background $`\mathrm{𝟏}𝑺`$ gas.

For a decay measurement with a very high signal-to-noise ratio, one could extract $`𝜶_\mathrm{𝟏}`$, $`𝜶_\mathrm{𝟐}`$, and $`𝑨_𝒐`$ simultaneously by a fit to the “one-plus-two” model. The parameters $`𝜶_\mathrm{𝟏}`$ and $`𝜶_\mathrm{𝟐}`$ are highly correlated, however. With the modest signal-to-noise ratio in our data, a consistent convergence of the fit is not assured. Thus, we employ an alternate method of analysis. For a given trap and $`\mathrm{𝟏}𝑺`$ density, the value for $`𝜶_\mathrm{𝟏}`$ is determined by extrapolating the best-fitting exponential decay rates to zero metastable density (Fig. 6). This leaves two free parameters in the one-plus-two model: $`𝑨_𝒐`$, the signal which would be observed at $`𝒕\mathbf{=}\mathrm{𝟎}`$, and $`𝜶_\mathrm{𝟐}`$, from which $`𝑲_\mathrm{𝟐}`$ can be calculated.

By integrating Eq. 1 over the entire trap, we obtain an expression for the total number of metastables lost per unit time due to two-body processes:

$$\mathbf{.}\dot{𝑵}_{\mathrm{𝟐}𝑺}\mathbf{|}_{\mathrm{𝐭𝐰𝐨}\mathbf{}\mathrm{𝐛𝐨𝐝𝐲}}\mathbf{=}\mathbf{}𝑲_\mathrm{𝟐}𝒏_{𝒐\mathbf{,}\mathrm{𝟐}𝑺}^\mathrm{𝟐}\mathbf{}𝒇_{\mathrm{𝟐}𝑺}^\mathrm{𝟐}\mathbf{(}𝐫\mathbf{)}𝒅^\mathrm{𝟑}𝐫\mathbf{.}$$
(4)
Here, $`𝑵_{\mathrm{𝟐}𝑺}`$ is the number of $`\mathrm{𝟐}𝑺`$ atoms in the trap, and $`𝒇_{\mathrm{𝟐}𝑺}\mathbf{(}𝐫\mathbf{)}\mathbf{=}𝒏_{\mathrm{𝟐}𝑺}\mathbf{(}𝐫\mathbf{)}\mathbf{/}𝒏_{𝒐\mathbf{,}\mathrm{𝟐}𝑺}`$ is the normalized spatial distribution function for metastables, where $`𝒏_{𝒐\mathbf{,}\mathrm{𝟐}𝑺}`$ is the peak $`\mathrm{𝟐}𝑺`$ density. Since $`𝒔\mathbf{=}\mathit{ϵ}𝑵_{\mathrm{𝟐}𝑺}`$, where $`\mathit{ϵ}`$ is the detection efficiency, and since

$$\mathbf{.}\dot{𝒔}\mathbf{|}_{\mathrm{𝐭𝐰𝐨}\mathbf{}\mathrm{𝐛𝐨𝐝𝐲}}\mathbf{=}\mathbf{}𝜶_\mathrm{𝟐}𝒔^\mathrm{𝟐}\mathbf{,}$$
(5)
it follows that

$$𝑲_\mathrm{𝟐}\mathbf{=}𝜻\mathit{ϵ}𝜶_\mathrm{𝟐}$$
(6)
in the static approximation, where $`𝜻`$ is a geometry factor defined by

$$𝜻\mathbf{=}\frac{\mathbf{}𝒇_{\mathrm{𝟐}𝑺}^\mathrm{𝟐}\mathbf{(}𝐫\mathbf{)}𝒅^\mathrm{𝟑}𝐫}{\mathbf{\left[}\mathbf{}𝒇_{\mathrm{𝟐}𝑺}\mathbf{(}𝐫\mathbf{)}𝒅^\mathrm{𝟑}𝐫\mathbf{\right]}^\mathrm{𝟐}}\mathbf{.}$$
(7)
To obtain Eq. 6, we have used the fact that $`𝒏_{𝒐\mathbf{,}\mathrm{𝟐}𝑺}\mathbf{=}𝑵_{\mathrm{𝟐}𝑺}\mathbf{/}\mathbf{}𝒇_{\mathrm{𝟐}𝑺}\mathbf{(}𝐫\mathbf{)}𝒅^\mathrm{𝟑}𝐫`$. The quantity $`𝜻`$, which has units of volume, is a measure of the spatial extent of the metastable cloud; $`𝜻`$ decreases as $`𝒇_{\mathrm{𝟐}𝑺}\mathbf{(}𝐫\mathbf{)}`$ becomes more sharply peaked.

Thus to determine $`𝑲_\mathrm{𝟐}`$ from our decay data it is necessary to know $`\mathit{ϵ}`$ and the shape of the metastable cloud. Our method for determining $`\mathit{ϵ}`$ was described in Sec. II. To calculate $`𝜻`$ for various trap configurations and temperatures, we employed a Monte Carlo simulation of $`\mathrm{𝟏}𝑺`$-$`\mathrm{𝟐}𝑺`$ excitation in our trap . The simulation calculates the excitation lineshape associated with atoms in a volume slice perpendicular to the trap axis. In addition to the position of the volume slice along the axis, the inputs to the simulation include the laser field geometry, the laser pulse length, the trap shape, the sample temperature, and the peak $`\mathrm{𝟏}𝑺`$ density. The lineshape is calculated by choosing atoms randomly from a thermal distribution, computing their trajectories during the laser pulse, and then finding the contribution of each atom to the spectrum by Fourier-transforming the time-varying laser field amplitude seen by the atom. The effects of the $`\mathrm{𝟏}𝑺`$-$`\mathrm{𝟐}𝑺`$ cold collision shift  are incorporated by allowing the resonance frequency to vary over the trajectory in accordance with the local $`\mathrm{𝟏}𝑺`$ density. After repeating the lineshape calculation for many slices covering the length of the $`\mathrm{𝟐}𝑺`$ cloud along the trap axis, the axial distribution of metastables at $`𝒕\mathbf{=}\mathrm{𝟎}`$ can be extracted for any given laser detuning. Under the assumption that the metastables establish an equilibrium radial distribution during the pulse, $`𝒇_{\mathrm{𝟐}𝑺}\mathbf{(}𝐫\mathbf{)}`$ for the entire trap is easily calculated from the axial distribution.

An additional numerical calculation corrects the lineshape simulation results to account for photoionization, the promotion of $`\mathrm{𝟐}𝑺`$ electrons to the continuum by the excitation laser. In a cold, dense sample, such as those prepared in Trap W, the metastables spend a relatively long time in the high-intensity UV laser field, and the fraction of metastables lost to photoionization can exceed 30% .

Because the lineshape simulation has given results in good agreement with experimental $`\mathrm{𝟏}𝑺`$-$`\mathrm{𝟐}𝑺`$ spectra for a wide range of trap shapes, densities, and temperatures, this method of calculating $`𝒇_{\mathrm{𝟐}𝑺}\mathbf{(}𝐫\mathbf{)}`$ is believed reliable. For each spatial distribution function, $`𝜻`$ is computed by numerical evaluation of Eq. 7. Values for $`𝜻`$ range from approximately $`\mathrm{𝟐}\mathbf{\times }\mathrm{𝟏𝟎}^\mathbf{}\mathrm{𝟑}`$ cm3 in Trap W to 0.2 cm3 in Trap Z. Depending on the laser detuning, an imperfect overlap of the counterpropagating laser beam focuses along the trap axis may cause $`𝜻`$ to be 50-100% larger. Uncertainty in the knowledge of the laser field geometry is the principal source of uncertainty in $`𝜻`$.

The experimental results for $`𝑲_\mathrm{𝟐}`$, calculated on the basis of Eq. 6, are listed in the “static approximation” column of Table I. For each standard trap, the accumulated decay curves were binned together according to laser detuning, and the resulting average curves were fit with the one-plus-two model to obtain $`𝜶_\mathrm{𝟐}`$ as a function of detuning. Next, the values for $`𝜶_\mathrm{𝟐}`$ were multiplied by the detection efficiency and the value of $`𝜻`$ appropriate for the detuning. The resulting values for $`𝑲_\mathrm{𝟐}`$ were found to be largely independent of laser detuning, indicating that variations in metastable cloud shape with detuning are reasonably accounted for by our lineshape simulation. For Traps W and X, the weighted averages across all detunings are reported in Table I. In the warmer samples, where the $`\mathrm{𝟐}𝑺`$ density is too low to observe two-body collisions, the weighted average for $`𝑲_\mathrm{𝟐}`$ is consistent with zero. For these, Table I reports upper limits for $`𝑲_\mathrm{𝟐}`$ based on linear addition of the uncertainties in $`𝜶_\mathrm{𝟐}`$, $`\mathit{ϵ}`$, and $`𝜻`$.

The error bars for Traps W and X are primarily due to systematic uncertainties. These errors are asymmetric due to asymmetric uncertainty in the detection efficiency and the fact that laser misalignments generally increase $`𝜻`$. The most important contributions to the overall error in $`𝑲_\mathrm{𝟐}`$ come from the detection efficiency and uncertainty in the axial overlap of the laser focuses. In Trap W, for example, $`\mathit{ϵ}\mathbf{=}\mathbf{1.8}\mathbf{\times }\mathrm{𝟏𝟎}^\mathbf{}\mathrm{𝟔}`$ with approximate relative uncertainties of 30% on the lower side and 60% on the upper side. The focus overlap uncertainty contributes 75% to the upper relative error. Other sources of uncertainty include the error in the determination of $`𝜶_\mathrm{𝟐}`$ $`\mathbf{(}\mathbf{\pm }\mathrm{𝟏𝟐}\mathbf{\%}\mathbf{)}`$ and contributions of a few percent each due to uncertainties in the beam radius, cold collision shift parameter, offset of the trap minimum and laser focuses, ground state density, and sample temperature. With some minor exceptions, the various uncertainties are assumed to be uncorrelated, and the relative errors we report for $`𝑲_\mathrm{𝟐}`$ are essentially quadrature sums of the contributing errors.

The static approximation results need to be refined to account for change in the metastable spatial distribution while the cloud decays. There is a flattening of the initial axial distribution due to diffusion and to two-body loss, which occurs preferentially where the $`\mathrm{𝟐}𝑺`$ density is highest. This flattening of the distribution in time causes the initial decay of the metastable population to be somewhat less steep than it would be if the cloud shape were completely static. Thus, the value of $`𝜶_\mathrm{𝟐}`$ determined by fitting a decay curve to the one-plus-two model is smaller than the effective value of $`𝜶_\mathrm{𝟐}`$ which accounts for dynamics when substituted into Eq. 6.

To calculate corrections to the static approximation analysis, a numerical simulation of the decaying metastable cloud was developed. For each standard sample, the simulation was repeated many times assuming different laser detunings. The starting point for each simulation run was a cloud with distribution $`𝒇_{\mathrm{𝟐}𝑺}\mathbf{(}𝐫\mathbf{)}`$, calculated by the method described above, and an initial total number of metastables $`𝑵_{\mathrm{𝟐}𝑺}\mathbf{(}\mathrm{𝟎}\mathbf{)}\mathbf{=}𝑨_𝒐\mathbf{/}\mathit{ϵ}`$, where $`𝑨_𝒐\mathbf{=}𝒔\mathbf{(}\mathrm{𝟎}\mathbf{)}`$ was determined from experimental data. The simulation then evolved both dynamic and static clouds forward in time from the initial condition. For the dynamic cloud, the axial distribution was adjusted locally at each time step for losses due to one-body processes (parametrized by $`𝜶_\mathrm{𝟏}`$), two-body processes (parametrized by a putative value for $`𝑲_\mathrm{𝟐}`$), and for diffusion of $`\mathrm{𝟐}𝑺`$ atoms in the $`\mathrm{𝟏}𝑺`$ background. For the static cloud, $`𝑵_{\mathrm{𝟐}𝑺}`$ was adjusted at each time step for one- and two-body losses integrated over the whole cloud, but the spatial distribution remained fixed. After evolving the clouds over 100 ms, the simulated static and dynamic metastable decay curves were fit by the one-plus-two model. The ratio $`𝜶_{\mathrm{𝟐}\mathbf{,}\mathrm{𝐬𝐭𝐚𝐭𝐢𝐜}}\mathbf{/}𝜶_{\mathrm{𝟐}\mathbf{,}\mathrm{𝐝𝐲𝐧𝐚𝐦𝐢𝐜}}`$ from the simulation was used as a correction factor for the experimentally determined $`𝜶_\mathrm{𝟐}`$ values. New average values for $`𝑲_\mathrm{𝟐}`$ were then calculated for each trap. It was found that the static approximation values of $`𝑲_\mathrm{𝟐}`$ must be multiplied by 1.35 and 1.30, respectively, to correct for dynamics in Traps W and X. By varying the inputs to the decay simulation, the uncertainty in the overall correction factors was estimated to be 10%. For Traps Y and Z, a conservative upper limit for the correction factor is 2. The final results for $`𝑲_\mathrm{𝟐}`$ are summarized in the last column of Table I.

Theoretical $`\mathrm{𝟐}𝑺`$-$`\mathrm{𝟐}𝑺`$ collision rates at low temperatures were calculated recently by Forrey and collaborators , extending their previous work for high temperatures . The calculations indicate that the most prevalent inelastic process in a cold metastable cloud is double excitation transfer, $`𝐇\mathbf{(}\mathrm{𝟐}𝑺\mathbf{)}\mathbf{+}𝐇\mathbf{(}\mathrm{𝟐}𝑺\mathbf{)}\mathbf{}𝐇\mathbf{(}\mathrm{𝟐}𝑷\mathbf{)}\mathbf{+}𝐇\mathbf{(}\mathrm{𝟐}𝑷\mathbf{)}`$. Next in importance are the two possible ionization processes: associative ionization, in which the molecular ion $`𝐇_\mathrm{𝟐}^\mathbf{+}`$ is formed, and Penning ionization, where the internal energy of one metastable atom causes the ionization of a second metastable. Since processes involving $`\mathrm{𝟐}𝑺`$ hyperfine state changes are mediated by comparatively weak magnetic dipole interactions, they are expected to be negligible in relation to excitation transfer and ionization. Thus, a theoretical estimate for $`𝑲_\mathrm{𝟐}`$ is obtained by summing the loss rate constants for double excitation transfer and ionization. For each process, the loss rate constant equals $`\mathrm{𝟐}\mathbf{}𝝈𝒗\mathbf{}`$, where $`𝝈`$ is the cross section for the process, $`𝒗`$ is the relative velocity of two atoms, the angle brackets indicate a thermal average, and the factor of 2 accounts for the fact that two metastables are lost in each collision. The theoretical curve for $`𝑲_\mathrm{𝟐}`$ as a function of temperature is shown alongside the experimental values in Fig. 7.

The uncertainty in the theoretical results is difficult to estimate. In its current state of development, the theory for $`\mathrm{𝟐}𝑺`$-$`\mathrm{𝟐}𝑺`$ collisions neglects hyperfine structure and assumes zero coupling between the relative orbital angular momentum of the nuclei and the angular momentum of the electrons. It is not yet certain whether this assumption is justified for spin-polarized metastables at low temperatures . Furthermore, the theory calculations assume zero magnetic field, while in the experiment the trapped metastables experience a field of a few gauss. Nevertheless, from Fig. 7 we can conclude that the present theory correctly predicts $`𝑲_\mathrm{𝟐}`$ at low temperatures to within an order of magnitude.

An intriguing aspect of the experimental results for $`𝑲_\mathrm{𝟐}`$ is the apparent temperature dependence between 87 $`𝝁`$K and 230 $`𝝁`$K. Although the absolute uncertainty is large for each temperature point, we believe that the temperature dependence is real. The error bars are primarily due to systematic uncertainties which affect both points in approximately the same way. When the correlated errors are removed, the ratio of $`𝑲_\mathrm{𝟐}`$ at 87 $`𝝁`$K to $`𝑲_\mathrm{𝟐}`$ at 230 $`𝝁`$K is found to be $`\mathbf{1.9}\mathbf{\pm }\mathbf{0.4}`$. This result contrasts with the theoretical prediction that $`𝑲_\mathrm{𝟐}`$ is virtually independent of temperature in this regime, consistent with Wigner threshold behavior. Further theoretical and experimental work will be necessary to understand the temperature dependence of inelastic $`\mathrm{𝟐}𝑺`$-$`\mathrm{𝟐}𝑺`$ collisions near $`𝑻\mathbf{=}\mathrm{𝟎}`$.

## V Upper Limit for the $`\mathrm{𝟏}𝑺`$-$`\mathrm{𝟐}𝑺`$ Loss Rate Constant

In addition to the two-body loss channels discussed above, metastables are lost from the trap by several one-body mechanisms. These include the natural decay of the $`\mathrm{𝟐}𝑺`$ state by spontaneous emission of two photons, quenching due to stray electric fields, and loss due to inelastic $`\mathrm{𝟏}𝑺`$-$`\mathrm{𝟐}𝑺`$ collisions. In a collision with a ground state atom, a metastable may be lost by a transition to either an untrapped $`\mathrm{𝟐}𝑺`$ state (hyperfine-changing collision) or to a $`\mathrm{𝟐}𝑷`$ state (excitation transfer). If inelastic $`\mathrm{𝟏}𝑺`$-$`\mathrm{𝟐}𝑺`$ collisions are sufficiently probable, they should result in a dependence of the observed one-body decay parameter $`𝜶_\mathrm{𝟏}`$ on the peak ground state density $`𝒏_{\mathrm{𝟏}𝑺\mathbf{,}𝒐}`$.

To search for evidence of inelastic $`\mathrm{𝟏}𝑺`$-$`\mathrm{𝟐}𝑺`$ collisions, the decay curves taken in the standard trap configurations were organized into consecutive time bins of 9.6 s duration. Next, a value for $`𝜶_\mathrm{𝟏}`$ was determined for each bin by extrapolating the decay rate to zero metastable signal, as in Fig. 6. Then, for each trap, $`𝜶_\mathrm{𝟏}`$ was plotted as a function of time over the 120 s period following initial preparation of the sample (see Fig. 8(a), for example). During this time, the ground state density decays significantly due to $`\mathrm{𝟏}𝑺`$-$`\mathrm{𝟏}𝑺`$ dipolar decay (Fig. 8(b)).

As shown in Fig. 8, the data suggests that some component of one-body loss decreases as the ground state density decreases. However, since the trend is not large relative to uncertainties, it is inconclusive whether significant loss due to $`\mathrm{𝟏}𝑺`$-$`\mathrm{𝟐}𝑺`$ collisions exists in our metastable clouds. We have not attempted to fit the time (or density) dependence of $`𝜶_\mathrm{𝟏}`$ to a model but instead have determined conservative upper limits for the metastable loss rate due to collisions with the ground state background.

An upper limit on the total inelastic contribution $`𝜶_{\mathrm{𝟏}\mathbf{,}\mathrm{𝐢𝐧𝐞𝐥}}`$ to the one-body decay rate is found by subtracting the natural decay rate $`\mathbf{(}𝜸_{\mathrm{𝐧𝐚𝐭}}\mathbf{=}\mathbf{8.2}`$ s$`{}_{}{}^{\mathbf{}\mathrm{𝟏}}\mathbf{)}`$ from the total one-body decay rate:

$$𝜶_{\mathrm{𝟏}\mathbf{,}\mathrm{𝐢𝐧𝐞𝐥}}\mathbf{<}𝜶_\mathrm{𝟏}\mathbf{}𝜸_{\mathrm{𝐧𝐚𝐭}}\mathbf{.}$$
(8)
Analogous to Eq. 1, we define a rate constant $`𝑲_{\mathrm{𝟏𝟐}}`$ for the metastable loss rate due to inelastic $`\mathrm{𝟏}𝑺`$-$`\mathrm{𝟐}𝑺`$ collisions, such that

$$\mathbf{.}\dot{𝒏}_{\mathrm{𝟐}𝑺}\mathbf{(}𝐫\mathbf{)}\mathbf{|}_{\mathrm{𝟏}𝑺\mathbf{}\mathrm{𝟐}𝑺}\mathbf{=}\mathbf{}𝑲_{\mathrm{𝟏𝟐}}𝒏_{\mathrm{𝟏}𝑺}\mathbf{(}𝐫\mathbf{)}𝒏_{\mathrm{𝟐}𝑺}\mathbf{(}𝐫\mathbf{)}\mathbf{.}$$
(9)
Using a static approximation for the metastable cloud, Eq. 9 can be integrated over the entire trap to find the relation between $`𝑲_{\mathrm{𝟏𝟐}}`$ and $`𝜶_{\mathrm{𝟏}\mathbf{,}\mathrm{𝐢𝐧𝐞𝐥}}`$:

$$𝑲_{\mathrm{𝟏𝟐}}\mathbf{=}\frac{𝜶_{\mathrm{𝟏}\mathbf{,}\mathrm{𝐢𝐧𝐞𝐥}}𝑽_{\mathrm{𝟐}𝑺}}{𝒏_{\mathrm{𝟏}𝑺\mathbf{,}𝒐}𝑸_{\mathrm{𝟏𝟐}}}\mathbf{,}$$
(10)
where $`𝑽_{\mathrm{𝟐}𝑺}\mathbf{=}\mathbf{}𝒇_{\mathrm{𝟐}𝑺}\mathbf{(}𝐫\mathbf{)}𝒅^\mathrm{𝟑}𝐫`$, and $`𝑸_{\mathrm{𝟏𝟐}}\mathbf{=}\mathbf{}𝒇_{\mathrm{𝟏}𝑺}\mathbf{(}𝐫\mathbf{)}𝒇_{\mathrm{𝟐}𝑺}\mathbf{(}𝐫\mathbf{)}𝒅^\mathrm{𝟑}𝐫`$. By combining an experimental upper limit for the initial value of $`𝜶_{\mathrm{𝟏}\mathbf{,}\mathrm{𝐢𝐧𝐞𝐥}}`$ in each trap, a numerically computed upper limit for the ratio $`𝑽_{\mathrm{𝟐}𝑺}\mathbf{/}𝑸_{\mathrm{𝟏𝟐}}`$, and the measured peak density, we have obtained the upper limits on $`𝑲_{\mathrm{𝟏𝟐}}`$ listed in Table II.

These upper limits are much larger than the theoretical rate constant for $`\mathrm{𝟏}𝑺`$-$`\mathrm{𝟏}𝑺`$ dipolar decay, $`𝒈\mathbf{=}\mathbf{1.2}\mathbf{\times }\mathrm{𝟏𝟎}^{\mathbf{}\mathrm{𝟏𝟓}}`$ cm3/s . Since the total rate constant for $`\mathrm{𝟏}𝑺`$-$`\mathrm{𝟐}𝑺`$ hyperfine-changing collisions is expected to be comparable to $`𝒈`$, our results admit the possibility that $`\mathrm{𝟏}𝑺`$-$`\mathrm{𝟐}𝑺`$ excitation transfer collisions are much more probable than hyperfine-changing collisions. A first theoretical calculation of the relevant cross sections in this temperature regime is currently in progress .

## VI Conclusion

In summary, we have completed the first determinations of loss rates due to inelastic collisions in a cold, trapped metastable H gas coexisting with a ground state background. The experimental values for the two-body loss rate $`𝑲_\mathrm{𝟐}`$ agree to order of magnitude $`\mathbf{(}\mathbf{}\mathrm{𝟏𝟎}^\mathbf{}\mathrm{𝟗}`$ cm3/s) with the present theory for $`\mathrm{𝟐}𝑺`$-$`\mathrm{𝟐}𝑺`$ collisions. Further refinements to the theory will be necessary to explain an apparent increase in $`𝑲_\mathrm{𝟐}`$ with decreasing temperature below 230 $`𝝁`$K. Our metastable decay measurements also show that the total rate constant for inelastic $`\mathrm{𝟏}𝑺`$-$`\mathrm{𝟐}𝑺`$ collisions is not larger than $`\mathbf{}\mathrm{𝟏𝟎}^{\mathbf{}\mathrm{𝟏𝟑}}`$ cm3/s at the temperatures probed below 1 mK.

The present results pertain to total $`\mathrm{𝟐}𝑺`$-$`\mathrm{𝟐}𝑺`$ and $`\mathrm{𝟏}𝑺`$-$`\mathrm{𝟐}𝑺`$ inelastic collision rates. In future experiments with an enhanced Lyman-$`𝜶`$ detection efficiency and more complete suppression of stray electric fields, it may be possible to measure the rates due to excitation transfer processes alone by directly observing the fluorescence of $`\mathrm{𝟐}𝑷`$ states generated in collisions.

With the ability to create clouds of $`\mathrm{𝟏𝟎}^\mathrm{𝟕}`$-$`\mathrm{𝟏𝟎}^\mathrm{𝟖}`$ trapped metastables and a preliminary understanding of their inelastic collisions, the prospects are bright for fruitful new spectroscopic experiments involving the $`\mathrm{𝟐}𝑺`$ state. Recently, our group has performed spectroscopy of the H $`\mathrm{𝟐}𝑺`$-$`\mathrm{𝟒}𝑷`$ transition in a magnetic trap for the first time. Absorption spectroscopy on this and other single-photon transitions provide new tools for studying interactions of hydrogen states at low temperatures. It will also soon be possible to excite narrow two-photon transitions from the $`\mathrm{𝟐}𝑺`$ state, setting the stage for precision frequency measurements which probe fundamental physics.

## Acknowledgments

We are grateful to Alexander Dalgarno, Robert C. Forrey, Dale G. Fried, Piotr Froelich, and Thomas C. Killian for helpful conversations. We also thank Walter Joffrain for his assistance in the laboratory. This research was supported by the National Science Foundation and the Office of Naval Research.


text:
# Inhomogeneous vortex matter

## Abstract

We present a generalization of the continuum theory of vortex matter for non-uniform superfluid density. This theory explains the striking regularity of vortex lattices observed in Bose-Einstein condensates, and predicts the frequencies of long-wavelength lattice excitations.

Dense lattices of quantized vortices in rotating Bose-Einstein condensates (BECs) MITSci ; JILAept  are strikingly more regular than finite vortex arrays in homogeneous superfluid CampbellZiff  (see Fig. 1), even though BEC densities vary greatly over the sample. This Letter generalizes the Feynman-Tkachenko Feynman ; T1  continuum theory of ‘vortex matter’ to cases in which the condensate density varies slowly on the scale of the lattice spacing. This theory explains the lattices’ surprising regularity, and find pronounced effects of nonuniform density on lattice excitations.

We consider a two-dimensional regular array of vortices, whether realized in a very oblate BEC, or as parallel vortex lines in a prolate one. Denoting the lattice length scale by $`b`$ and introducing dimensionless complex co-ordinates $`bz=x+iy`$, a regular lattice has the positions $`z_{jk}`$ of parallel vortex lines given by $`z_{jk}=z_{jk}^0k+\tau j`$ for $`\tau =\tau _1+i\tau _2`$ ($`\tau _2>0`$ and $`\tau _i`$ real). Much is known about vortex lattices in superfluids of constant density $`\rho `$, simply from incompressible hydrodynamics; and it is all simplified by using dimensionless variables, expressing time, velocity, and energy in lattice units $`Mb^2/\mathrm{}`$, $`\mathrm{}/(Mb)`$, $`\mathrm{}^2/(Mb^2)`$ respectively, for $`M`$ the mass of the particles composing the superfluid. The regular triangular case $`\tau =(1+i\sqrt{3})/2`$ is the ground state T1  of a sample rotating at dimensionless rate $`\mathrm{\Omega }=\pi /\tau _2`$. The irrotational velocity field $`vv_x+iv_y`$ consists of a fine field $`v_f`$, which is periodic on the lattice scale, plus a coarse field $`v_c`$ obeying Feynman’s criterion Feynman

$$_zv_c_{\overline{z}}\overline{v}_c=2\pi i\rho _V$$
(1)
where for a regular lattice the vortex density $`\rho _V`$ is $`1/\tau _2`$. (In our complex notation 2$`_zA=b(\stackrel{}{}\stackrel{}{A}+i\stackrel{}{}\times \stackrel{}{A})`$ for any $`A`$. We assume counter-clockwise rotation.) And long-wavelength excitations of the lattice, $`z_{jk}(t)z_{jk}^0+D(z_{jk}^0,\overline{z}_{jk}^0;t)`$ obey the wave equations T1

$`i_t\left(_zD_{\overline{z}}\overline{D}\right)`$ $`=`$ $`2\mathrm{\Omega }\left(_zD+_{\overline{z}}\overline{D}\right)`$ (2)
$`i_t\left(_zD+_{\overline{z}}\overline{D}\right)`$ $`=`$ $`{\displaystyle \frac{1}{2}}_{\overline{z}z}\left(_zD_{\overline{z}}\overline{D}\right).`$ (3)
(Extensions to three-dimensional vortex matter Baym ; Fetter  become considerably more complicated.)

Although current dilute gaseous BECs are compressible fluids governed by the Gross-Pitaevskii equation (GPE) gp , much of BEC vortex physics can be cast into a simpler hydrodynamic form. In our dimensionless variables, the GPE in a frame rotating about the co-ordinate origin may be written

$`_t\rho `$ $`=`$ $`\mathrm{\Omega }_\varphi \rho \left[_z\left(\rho v\right)+_{\overline{z}}\left(\rho \overline{v}\right)\right]`$ (4)
$`_tv`$ $`=`$ $`_{\overline{z}}\left[|vi\mathrm{\Omega }z|^2+2V+2g\rho {\displaystyle \frac{_{\overline{z}z}\sqrt{\rho }}{\sqrt{\rho }}}\right]`$ (5)
$`V`$ $``$ $`V_{trap}{\displaystyle \frac{\mathrm{\Omega }^2}{2}}|z|^2`$
where $`v`$ is still the lab-frame velocity, $`\varphi `$ is the polar angular co-ordinate, $`V_{trap}(z,\overline{z})`$ is the trap potential, and $`g`$ is the dimensionless 2D coupling, determined by atomic and trap parameters. Except very near vortex cores, the rotating-frame velocity $`vi\mathrm{\Omega }z`$ is of order unity; and $`\rho ^{1/2}_{\overline{z}z}\sqrt{\rho }`$ is no larger, except in cores and sample edges, which may be treated separately as boundary layers. In current experiments the healing length $`\xi \left(g\rho \right)^{1/2},`$ which sets the vortex core size, is everywhere else much smaller than $`b`$.

So, outside vortex cores, the leading order results in the co-rotating frame are

$`\rho `$ $`=`$ $`\text{const.}{\displaystyle \frac{V}{g}}`$ (6)
$`\mathrm{\Omega }_\varphi \rho `$ $`=`$ $`_z\left(\rho v\right)+_{\overline{z}}\left(\rho \overline{v}\right).`$ (7)
Thus $`\rho `$ depends on $`v`$ only through the centrifugal modification of the trap potential. The velocity field, $`v`$ is determined by (7), plus the condition of irrotationality except at the quantized vortex cores:

$$_zv_{\overline{z}}\overline{v}=2\pi i\underset{jk}{}\delta ^2\left(zz_{jk}\right).$$
(8)
This, with $`\stackrel{}{}\stackrel{}{v}=0`$ instead of (7), is the starting point for Ref. T1 . So for constant $`V`$, and hence uniform $`\rho `$, the results of the incompressible case also apply to BECs.

If $`\rho `$ varies slowly in space, will not inhomogeneous effects be small? Not obviously: like $`\rho `$ itself, the lattice shape might vary slowly, but change greatly over the whole sample. Indeed, even for constant $`\rho `$ in the ground state of a finite vortex array, Campbell and Ziff CampbellZiff  found gradual but cumulatively large distortions. But there is a basic problem in extending their analysis to non-uniform $`\rho `$.

Investigations of homogeneous vortex matter have generally relied on the exact single-vortex solution $`v=v_1(\overline{z}\overline{z}_{jk})`$ to (7) and (8) for constant $`\rho `$, which one may simply sum over the vortex labels $`j,k`$, because (7) and (8) are linear in $`v`$. For non-constant $`\rho `$ the familiar $`v_1(z)=i/\overline{z}`$ satisfies (8) but not (7), and so we do not have the exact single-vortex solution. Perturbative approximations about $`v_1(z)=i/\overline{z}`$ as an ansatz break down at distances from the core beyond the length scale of the density variation JRA . So the few-vortex problem in inhomogeneous BECs becomes analytically intractable. For sufficiently dense lattices, however, inhomogeneous vortex matter yields to a different approach.

By a ‘dense’ vortex lattice, we mean $`\rho =\rho (\epsilon z,\epsilon \overline{z})`$ for small $`\epsilon `$. (For a round harmonic trap with Thomas-Fermi radius $`R`$, $`\epsilon =b/R`$, giving $`\epsilon 0.1`$ in current experiments.) We can therefore perturb in $`\epsilon `$; but to distinguish smallness from slowness, we must use multiple scale analysis (MSA) MSA . This formalism produces coarse-scale equations of motion for $`D`$, from which all lattice-scale physics has been eliminated, in the same sense that high-frequencies are eliminated by adiabatic methods. These will be our generalizations of (7) and (8).

The application of MSA leads to a rather involved derivation the details of which will be reported elsewhere; here we outline its steps and report its conclusions. We begin by satisfying (7) identically by defining

$$v=i\mathrm{\Omega }z+\frac{i}{\rho }_{\overline{z}}\left(\rho F\right)$$
(9)
for real $`F`$. We use vortex-centered co-ordinates:

$$z=z^{}+D(z^{},\overline{z}^{},t),$$
(10)
so that $`z_{jk}^{}=z_{jk}^0`$ regardless of $`D`$. We do this so that in the $`z^{}`$ co-ordinates we always have a regular lattice, whose symmetries we can exploit, even though the physical lattice may be distorted by excitations, or by a static $`D`$ field induced by inhomogeneous $`\rho `$. In the non-Cartesian $`z^{}`$ co-ordinates, Eqn. (8) becomes

$$\begin{array}{c}\mathrm{Re}\left(_z^{}\frac{_{\overline{z}^{}}2\left(_{\overline{z}^{}}D\right)_z^{}}{\rho }\rho F\right)=\pi \underset{jk}{}\delta ^2\left(z^{}z_{jk}^0\right)\hfill \\ \hfill \mathrm{\Omega }\left(1+_z^{}D+_{\overline{z}^{}}\overline{D}+\left|\begin{array}{cc}_z^{}D& _z^{}\overline{D}\\ _{\overline{z}^{}}D& _{\overline{z}^{}}\overline{D}\end{array}\right|\right)\end{array}$$
(11)
This shows explicitly that only gradients in $`D`$ affect $`F`$.

MSA then embeds the physical $`z^{}`$-plane in a fictitious 4-space of complex co-ordinates $`\zeta ,Z`$, as the subspace $`(\zeta ,Z)=(z^{},\epsilon z^{})`$. This provides $`_z^{}_\zeta +\epsilon _Z`$, etc., and proceeding perturbatively in $`\epsilon `$, we are able to write explicit solutions (in terms of rapidly converging series) for the $`\zeta `$-dependence of $`F`$ at every order (we need to go to third). As usual with MSA, the ‘gauge freedom’ in how functions depend explicitly on the two extra dimensions is used to remove solutions growing secularly with $`\zeta `$, by constraining the purely $`Z`$-dependent part of $`F`$ (which we denote by $`F_c`$). Once we restrict back to physical two-space by setting $`Z=\epsilon z^{},`$ $`\zeta =z^{}`$, and return from vortex-fixed $`z^{}`$ to Cartesian $`z`$, we recognize the constraint on $`F_c`$ as just what is needed to maintain Feynman’s condition (1), for $`\rho _V`$ as perturbed by $`D`$:

$$_z\frac{_{\overline{z}}\left(\rho F_c\right)}{\rho }+_{\overline{z}}\frac{_z\left(\rho F_c\right)}{\rho }=2\mathrm{\Omega }\left(_zD+_{\overline{z}}\overline{D}\right),$$
(12)
where we drop determinant terms quadratic in $`D`$ (because we will have $`D`$ order $`\epsilon `$ or smaller).

Having solved for $`F,`$ and hence $`v,`$ in terms of explicit lattice-periodic functions and $`v_c`$, we know the local fluid velocity near each vortex. This fixes the instantaneous vortex translational velocity field $`\dot{D}`$; but the fixing is not trivial. Since the hydrodynamic approximation to the GPE breaks down within $`|zz_{jk}|`$ $`\xi /b`$, in these small regions we must solve the time-dependent GPE using a different perturbation theory, based on Taylor-expanding $`V`$ about $`z_{jk}`$. Matching the hydrodynamic and core solutions smoothly together (see Refs. JRA ; RP ) finally yields, to leading order in $`\epsilon `$,

$`i\rho \dot{D}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\left[_{\overline{z}}\left(\rho _zD+\rho _{\overline{z}}\overline{D}\right)_z\left(\rho _{\overline{z}}D\right)\right]_{\overline{z}}\left(\rho F_c\right)`$ (13)
$`+\left(\mathrm{ln}{\displaystyle \frac{b}{2\pi \xi }}+1.17\right)_{\overline{z}}\rho 0.20_z{\displaystyle \frac{\left(_{\overline{z}}\rho \right)^2}{\rho }}`$
which is expressed in the co-rotating frame.

Eqns. (12) and (13) are our main results. The numerical co-efficients in (13) include some numerically evaluated contributions from the nonlinear core regions (compare with RP ), and also functions of $`\tau `$, generally related to $`\theta `$-functions, evaluated for the triangular case $`\tau =(1+i\sqrt{3})/2`$. For general $`\tau `$ (i.e. for lattices other than the regular triangular), (13) would have several additional terms, such as $`B(\tau )_{\overline{z}\overline{z}}D`$, where $`B(\tau )`$ is another rapidly converging series. Modular covariance (a general type of lattice symmetry) of the extra co-efficients like $`B(\tau )`$ constrains them to vanish when $`\tau =\frac{1+i\sqrt{3}}{2}`$.

MSA implies that the lattice scale $`b`$ is to the vortex matter equations much as the healing length $`\xi `$ is to the hydrodynamic equations that underly them. Thus Eqns. (12) and (13) should be accurate except at distances of order $`b`$ or less from of the edge of a vortex array. But is this claim really compatible with the results of Campbell and Ziff CampbellZiff  for finite vortex arrays in infinite homogeneous superfluid? Setting $`\rho `$ constant, to leading order in $`\epsilon ^2`$ the pair (12) (13) reduces to (2) and (3), and setting $`\dot{D}\dot{D}_0=0`$, we find a wealth of solutions to these fourth order equations:

$$D_0=\underset{m=1}{\overset{\mathrm{}}{}}\left(a_m\overline{z}^{m1}+(m+1)b_m\overline{z}^mzb_m^{}z^{m+1}\right)$$
(14)
for arbitrary complex constants $`a_m,b_m`$. A priori it is unclear what boundary condition $`D`$ should respect at the array edge; however if we assume that the edge should be a uniform circle of vortices, fitting leads to a unique combination of multipolar distortions with $`m=6,12,18\mathrm{}`$. Figure 2 shows that stopping at only $`m=12`$ gives quite good agreement with the first excited state for 217 vortices found in Ref. CampbellZiff . (The ground state differs only in the outermost ring.)

If we set $`\rho =\rho _0(1\epsilon ^2|z|^2)`$ for a BEC in a round harmonic trap, though, the only distortion forced on us by the inhomogeneity is the very mild

$$D_0=\frac{\sqrt{3}}{4\pi }\frac{\epsilon ^2z}{1\epsilon ^2|z|^2}\left(\mathrm{ln}\frac{b}{2\pi \xi }+1.17\right)+𝒪\left(\epsilon ^4\right)$$
(15)
This scarcely visible radial shift of each vortex is shown in Fig. 1, for $`\epsilon =1/6=\xi (0)/b`$. For vortices very close to the TF surface where formally $`\xi \mathrm{}`$, (15) spuriously predicts large inward displacements. (Four such vortices have been excised from Fig. 1a.) Apart from this failure in the lattice-edge boundary layer of thickness $`b`$, the accord with experiment is excellent.

Using (9) and (13), the $`v_c`$ associated with this $`D_0`$ is

$$v_{c0}=iz\left(\mathrm{\Omega }\frac{\epsilon ^2[\mathrm{ln}\frac{b}{2\pi \xi }+1.17]}{1\epsilon ^2|z|^2}\right).$$
(16)
This purely azimuthal flow is slightly less than rigid body rotation at $`\mathrm{\Omega }`$, and the magnitude of the backflow increases with radius. This agrees qualitatively with the numerical results shown in Fig. 5 of Ref. David .

Any static distortions forced by boundary conditions, like those of Fig. (2), would appear as zero-frequency modes among the collective excitations; so to these we now turn. We write $`D=D_0+_{\overline{z}}(P+iQ)`$, without loss of generality, for real $`P,Q`$. Then $`\mathrm{Re}\{_z[`$(13)$`/\rho ]\}`$ and (12) yield

$$_{\overline{z}z}\dot{Q}=2\mathrm{\Omega }_{\overline{z}z}P\times \left[1+𝒪(\epsilon ^2)\right].$$
(17)
Since only $`D`$ is physical, any terms in $`P`$ annihilated by $`_{\overline{z}z}`$ are of form $`f(z)+f^{}(\overline{z})`$ and so can be absorbed in $`Q`$ as $`i[f^{}(\overline{z})f(z)]`$. And since one can easily show that Laplacian-free terms in $`Q`$ must be time-independent to leading order in $`\epsilon ^2`$, we can set $`P=\frac{1}{2\mathrm{\Omega }}\dot{Q}`$.

Introducing polar co-ordinates $`re^{i\varphi }=\epsilon z`$, so that $`\rho =\rho _0(1r^2)`$, we can write $`Q=q_{mn}(r)\mathrm{cos}\left(m\varphi \omega _{mn}t\alpha \right)`$ (for arbitrary constant $`\alpha `$ and angular and radial quantum numbers $`m`$ and $`n`$). Considering first the case $`m=0,`$ where $`_\varphi F_c=0`$, the imaginary part of $`e^{i\varphi }\times `$(13) implies

$$\frac{4\omega _{0n}^2}{\epsilon ^2\mathrm{\Omega }}_rq_{0n}=\frac{\left(_r+2r^1\right)}{\left(1r^2\right)}\left[\left(1r^2\right)\left(_rr^1\right)_rq_{0n}\right]$$
(18)
plus order $`\epsilon .`$ Frobenius analysis shows that only one solution to this second order equation for $`_rq_{0n}`$ is finite at $`r0`$, and imposing finiteness at $`r1`$ as well fixes the discrete spectrum. (Modified behavior within the boundary layer $`r1\epsilon `$ will be able to reconcile our coarse-scale $`D`$ and $`v_c`$ with microscopic boundary conditions, as long as our functions remain regular; compare the hydrodynamic derivation of collective modes in the vortex-free BEC Sandro1 .) Numerical search yields

$$\omega _{0n}=\epsilon \mathrm{\Omega }\{0,1.43,2.32,3.18,4.03,\mathrm{}\}.$$
(19)
The zero mode is global rotation of the lattice, $`q_{00}=r^2`$.

These results are unsurprising for Tkachenko waves in a finite cylindrical system; but for $`m0`$, the dynamics is very different. To eliminate $`F_c`$ we must take $`\mathrm{Im}[_z`$(13)$`]`$, and the non-constant $`\rho `$ leaves a first order time derivative on the left side, implying $`\omega _{mn}`$ of order $`\epsilon ^2`$ instead of $`\epsilon .`$ Our leading order equation is thus first order in $`t`$ but fourth in $`r`$:

$`{\displaystyle \frac{16m\omega _{mn}}{\epsilon ^2}}q_{mn}`$ $`=`$ $`\left[\left(_{rr}+r^1_rm^2r^2\right)4r_r\right]`$ (20)
$`\times \left(_{rr}+r^1_rm^2r^2\right)q_{mn}.`$
This equation is quite singular (Frobenius analysis shows that of the four solutions only one is finite at both $`r=0,1`$); but its differential operator is Hermitian and self-dual, and its regular solution has a rapidly converging Frobenius series in $`r.`$ From (13) we see that the radial component of $`v_c`$ blows up at $`r=1`$ unless $`\left(_{rr}r^1_rm^2r^2\right)q_{mn}(1)=0`$. This boundary condition on the regular solutions $`q_{mn}(r)`$ fixes the discrete spectrum, which may be found by summing Frobenius series numerically:

| $`\frac{\omega _{mn}}{\epsilon ^2\mathrm{\Omega }}`$ | $`m=1`$ | 2 | 3 | 4 | 5 | 6 |
| --- | --- | --- | --- | --- | --- | --- |
| $`n=0`$ | 0 | 0.365 | 0.900 | 1.60 | 2.46 | 3.49 |
| 1 | 2.93 | 5.02 | 7.53 | 10.5 | 13.8 | 17.5 |
| 2 | 32.0 | 31.3 | 35.6 | 41.6 | 48.8 | 56.8 |
| 3 | 130. | 105. | 106. | 113. | 124. | 137. |

and so on. The translational zero mode is $`q_{10}=r`$. These eigenvalues of a fourth order equation increase rapidly with radial quantum number $`n`$, and as the co-efficients reach order $`\epsilon ^1`$ (20) becomes invalid, and WKB-like Tkachenko waves will emerge instead.

The only zero-frequency solution satisfying the boundary conditions is the rigid translation; and so (15) is the full static distortion. This reflects the fact that $`\rho `$ is so small at the edges of the lattice that no boundary energies are large enough to influence the bulk lattice.

Only positive $`m`$ have been reported, because $`\omega _{m,n}=\omega _{mn}`$ and replacing $`mm`$ leaves our ansatz for $`Q`$ unchanged. What this means is that the nonuniform $`\rho `$ has drastically split the degeneracy of the two modes that would, for constant $`\rho `$, be proportional to $`e^{\pm im\varphi }`$, with frequency of order $`\epsilon `$. The linear combination proportional to $`\mathrm{cos}(m\varphi \omega _{mn}t)`$ propagates much more slowly, in the co-rotating frame; evidently the orthogonal combination, which we have not examined, propagates much more quickly. (The distortion patterns rotate about the origin; the vortices follow elliptical orbits about their equilibrium positions: see Fig. 3.) Thus the lowest frequency lattice modes have a first order dynamics, and only half as many distinct modes as for constant $`\rho `$.

Finally, note that for the quadrupole mode $`\omega _{20}0`$ (the zero eigenvalue solution $`\stackrel{~}{q}_{20}=r^2`$ does not satisfy the boundary condition). Since it is this mode which would distort the equilateral triangular lattice into the moderately different regular lattices that are also dynamically stable on short wavelengths T1 , we conclude that although stable in bulk those lattices are frustrated in the finite system, and cannot even be stationary without slow but cumulatively large distortions. Reviewing our calculations in this context, it is clear that the only reason the equilateral triangular lattice does not suffer a similar fate is the vanishing, due to lattice symmetry, of several awkward terms from (13). Once the threat of cumulative distortion is lifted, it is not surprising that merely local distortion is of order $`\epsilon ^2`$. So the regularity of the observed BEC vortex lattices is ultimately due to their triangular structure. An engineer would attribute this to triangular rigidity, and a mathematician to the fact that the triangular lattice is the $`Z_3`$ fixed point of the modular group. Physicists are entitled to arbitrary linear combinations of the two explanations.

This work has been supported by the US National Science Foundation through a grant for the Center for Ultracold Atoms and through a grant for the Institute for Theoretical Atomic and Molecular Physics at Harvard University and the Smithsonian Astrophysical Observatory. MC acknowledges partial support from the DOE EPSCoR grant DE-FG01-00ER45832, Research Corporation Cottrell Science Award #CC5285, and a YSU Research Professorship Award.