text:
# A New Proof of the Integral Localization Formula for Equivariantly Closed Differential Forms

## 1 Introduction

Equivariant forms were introduced in 1950 by Henri Cartan. There are many good texts on this subject including \[BGV\] and \[GS\].

Let $`G`$ be a compact Lie group acting on a compact oriented manifold $`M`$, let $`𝔤`$ be the Lie algebra of $`G`$, and let $`\alpha :𝔤𝒜(M)`$ be an equivariantly closed form on $`M`$. For $`X𝔤`$, we denote by $`M_0(X)`$ the set of zeroes of the vector field on $`M`$ induced by the infinitesimal action of $`X`$. We assume that $`M_0(X)`$ is discrete. Then Theorem 7.11 in \[BGV\] says that the integral of $`\alpha (X)`$ can be expressed as a sum over the set of zeroes $`M_0(X)`$ of certain local quantities of $`M`$ and $`\alpha `$:

$$_M\alpha (X)=\underset{pM_0(X)}{}\text{local invariant of }M\text{ and }\alpha \text{ at }p.$$

This is the essence of the integral localization formula for equivariantly closed differential forms (Theorem 2).

This localization formula for compact Lie groups $`G`$ is very well known, but the problem with previously existed proofs is that they rely heavily on compactness of the acting group $`G`$ and none of them generalizes to a noncompact group setting. In the future, the author intends to adapt this proof and obtain a localization formula which can be applied to some nontrivial integrals that cannot be reduced to the compact group case.

The proof which appears here is a significant modification of the localization argument which appeared in my Ph.D. thesis \[L1\]. This thesis contains a result on localizations of actions of noncompact groups, which is probably the only existing one on this subject. Article \[L2\] gives a very accessible introduction to \[L1\] and explains key ideas used there by way of examples.

I would like to thank Sergey Vasilyev and Michele Vergne for pointing out some very serious errors in an earlier version of this article.

## 2 Equivariant Forms

In this article we use the same notations as in \[BGV\].

Let $`M`$ be a $`𝒞^{\mathrm{}}`$-manifold of dimension $`n`$ with an action of a (possibly noncompact) Lie group $`G`$, and let $`𝔤`$ be the Lie algebra of $`G`$. The group $`G`$ acts on $`𝒞^{\mathrm{}}(M)`$ by the formula $`(g\phi )(x)=\phi (g^1x)`$. For $`X𝔤`$, we denote by $`X_M`$ the vector field on $`M`$ given by (notice the minus sign)

$$(X_M\phi )(x)=\frac{d}{d\epsilon }\phi \left(\mathrm{exp}(\epsilon X)x\right)|_{\epsilon =0}.$$

Let $`𝒜(M)`$ denote the (graded) algebra of smooth differential forms on $`M`$, and let $`𝒞^{\mathrm{}}(𝔤)\widehat{}𝒜(M)`$ denote the algebra of all smooth $`𝒜(M)`$-valued functions on $`𝔤`$. The group $`G`$ acts on an element $`\alpha 𝒞^{\mathrm{}}(𝔤)\widehat{}𝒜(M)`$ by the formula

$$(g\alpha )(X)=g(\alpha (g^1X))\text{for all }gG\text{ and }X𝔤\text{.}$$

Let $`𝒜_G^{\mathrm{}}(M)=(𝒞^{\mathrm{}}(𝔤)\widehat{}𝒜(M))^G`$ be the subalgebra of $`G`$-invariant elements. An element $`\alpha `$ of $`𝒜_G^{\mathrm{}}(M)`$ satisfies the relation $`\alpha (gX)=g\alpha (X)`$ and is called and equivariant differential form.

We define the equivariant exterior differential $`d_𝔤`$ on $`𝒞^{\mathrm{}}(𝔤)\widehat{}𝒜(M)`$ by the formula

$$(d_𝔤\alpha )(X)=d(\alpha (X))\iota (X_M)(\alpha (X)),$$

where $`\iota (X_M)`$ denotes contraction by the vector field $`X_M`$. This differential $`d_𝔤`$ preserves $`𝒜_G^{\mathrm{}}(M)`$, and $`(d_𝔤)^2\alpha =0`$ for all $`\alpha 𝒜_G^{\mathrm{}}(M)`$. The elements of $`𝒜_G^{\mathrm{}}(M)`$ such that $`d_𝔤\alpha =0`$ are called equivariantly closed forms.

###### Example 1

Let $`T^{}M`$ be the cotangent bundle of $`M`$, and let $`\sigma `$ denote the canonical symplectic form on $`T^{}M`$. It is defined, for example, in \[KaScha\], Appendix A2. The action of the Lie group $`G`$ on $`M`$ naturally extends to $`T^{}M`$. Then we always have a canonical equivariantly closed form on $`T^{}M`$, namely, $`\mu +\sigma `$. Here $`\mu :𝔤𝒞^{\mathrm{}}(T^{}M)`$ is the moment map defined by:

$$\mu (X):\xi \xi ,X_M,X𝔤,\xi T^{}M.\mathrm{}$$

If $`\alpha `$ is a non-homogeneous equivariant differential form, $`\alpha _{[k]}`$ denotes the homogeneous component of degree $`k`$. If $`M`$ is a compact oriented manifold, we can integrate equivariant differential forms over $`M`$, obtaining a map

$$_M:𝒜_G^{\mathrm{}}(M)𝒞^{\mathrm{}}(𝔤)^G,$$

defined by the formula $`(_M\alpha )(X)=_M\alpha (X)_{[n]}`$, where $`n=dimM`$.

Finally, notice that if $`\alpha `$ is an equivariantly closed form whose top homogeneous component has degree $`k`$, then $`\alpha (X)_{[k]}`$ is closed with respect to the ordinary exterior differential.

## 3 Localization Formula

From now on we assume that the Lie group $`G`$ is compact.

We recall some more notations from \[BGV\]. Let $`M_0(X)`$ be the set of zeroes of the vector field $`X_M`$. We will state and prove the localization formula in the important special case where $`X_M`$ has isolated zeroes. Here, at each point $`pM_0(X)`$, the infinitesimal action of $`X`$ on $`𝒞^{\mathrm{}}(M)`$ gives rise to a linear transformation $`L_p`$ on $`T_pM`$.

Because the Lie group $`G`$ is compact, the transformation $`L_p`$ is invertible and has only imaginary eigenvalues. Thus the dimension of $`M`$ is even and there exists an oriented basis $`\{e_1,\mathrm{},e_n\}`$ of $`T_pM`$ such that for $`1il=n/2`$,

$$L_pe_{2i1}=\lambda _{p,i}e_{2i},L_pe_{2i}=\lambda _{p,i}e_{2i1}.$$

We have $`det(L_p)=\lambda _{p,1}^2\lambda _{p,2}^2\mathrm{}\lambda _{p,l}^2`$, and it is natural to take the following square root (dependent only on the orientation of the manifold):

$$\stackrel{1/2}{det}(L_p)=\lambda _{p,1}\mathrm{}\lambda _{p,l}.$$

For convenience, we restate Theorem 7.11 from \[BGV\].

###### Theorem 2

Let $`G`$ be a compact Lie group with Lie algebra $`𝔤`$ acting on a compact oriented manifold $`M`$, and let $`\alpha `$ be an equivariantly closed differential form on $`M`$. Let $`X𝔤`$ be such that the vector field $`X_M`$ has only isolated zeroes. Then

$$_M\alpha (X)=(2\pi )^l\underset{pM_0(X)}{}\frac{\alpha (X)(p)}{det^{1/2}(L_p)},$$

where $`l=dim(M)/2`$, and by $`\alpha (X)(p)`$ we mean the value of the function $`\alpha (X)_{[0]}`$ at the point $`pM`$.

## 4 Proof of the Localization Formula

Proof of Theorem 2. Let $`\pi :T^{}MM`$ denote the projection map. We regard $`M`$ as a submanifold of $`T^{}M`$ via the zero section inclusion. Recall the canonical equivariantly closed form $`\mu +\sigma `$ on $`T^{}M`$ constructed in Example 1. We can consider the form

$$\omega (X)=e^{\mu (X)+\sigma }\pi ^{}\left(\alpha (X)\right),X𝔤.$$
(1)
It is an equivariantly closed form on $`T^{}M`$ for the reason that it is “assembled” from equivariantly closed forms. Moreover, its restriction to $`M`$ is just $`\alpha (X)`$. We will see later that, in a way, $`\omega `$ is the most natural equivariant extension of $`\alpha `$ to $`T^{}M`$.

Fix a $`G`$-invariant Riemannian metric on $`T^{}M`$ and let $`.`$ be the norm. It induces a vector bundle diffeomorphism $`TMT^{}M`$ and a norm on $`TM`$ which we also denote by $`.`$. Let $`s:MT^{}M`$ denote the image of the vector field $`X_M`$ under this diffeomorphism. Observe that

$$\mu (X)(s(m))=X_M^20,\text{for all }mM,$$

and that the set

$$\{mM;\mu (X)(s(m))=0\}=\{mM;X_M=0\}=M_0(X)$$

is precisely the set of zeroes of $`X_M`$.

Pick an $`R>0`$ and let $`t`$ increase from 0 to $`R`$. Then the family of sections $`\{ts\}_{t[0,R]}`$ provides a deformation of the initial cycle $`M`$ into $`Rs`$. In other words, there exists a chain $`C_R`$ in $`T^{}M`$ of dimension $`(dimM+1)`$ such that

$$C_R=MRs.$$
(2)
Note that the support of $`C_R`$, $`|C_R|`$, lies inside the set

$$\{\xi T^{}M;\mu (X)(\xi )=\xi ,X_M0\}.$$
(3)
This will ensure good behavior of the integrand when we let $`R`$ tend to infinity.

Observe that if we let $`R\mathrm{}`$, then we obtain a conic Borel-Moore chain $`C_{\mathrm{}}`$ such that

$$C_{\mathrm{}}=M\underset{pM_0(X)}{}T_p^{}M,$$

where each cotangent space $`T_p^{}M`$ is provided with appropriate orientation. But we will not use this observation in the proof because one cannot interchange the order of integration and taking limit as $`R\mathrm{}`$.

###### Lemma 3

$`\left(d\omega (X)_{[n]}\right)|_{|C_R|}=0`$.

###### Remark 4

The form $`\omega (X)_{[n]}`$ itself need not be closed; this lemma only says that $`\omega (X)_{[n]}`$ becomes closed when restricted to $`|C_R|`$.

This lemma is the key moment in the proof and it is the only place where we use that the form $`\alpha `$ is equivariantly closed and that our Riemannian metric is $`G`$-invariant. Its proof is quite technical.

Proof. The form $`\omega (X)_{[n]}`$ can be written more explicitly as

$$\begin{array}{c}\omega (X)_{[n]}=\left(e^{\mu (X)+\sigma }\alpha (X)\right)_{[n]}\hfill \\ \hfill =e^{\mu (X)}\left((e^\sigma )_{[0]}\alpha (X)_{[n]}+(e^\sigma )_{[2]}\alpha (X)_{[n2]}+(e^\sigma )_{[4]}\alpha (X)_{[n4]}+\mathrm{}\right)\\ \hfill =e^{\mu (X)}\left(\frac{1}{0!}\sigma ^0\alpha (X)_{[n]}+\frac{1}{1!}\sigma ^1\alpha (X)_{[n2]}+\frac{1}{2!}\sigma ^2\alpha (X)_{[n4]}+\mathrm{}\right)\\ \hfill =e^{\mu (X)}\underset{k=0}{\overset{l}{}}\frac{1}{k!}\sigma ^k\alpha (X)_{[n2k]}.\end{array}$$

It is sufficient to prove that the form $`\omega (X)_{[n]}|_{|C_R|}`$ is closed away from zeroes of $`X_M`$, i.e. on the set $`|C_R|M_0(X)`$. So let us pick $`\xi |C_R|M_0(X)`$, say $`\xi T_q^{}M`$, $`qM_0(X)`$. Because $`X_M(q)0`$ we can trivialize vector field $`X_M`$ near $`q`$. That is, there is a neighborhood $`U_q`$ of $`M`$ containing point $`q`$ and positively oriented local coordinates $`(x_{q,1},\mathrm{},x_{q,n})`$ defined on $`U_q`$ such that $`X_M`$ is expressed by

$$X_M=\frac{}{x_{q,1}}.$$
(4)
Define a function $`r:|C_R|`$ by $`r(c)=c,\mu (X)`$, $`c|C_R|T^{}M`$. Then $`(x_{q,1},\mathrm{},x_{q,n},r)`$ defined on $`\pi ^1(U_q)|C_R|`$ is a system of coordinates on $`|C_R|`$ containing $`\xi `$.

We expand the coordinate system $`(x_{q,1},\mathrm{},x_{q,n})`$ of $`M`$ defined on $`U_q`$ to the standard coordinate system $`(x_{q,1},\mathrm{},x_{q,n},\zeta _{q,1},\mathrm{},\zeta _{q,n})`$ of $`T^{}M`$ defined on $`T^{}U_q`$ so that every element of $`T^{}U_q`$ is expressed in these coordinates as $`(x_{q,1},\mathrm{},x_{q,n},\zeta _{q,1}dx_{q,1}+\mathrm{}+\zeta _{q,n}dx_{q,n})`$. In these coordinates,

$$\sigma =dx_{q,1}d\zeta _{q,1}+\mathrm{}+dx_{q,n}d\zeta _{q,n}.$$

Using this coordinate system $`(x_{q,1},\mathrm{},x_{q,n},\zeta _{q,1},\mathrm{},\zeta _{q,n})`$ we can regard $`X_M`$ as a horizontal vector field on $`T^{}U_q`$. Because we chose our Riemannian metric to be $`G`$-invariant, the vector field $`X_M`$ is tangent to $`|C_R|`$.

Using our assumption that the form $`\alpha `$ is equivariantly closed we obtain:

$$\begin{array}{c}d\left(\omega (X)_{[n]}|_{|C_R|}\right)=d\left(e^r\underset{k=0}{\overset{l}{}}\frac{1}{k!}\sigma ^k\alpha (X)_{[n2k]}\right)|_{|C_R|}\hfill \\ \hfill =dr\left(e^r\underset{k=0}{\overset{l}{}}\frac{1}{k!}\sigma ^k\alpha (X)_{[n2k]}\right)|_{|C_R|}+e^rd\left(\underset{k=0}{\overset{l}{}}\frac{1}{k!}\sigma ^k\alpha (X)_{[n2k]}\right)|_{|C_R|}\\ \hfill =e^r\left(\underset{k=0}{\overset{l}{}}\frac{1}{k!}dr\sigma ^k\alpha (X)_{[n2k]}+\underset{k=0}{\overset{l}{}}\frac{1}{k!}\sigma ^kd\alpha (X)_{[n2k]}\right)|_{|C_R|}\\ \hfill =e^r\left(\underset{k=0}{\overset{l}{}}\frac{1}{k!}dr\sigma ^k\alpha (X)_{[n2k]}+\underset{k=1}{\overset{l}{}}\frac{1}{k!}\sigma ^k\iota (X_M)\left(\alpha (X)_{[n2k+2]}\right)\right)|_{|C_R|}\\ \hfill =e^r\left(\underset{k=0}{\overset{l}{}}\frac{1}{k!}dr\sigma ^k\alpha (X)_{[n2k]}\underset{k=1}{\overset{l+1}{}}\frac{1}{k!}\iota (X_M)\left(\sigma ^k\right)\alpha (X)_{[n2k+2]}\right)|_{|C_R|}\end{array}$$

because $`\sigma ^{l+1}=0`$ and $`\left(\sigma ^k\alpha (X)_{[n2k+2]}\right)|_{|C_R|}=0`$ for dimension reasons, and so

$$\begin{array}{c}0=\iota (X_M)\left(\sigma ^k\alpha (X)_{[n2k+2]}\right)|_{|C_R|}\hfill \\ \hfill =\iota (X_M)\left(\sigma ^k\right)\alpha (X)_{[n2k+2]}|_{|C_R|}+\sigma ^k\iota (X_M)\left(\alpha (X)_{[n2k+2]}\right)|_{|C_R|}.\end{array}$$

Hence:

$$\begin{array}{c}d\left(\omega (X)_{[n]}|_{|C_R|}\right)=e^r(\underset{k=0}{\overset{l}{}}\frac{1}{k!}dr\sigma ^k\alpha (X)_{[n2k]}\hfill \\ \hfill \underset{k=1}{\overset{l+1}{}}\frac{1}{(k1)!}\iota (X_M)\left(\sigma \right)\sigma ^{k1}\alpha (X)_{[n2k+2]})|_{|C_R|}\\ \hfill =e^r\left(\underset{k=0}{\overset{l}{}}\frac{1}{k!}dr\sigma ^k\alpha (X)_{[n2k]}\underset{k=0}{\overset{l}{}}\frac{1}{k!}d\zeta _{q,1}\sigma ^k\alpha (X)_{[n2k]}\right)|_{|C_R|}\\ \hfill =e^r\left(\underset{k=0}{\overset{l}{}}\frac{1}{k!}(drd\zeta _{q,1})\sigma ^k\alpha (X)_{[n2k]}\right)|_{|C_R|}=0\end{array}$$

because equation (4) implies that in these coordinates, for $`c|C_R|`$,

$$\zeta _{q,1}(c)=c,X_M=c,\mu (X)=r(c),$$

hence $`dr=d\zeta _{q,1}`$. This finishes the proof of Lemma 3. $`\mathrm{}`$

Therefore,

$$_M\alpha (X)=_M\omega (X)=_{Rs}\omega (X).$$
(5)
For each $`pM_0(X)`$, there exists a positively oriented orthonormal basis $`\{e_1,\mathrm{},e_n\}`$ of $`T_pM`$ such that for $`1il=n/2`$,

$$L_pe_{2i1}=\lambda _{p,i}e_{2i},L_pe_{2i}=\lambda _{p,i}e_{2i1}.$$

Using the exponential map we can find a neighborhood $`U_p`$ of $`M`$ containing point $`p`$ and local coordinates $`(x_{p,1},\mathrm{},x_{p,n})`$ defined on $`U_p`$ and centered at $`p`$ such that

$$\frac{}{x_{p,1}}(p)=e_1,\mathrm{},\frac{}{x_{p,n}}(p)=e_n$$

(in particular, the coordinate system $`(x_{p,1},\mathrm{},x_{p,n})`$ is positively oriented) and the vector field $`X_M`$ is expressed by

$$\begin{array}{c}X_M=\lambda _{p,1}\left(x_{p,2}\frac{}{x_{p,1}}x_{p,1}\frac{}{x_{p,2}}\right)+\mathrm{}+\lambda _{p,l}\left(x_{p,2l}\frac{}{x_{p,2l1}}x_{p,2l1}\frac{}{x_{p,2l}}\right).\hfill \end{array}$$

Since $`M_0(X)`$ is discrete, we can make these neighborhoods $`U_p`$’s small enough so that they do not overlap.

Find a sufficiently small $`\epsilon >0`$ such that, for each $`pM_0(X)`$, the cube

$$\{(t_1,\mathrm{},t_n);|t_1|\epsilon ,\mathrm{},|t_n|\epsilon \}$$

is contained in the image of $`(x_{p,1},\mathrm{},x_{p,n})`$. For each $`pM_0(X)`$, we define a subset of $`M`$ containing $`p`$

$$I_p=(x_{p,1},\mathrm{},x_{p,n})^1\left(\{(t_1,\mathrm{},t_n);|t_1|\epsilon ,\mathrm{},|t_n|\epsilon \}\right).$$

In formula (5) we let $`R`$ tend to infinity. Recall the definition of the equivariant form $`\omega `$ (formula (1)). Because of the presence of the term $`e^{\mu (X)+\sigma }`$ there and because of the exponential decay on the set (3) we get:

$$\begin{array}{c}_M\alpha (X)=\underset{R\mathrm{}}{lim}_{Rs}\omega (X)=\underset{R\mathrm{}}{lim}_{Rs_{pM_0(X)}\pi ^1I_p}\omega (X)\hfill \\ \hfill =\underset{pM_0(X)}{}\underset{R\mathrm{}}{lim}_{Rs\pi ^1I_p}\omega (X).\end{array}$$

That is, when we integrate over $`Rs`$, only the covectors whose basepoint lies in some $`I_p`$ (i.e. “near” $`M_0(X)`$) count. Thus we see that the integral is localized around the set of zeroes $`M_0(X)`$. It remains to calculate the contribution of each point $`p`$:

$$\underset{R\mathrm{}}{lim}_{Rs\pi ^1I_p}\omega (X).$$

Recall that

$$\left\{\frac{}{x_{p,1}}=e_1,\mathrm{},\frac{}{x_{p,n}}=e_n\right\}$$

is an orthonormal basis of $`T_p^{}M`$. Hence the Riemannian metric around each $`p`$ can be written in terms of the frame on $`U_p`$

$$(\frac{}{x_{p,1}},\mathrm{},\frac{}{x_{p,n}})$$

as $`I+E_p(m)`$, where $`I`$ is the $`n\times n`$ identity matrix and $`E_p(m)`$ is a symmetric $`n\times n`$ matrix which depends on $`m`$ and vanishes at the point $`p`$.

We expand each coordinate system $`(x_{p,1},\mathrm{},x_{p,n})`$ of $`M`$ defined on $`U_p`$ to the standard coordinate system $`(x_{p,1},\mathrm{},x_{p,n},\zeta _{p,1},\mathrm{},\zeta _{p,n})`$ of $`T^{}M`$ defined on $`T^{}U_p`$ so that every element of $`T^{}U_p`$ is expressed in these coordinates as $`(x_{p,1},\mathrm{},x_{p,n},\zeta _{p,1}dx_{p,1}+\mathrm{}+\zeta _{p,n}dx_{p,n})`$.

In these coordinates,

$$\begin{array}{c}s=(x_{p,1},\mathrm{},x_{p,n},\lambda _{p,1}(x_{p,2}dx_{p,1}+x_{p,1}dx_{p,2})+\mathrm{}+\lambda _{p,l}(x_{p,2l}dx_{p,2l1}+x_{p,2l1}dx_{p,2l})+O(x^2)),\hfill \end{array}$$

where $`O(x^2)`$ denotes terms of quadratic order in $`x_{p,1},\mathrm{},x_{p,n}`$. We also have:

$$\sigma =dx_{p,1}d\zeta _{p,1}+\mathrm{}+dx_{p,n}d\zeta _{p,n},$$

and

$$\begin{array}{c}\mu (X):(x_{p,1},\mathrm{},x_{p,n},\zeta _{p,1},\mathrm{},\zeta _{p,n})\hfill \\ \hfill \lambda _{p,1}(x_{p,2}\zeta _{p,1}x_{p,1}\zeta _{p,2})+\mathrm{}+\lambda _{p,l}(x_{p,2l}\zeta _{p,2l1}x_{p,2l1}\zeta _{p,2l}).\end{array}$$

Let us choose a coordinate system for $`Rs\pi ^1I_p`$, namely

$$(y_{p,1}=x_{p_1}\pi ,\mathrm{},y_{p,n}=x_{p,n}\pi ).$$

It is positively oriented. Then in these coordinates

$$dx_{p,1}=dy_{p,1},\mathrm{},dx_{p,n}=dy_{p,n},$$

$$\begin{array}{c}d\zeta _{p,1}=R\lambda _{p,1}dy_{p,2}+RO(y),d\zeta _{p,2}=R\lambda _{p,1}dy_{p,1}+RO(y),\mathrm{},\hfill \\ \hfill d\zeta _{p,2l1}=R\lambda _{p,l}dy_{p,2l}+RO(y),d\zeta _{p,2l}=R\lambda _{p,l}dy_{p,2l1}+RO(y),\end{array}$$

$$\mu (X)|_{Rs}=R\lambda _{p,1}^2(y_{p,1}^2+y_{p,2}^2)\mathrm{}R\lambda _{p,l}^2(y_{p,2l1}^2+y_{p,2l}^2)+RO(y^3).$$

Thus $`\sigma |_{Rs}`$ in these coordinates becomes

$$\sigma |_{Rs}=2R\lambda _{p,1}dy_{p,1}dy_{p,2}\mathrm{}2R\lambda _{p,l}dy_{p,2l1}dy_{p,2l}+RO(y).$$

Hence

$$(e^\sigma )_{[n]}|_{Rs}=(2R)^l\lambda _{p,1}\mathrm{}\lambda _{p,l}dy_{p,1}\mathrm{}dy_{p,n}+R^lO(y).$$

On the other hand, the expression of $`\pi ^{}\alpha (X)|_{Rs}`$ in these coordinates does not depend on $`R`$. We have:

$$\begin{array}{c}\omega (X)_{[n]}\hfill \\ \hfill =\left(e^{\mu (X)+\sigma }\right)_{[n]}\alpha (X)_{[0]}+\left(e^{\mu (X)+\sigma }\right)_{[n2]}\alpha (X)_{[2]}+\mathrm{}+\left(e^{\mu (X)+\sigma }\right)_{[0]}\alpha (X)_{[n]}.\end{array}$$

Hence

$$\begin{array}{c}\omega (X)_{[n]}|_{Rs}=\left(e^{\mu (X)+\sigma }\right)_{[n]}\alpha (X)_{[0]}|_{Rs}\hfill \\ \hfill +e^{\mu (X)}\left(\text{terms depending on }R\text{ through powers of }R\text{ strictly less than }\frac{n}{2}\right)\\ \hfill =e^{\mu (X)}(e^\sigma )_{[n]}\alpha (X)_{[0]}|_{Rs}+e^{\mu (X)}o(R^{n/2})=\\ \hfill e^{R\left(\lambda _{p,1}^2(y_{p,1}^2+y_{p,2}^2)+\mathrm{}+\lambda _{p,l}^2(y_{p,2l1}^2+y_{p,2l}^2)+O(y^3)\right)}\left((e^\sigma )_{[n]}\alpha (X)_{[0]}|_{Rs}+o(R^{n/2})\right).\end{array}$$

Recall that the integration takes place over $`I_p`$, hence the range of variables is $`|y_{p,1}|\epsilon ,\mathrm{},|y_{p,n}|\epsilon `$. We perform a change of coordinates

$$z_{p,1}=\sqrt{R}y_{p,1},\mathrm{},z_{p,n}=\sqrt{R}y_{p,n},|z_{p,1}|\epsilon \sqrt{R},\mathrm{},|z_{p,n}|\epsilon \sqrt{R}.$$

In these new coordinates,

$$\begin{array}{c}\omega (X)_{[n]}|_{Rs}=e^{\left(\lambda _{p,1}^2(z_{p,1}^2+z_{p,2}^2)+\mathrm{}+\lambda _{p,l}^2(z_{p,2l1}^2+z_{p,2l}^2)+O(z^3)/O(\sqrt{R})\right)}\hfill \\ \hfill \left((e^\sigma )_{[n]}\alpha (X)_{[0]}|_{Rs}+1/O(\sqrt{R})\right),\end{array}$$

and

$$(e^\sigma )_{[n]}|_{Rs}=(2)^l\lambda _{p,1}\mathrm{}\lambda _{p,l}dz_{p,1}\mathrm{}dz_{p,n}+O(z)/O(\sqrt{R}).$$

As $`R`$ tends to infinity, the integrand tends pointwise to

$$\begin{array}{c}(2)^l\lambda _{p,1}\mathrm{}\lambda _{p,l}e^{\lambda _{p,1}^2(z_{p,1}^2+z_{p,2}^2)\mathrm{}\lambda _{p,l}^2(z_{p,n1}^2+z_{p,n}^2)}\hfill \\ \hfill \alpha (X)_{[0]}(p)dz_{p,1}\mathrm{}dz_{p,n}.\end{array}$$

Hence by the Lebesgue dominated convergence theorem,

$$\begin{array}{c}_M\alpha (X)=\underset{pM_0(X)}{}\underset{R\mathrm{}}{lim}_{Rs\pi ^1I_p}\omega (X)\hfill \\ \hfill =\underset{pM_0(X)}{}_^n(2)^l\lambda _{p,1}\mathrm{}\lambda _{p,l}\alpha (X)_{[0]}(p)\\ \hfill e^{\lambda _{p,1}^2(z_{p,1}^2+z_{p,2}^2)\mathrm{}\lambda _{p,l}^2(z_{p,n1}^2+z_{p,n}^2)}dz_{p,1}\mathrm{}dz_{p,n}\\ \hfill =\underset{pM_0(X)}{}(2\pi )^l\frac{1}{\lambda _{p,1}\mathrm{}\lambda _{p,l}}\alpha (X)_{[0]}(p)=(2\pi )^l\underset{pM_0(X)}{}\frac{\alpha (X)(p)}{det^{1/2}(L_p)}.\end{array}$$

This finishes our proof of Theorem 2. $`\mathrm{}`$

matvei@math.umass.edu

DEPARTMENT OF MATHEMATICS AND STATISTICS, UNIVERSITY OF MASSACHUSETTS, LEDERLE GRADUATE RESEARCH TOWER, 710 NORTH PLEASANT STREET, AMHERST, MA 01003


format: text

text:
# Bose-Einstein condensation dynamics in three dimensions by the pseudospectral and finite-difference methods

## 1 Introduction

The experimental realization  of Bose-Einstein condensates (BECs) in dilute weakly-interacting trapped bosonic atoms at ultra-low temperature initiated intense theoretical effort to describe the properties of the condensate . The properties of a condensate at zero temperature are usually described by the time-dependent, nonlinear, mean-field Gross-Pitaevskii (GP) equation . The effect of the interatomic interaction leads to a nonlinear term in the GP equation which complicates the solution procedure. Also, to simulate the proper experimental situation one should be prepared to deal with an anisotropic trap .

A numerical study of the time-dependent GP equation is of interest, as this can provide solution to many stationary and time-evolution problems. The time-independent GP equation yields only the solution of stationary problems. As our principal interest is in time evolution problems, we shall only consider the time-dependent GP equation in this paper. There are many numerical methods for the solution of the GP equation .

Here we suggest a pseudospectral time-iteration method  for the solution of the three-dimensional GP equation with an anisotropic harmonic trap and contrast it with the finite-difference method . In the pseudospectral method the unknown wave function is expanded in terms of a set of $`N`$ orthogonal polynomials. When this expansion is substituted into the GP equation, the (space) differential operators operate on a set of known polynomials and generate a differentiation matrix operating on the unknown coefficients. Consequently, the time-dependent partial-differential nonlinear GP equation in space and time variables is reduced to a set of $`N`$ coupled ordinary differential equations (ODEs) in time which is solved by a fourth-order adaptive step-size control Runge-Kutta method  using successive time iteration. In the pseudospectral method we use the Hermite polynomials to expand the wave function. In references  pseudospectral methods have been employed for the solution of the GP equation, where a variable step forth order Runge-Kutta time propagator was used, as in the present work. In  a pseudospectral fourier-sine basis was used for finite traps, and a corresponding complex pseudospectral basis was used for systems with periodic boundary conditions. However, the present study seems to be the first systematic study in dealing with the GP equation in three space dimensions using the pseudospectral and finite-difference approaches.

In the finite-difference method the time iteration is implemented by the split-step Crank-Nicholson scheme. The first approach will be termed pseudospectral-Runge-Kutta (PSRK) method and the second finite-difference-Crank-Nicholson (FDCN) method in the following. Here pseudospectral and finite-difference refer to the space part and Runge-Kutta and Crank-Nicholson refer to the time part. However, it should be noted that the temporal and the spatial parts in the GP equation can be dealt with in independent fashions. One can combine the pseudospectral or finite difference approach with a specific algorithm for time stepping which can be Runge-Kutta or Crank-Nicholson (with or without time-splitting) scheme, or any other. In fact it is quite common to implement some version of the pseudospectral method with Crank-Nicholson time stepping when solving the Navier-Stokes equation.

The PSRK method is illustrated by calculating the wave function and energy eigenvalue of the GP equation for different nonlinearity and trap symmetries in three dimensions as well as by studying a resonance dynamics. We find that in both cases the PSRK method presented here turns out to be a practical and efficient one for the solution of the time-dependent GP equation.

Resonance is an interesting feature of an oscillation under the action of an external periodic force manifesting in a large amplitude, when the frequency of the external force equals a multiple of the natural frequency of oscillation. Although, the phenomenon of resonance is well-understood in linear problems, in nonlinear dynamics it is far more nontrivial. Hence, it is worthwhile to study the dynamics of resonance using the nonlinear GP equation. The possibility of generating a resonance in a BEC subject to an oscillating trapping potential has been explored previously . Here, using the PSRK approach we study the resonance dynamics of a three-dimensional BEC subject to a periodic variation in the scattering length. One can have resonant oscillation in the $`x`$, $`y`$, and $`z`$ directions, when the frequency of oscillation of scattering length equals an even multiple of the trapping frequency in the respective direction. Such a variation in the scattering length is possible near a Feshbach resonance by manipulating an external magnetic field . There have been recent studies on this topic in one  and two space  variables.

Next we consider the FDCN method for the numerical solution of the time-dependent GP equation . In this approach the time-dependent GP equation is first discretized in space and time using a specific rule  and the resultant set of equations is solved by time iteration with an initial input solution . This procedure leads to good result in the effective one-  and two-space-variable  cases and we extend it here to the full three-dimensional case. In the two-space-variable case a three-step procedure  is first used to separate the Hamiltonian into three parts before applying the Crank-Nicholson scheme in two directions. In the three-space variable case we consider a four-step procedure where the Hamiltonian is broken up into four parts. We illustrate the FDCN method in three dimensions for the solution of the GP equation with a three-dimensional standing-wave optical-lattice potential.

In section 2 we describe briefly the three-dimensional, time-dependent GP equation with an anisotropic trap. The PSRK and FDCN methods are described in section 3. In section 4 we report the numerical results for the wave function and energy for different symmetries and nonlinearities as well as an account of our study of resonance in the anisotropic case due to an oscillating scattering length. In section 5 we present the solution of the GP equation for an optical-lattice potential using the FDCN method and finally, in section 6 we present a discussion of our study.

## 2 Nonlinear Gross-Pitaevskii Equation

At zero temperature, the time-dependent Bose-Einstein condensate wave function $`\mathrm{\Psi }(𝐫;\tau )`$ at position $`𝐫`$ and time $`\tau `$ may be described by the following mean-field nonlinear GP equation

$`\left[{\displaystyle \frac{\mathrm{}^2^2}{2m}}+V(𝐫)+gN_0|\mathrm{\Psi }(𝐫;\tau )|^2i\mathrm{}{\displaystyle \frac{}{\tau }}\right]\mathrm{\Psi }(𝐫;\tau )=0.`$
Here $`m`$ is the mass and $`N_0`$ the number of atoms in the condensate, $`g=4\pi \mathrm{}^2a/m`$ the strength of interatomic interaction, with $`a`$ the atomic scattering length. The normalization condition of the wave function is $`𝑑𝐫|\mathrm{\Psi }(𝐫;\tau )|^2=1.`$

The three-dimensional trap potential is given by $`V(𝐫)=\frac{1}{2}m(\omega _0^2\overline{x}^2+\omega _y^2\overline{y}^2+\omega _z^2\overline{z}^2)`$, where $`\omega _x\omega _0`$, $`\omega _y`$, and $`\omega _z`$ are the angular frequencies in the $`x`$, $`y`$ and $`z`$ directions, respectively, and $`𝐫(\overline{x},\overline{y},\overline{z})`$ is the radial vector. The wave function can be written as $`\mathrm{\Psi }(𝐫;\tau )=\psi (\overline{x},\overline{y},\overline{z},\tau )`$. After a transformation of variables to dimensionless quantities defined by $`x=\sqrt{2}\overline{x}/l`$, $`y=\sqrt{2}\overline{y}/l`$, $`z=\sqrt{2}\overline{z}/l`$, $`t=\tau \omega _0,`$ $`l\sqrt{(\mathrm{}/m\omega _0)}`$ and $`\varphi (x,y,z,t)=\psi (\overline{x},\overline{y},\overline{z},\tau )(l^3/\sqrt{8})^{1/2}`$, the GP equation becomes

$`\left[{\displaystyle \frac{^2}{x^2}}{\displaystyle \frac{^2}{y^2}}{\displaystyle \frac{^2}{z^2}}+{\displaystyle \frac{x^2+\kappa ^2y^2+\nu ^2z^2}{4}}+𝒩\left|\varphi (x,y,z,t)\right|^2i{\displaystyle \frac{}{t}}\right]\varphi (x,y,z,t)=0,`$ (2.1)
where $`𝒩=8\sqrt{2}\pi n`$, $`\kappa =\omega _y/\omega _0`$ and $`\nu =\omega _z/\omega _0`$ with nonlinearity $`n=N_0a/l`$. The normalization condition for the wave function is

$`{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑x{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑y{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑z|\varphi (x,y,z,t)|^2=1.`$ (2.2)
## 3 Numerical Methods

### 3.1 Pseudospectral Runge-Kutta (PSRK) method

First we describe the PSRK method  for the one-dimensional GP equation in some detail and then indicate the necessary changes for the three-dimensional case. The one-dimensional GP equation is obtained by eliminating kinetic energy (derivative) terms in $`y`$ and $`z`$, setting $`\kappa =\nu =0`$ and eliminating the $`y`$ and $`z`$ dependence of $`\varphi `$ in (2.1), e.g.,

$`\left[{\displaystyle \frac{^2}{x^2}}+{\displaystyle \frac{x^2}{4}}+𝒩\left|\varphi (x,t)\right|^2i{\displaystyle \frac{}{t}}\right]\varphi (x,t)=0,`$ (3.1)
with the normalization $`_{\mathrm{}}^{\mathrm{}}𝑑x|\varphi (x,t)|^2=1.`$

In this method the unknown function $`\varphi (x,t)\varphi (x)`$ is expanded in terms of a set of $`N`$ known interpolating orthogonal functions $`\{f_j(x)\}_{j=0}^{N1}`$ as follows

$$\varphi (x)p_{N1}(x)=\underset{j=0}{\overset{N1}{}}\frac{\alpha (x)}{\alpha (x_j)}f_j(x)\varphi _j,$$
(3.2)
where $`\{x_j\}_{j=0}^{N1}`$ is a set of distinct interpolation nodes, $`\varphi _j\varphi (x_j)`$, $`\alpha (x)`$ is a weight function, and the functions $`\{f_j(x)\}_{j=0}^{N1}`$ satisfy $`f_j(x_k)=\delta _{jk}`$ (the Kronecker delta) and involve orthogonal polynomials of degree $`(N1)`$, so that $`\varphi (x_k)=p_{N1}(x_k),k=0,1,\mathrm{},N1.`$ In this work the interpolating functions $`\{f_j(x)\}_{j=0}^{N1}`$ are the Hermite polynomials $`H_j(x)`$: $`f_j(x)=H_j(x)`$. However, one could use other polynomials, such as, Chebyshev, Laguerre, and Legendre. One could also consider a Fourier (spectral) expansion of the wave function in terms of periodic cosine and sine functions. The Hermite polynomials are the eigenfunctions of the linearized GP equation and hence already satisfy the boundary conditions of the wave function of the GP equation . Consequently, by choosing the Hermite polynomials in the expansion, (3.2) satisfies the proper boundary conditions by construction.

For obtaining the numerical solution, the GP equation (3.1) is defined and solved on the set of grid points $`x_j`$. The solution at any other point is obtained by using the interpolation formula (3.2), or any other convenient interpolation rule. The advantage of the expansion (3.2) is that when it is substituted in the GP equation, the space derivatives operate only on the known analytic functions $`\alpha (x)`$ and $`f_j(x)`$, so that one can define a matrix for the second-order space derivative :

$$D_{k,j}^2=\frac{d^2}{dx^2}\left[\frac{\alpha (x)}{\alpha (x_j)}f_j(x)\right]_{x=x_k},$$
(3.3)
and the numerical differentiation process may therefore be performed as the matrix-vector product: $`_{j=0}^{N1}D_{k,j}^2\varphi _j`$. Consequently, the partial differential equation (3.1) is reduced to a set of coupled ODEs in the time variable $`t`$ involving $`\varphi _j,j=0,1,\mathrm{},(N1)`$. In this way we obtain a set of ODEs by considering the original equations on a suitable set of discretization points (the roots of Hermite polynomials). One could have also used a Galerkin procedure  and projected the equations, essentially by integrating them against Hermite polynomials for that purpose. However, we do not explore this possibility in this paper.

For solving the set of ODEs we use the adaptive step-size control based on the embedded Runge-Kutta formulas due to Fehlberg , which gives a definite clue about how to modify the step size in order to achieve a desired accuracy in a controlled way. For orders $`M`$ higher than four of the Runge-Kutta formula, evaluation of more than $`M`$ functions (though never more than $`M+2`$) is required. This makes the classic fourth order method requiring the evaluation of four functions more economic. Fehlberg suggested a fourth-order and a fifth-order method each requiring the evaluation of six functions. The difference between the results of these two gives the error $`\delta `$ in the fourth-order method with a step size $`h`$, where $`\delta `$ scales as $`h^5:`$ $`\delta h^5`$. This scaling immediately gives the factor by which the step size $`h`$ should be reduced, so that a desired $`\delta `$ can be obtained. The detailed fourth-order and fifth-order Runge-Kutta formulas of Fehlberg are given in . We use these formulas with the constants given by Cash and Karp also tabulated in . For the present problem we find that the use of Cash-Karp constants in the Fehlberg formulas leads to more accurate results than the original constants due to Fehlberg.

Next we specialize to the case of Hermite polynomials used to generate the solution. Hermite polynomials are very convenient in this case  as the solutions of the linear GP equation (3.1) with $`𝒩=0`$ are Gaussian-type functions. They also satisfy the correct boundary conditions of the wave functions. In this case $`x_j`$ are the roots of $`H_{N1}(x)`$: $`H_{N1}(x_j)=0,j=0,1,\mathrm{},N1.`$ The roots can be found by diagonalizing a tridiagonal symmetric Jacobi matrix as described in  or otherwise. The weight functions are taken as $`\alpha (x)=\mathrm{exp}(x^2/2)`$, such that expansion (3.2) becomes

$$p_{N1}(x)=\underset{j=0}{\overset{N1}{}}\frac{\mathrm{exp}(x^2/2)}{\mathrm{exp}(x_j^2/2)}f_j(x)\varphi _j,$$
(3.4)
where $`f_j(x)`$ are taken as

$$f_j(x)=\frac{H_{N1}(x)}{H_{N1}^{}(x_j)(xx_j)},$$
(3.5)
where prime refers to derivative with respect to $`x`$. Consequently, the differentiation matrix can be obtained from

$$D_{k,j}^2=\frac{1}{\alpha (x_j)H_{N1}^{}(x_j)}\frac{d^2}{dx^2}\left[\frac{\mathrm{exp}(x^2/2)H_{N1}(x)}{(xx_j)}\right]_{x=x_k}$$
(3.6)
and calculated using an algorithm described in .

Using the differentiation matrix (3.6), the GP equation is discretized. The grid points are the roots of the Hermite polynomial $`H_{N1}(x_j)=0`$. However, the actual $`x_j`$ values employed are obtained by scaling these roots by a constant factor so that most of the roots fall in the region where the condensate wave function is sizable and only a few points are located in the region where the wave function is negligible. For the spherically symmetric case $`\omega _0=\omega _y=\omega _z,`$ the discretization mesh in the three directions are identical. For the anisotropic cases, in general, the discretization points in the three directions are different from each other.

Though the passage from the one-dimensional to three-dimensional PSRK method is formally straightforward, it involves nontrivial computational steps. The unknown function $`\varphi (x,y,z,t)\varphi (x,y,z)`$ is expanded in terms of a set of $`N`$ known interpolating orthogonal functions $`\{f_j(x)\}_{j=0}^{N1}`$ as follows

$`\varphi (x,y,z)p_{N1}(x,y,z)={\displaystyle \underset{i=0}{\overset{N1}{}}}{\displaystyle \underset{j=0}{\overset{N1}{}}}{\displaystyle \underset{k=0}{\overset{N1}{}}}{\displaystyle \frac{\alpha (x)}{\alpha (x_i)}}{\displaystyle \frac{\alpha (y)}{\alpha (y_j)}}{\displaystyle \frac{\alpha (z)}{\alpha (z_k)}}f_i(x)f_j(y)f_k(z)\varphi _{ijk},`$ (3.7)
where $`\{x_j\}_{j=0}^{N1}`$ is a set of distinct interpolation nodes, $`\varphi _{ijk}\varphi (x_i,x_j,x_k)`$, and the functions $`\alpha `$ and $`f`$ are defined as in the one-dimensional case, so that $`\varphi (x_i,y_j,z_k)=p_{N1}(x_i,y_j,z_k)`$. The differentiation matrices along $`x`$, $`y`$ and $`z`$ directions can be defined via (3.6) as in the one-dimensional case. Consequently, the partial differential GP equation (2.1) in three space variables is transformed to a set of ODEs in time variable on the grid points $`x_i`$, $`y_j`$ and $`z_k`$, which is solved by the adaptive step-size controlled Runge-Kutta method. The wave function at any point is then calculated using the interpolation formula (3.7).

### 3.2 Finite-difference Crank-Nicholson (FDCN) method

The GP equation (2.1) has the form of the following nonlinear Schrödinger equation

$$i\frac{\varphi }{t}=H\varphi ,$$
(3.8)
where the Hamiltonian $`H`$ contains the different linear and nonlinear terms including the spatial derivatives. We solve this equation by time iteration after discretization in space and time using the finite difference scheme . This procedure leads to a set of algebraic equations. In the present split-step method the iteration is conveniently performed in four steps by breaking up the full Hamiltonian into different derivative and nonderivative parts: $`H=H_1+H_2+H_3+H_4`$, where

$`H_1`$ $`=`$ $`{\displaystyle \frac{x^2+\kappa ^2y^2+\nu ^2z^2}{4}}+𝒩|\varphi (x,y,z,t)|^2,`$ (3.9)
$`H_2`$ $`=`$ $`{\displaystyle \frac{^2}{x^2}},H_3={\displaystyle \frac{^2}{y^2}},H_4={\displaystyle \frac{^2}{z^2}}.`$ (3.10)
The time variable is discretized as $`t_n=n\mathrm{\Delta }`$ where $`\mathrm{\Delta }`$ is the time step. The solution is advanced first over the time step $`\mathrm{\Delta }`$ at $`t_n`$ by solving (3.8) with $`H=H_1`$ to produce an intermediate solution $`\varphi ^{n+1/4}`$ from $`\varphi ^n`$, where $`\varphi ^n`$ is the discretized wave function at $`t_n`$. This propagation is performed essentially exactly for small $`\mathrm{\Delta }`$ via

$$\varphi ^{n+1/4}=𝒪_{\text{nd}}(H_1)\varphi ^n\mathrm{exp}(i\mathrm{\Delta }H_1)\varphi ^n,$$
(3.11)
where $`𝒪_{\text{nd}}(H_1)`$ denotes time evolution with $`H_1`$ and the suffix ‘nd’ refers to non-derivative terms. Next we perform the time propagation corresponding to the operators $`H_i,i=2,3,4`$ successively via the following semi-implicit Crank-Nicholson schemes :

$`\varphi ^{n+2/4}=𝒪_{\text{CN}}(H_2)\varphi ^{n+1/4}{\displaystyle \frac{1i\mathrm{\Delta }H_2/2}{1+i\mathrm{\Delta }H_2/2}}\varphi ^{n+1/4}`$ (3.12)
$`\varphi ^{n+3/4}=𝒪_{\text{CN}}(H_3)\varphi ^{n+2/4}{\displaystyle \frac{1i\mathrm{\Delta }H_3/2}{1+i\mathrm{\Delta }H_3/2}}\varphi ^{n+2/4}`$ (3.13)
$`\varphi ^{n+1}=𝒪_{\text{CN}}(H_4)\varphi ^{n+3/4}{\displaystyle \frac{1i\mathrm{\Delta }H_4/2}{1+i\mathrm{\Delta }H_4/2}}\varphi ^{n+3/4},`$ (3.14)
where $`𝒪_{\text{CN}}(H_i)`$ denotes time evolution with $`H_i`$ and the suffix ‘CN’ refers to Crank-Nicholson. Hence the final solution at time $`t_{n+1}`$ is obtained from

$`\varphi ^{n+1}=𝒪_{\text{CN}}(H_4)𝒪_{\text{CN}}(H_3)𝒪_{\text{CN}}(H_2)𝒪_{\text{nd}}(H_1)\varphi ^n.`$ (3.15)
The details of the Crank-Nicholson discretization scheme can be found in . The advantages of the above split-step method are the following. First, the error involved in splitting the Hamiltonian is proportional to $`\mathrm{\Delta }^2`$ and can be neglected for small $`\mathrm{\Delta }`$. A considerable fraction $`H_1`$ of the Hamiltonian is treated fairly accurately without mixing with the Crank-Nicholson propagation. This method can deal with a large nonlinear term accurately and lead to stable and accurate converged result.

### 3.3 Calculational Details

In both the FDCN and PSRK methods the time iteration is started with the following normalized ground-state solution of the linear GP equation (2.1) with $`𝒩=0`$:

$`\varphi (x,y,z)=\left[{\displaystyle \frac{\kappa \nu }{8\pi ^3}}\right]^{1/4}\mathrm{exp}[(x^2+\kappa y^2+\nu z^2)/4].`$ (3.16)
The norm of the wave function is conserved after each iteration due to the unitarity of the time evolution operator. However, it is of advantage to reinforce numerically the proper normalization of the wave function after several (100) time iterations in order to improve the precision of the result. Typical time step employed in the calculation is $`\mathrm{\Delta }=0.001`$. During the iteration the coefficient $`n=N_0a/l`$ of the nonlinear term is increased from 0 at each step by $`\mathrm{\Delta }_1=0.001`$ until the final value of nonlinearity $`n`$ is attained. This corresponds to the final solution. Then several thousand time iterations of the equation were performed until a stable result is obtained.

For large nonlinearity, the Thomas-Fermi (TF) solution of the GP equation is a better approximation to the exact result  than the harmonic oscillator solution (3.16). In that case it might be advisable to use the TF solution as the initial trial input to the GP equation with full nonlinearity and consider time iteration of this equation without changing the nonlinearity. This time iteration is to be continued until a converged solution is obtained. However, in all the calculations reported in this paper only (3.16) is used as trial input. The use of initial TF solution did not lead to satisfactory result.

## 4 Numerical Results with the PSRK Method

### 4.1 Wave Function and Energy

The present method relies on time evolution and is suitable for both stationary and time-evolution problems. The stationary problems are governed by a wave function with trivial time dependence $`\varphi (x,y,z,t)=\varphi (x,y,z)\mathrm{exp}(i\mu t)`$, where $`\mu `$ is a real energy parameter. Thus the stationary wave function $`\varphi `$ and the parametric energy $`\mu `$ (the chemical potential) can be extracted from the evolution of the time-dependent GP equation over a macroscopic interval of time . Here we present results for the chemical potential of several BECs in three dimensions for spherically symmetric, axially symmetric and anisotropic cases. These traps with different geometries for <sup>23</sup>Na have been employed in experiments as well as in a time-independent solution of the three-dimensional GP equation . The completely anisotropic trap employs the parameters $`\omega _0^A=354\pi `$ rad/s, $`\kappa =\sqrt{2},\nu =2`$ as in the experiment of Kozuma et al. . The cylindrically symmetric trap parameters are $`\omega _0^C=33.86\pi `$ rad/s, $`\kappa =\nu =13.585`$ as in . The spherically symmetric trap parameters are $`\omega _0^S=87`$ rad/s, $`\kappa =\nu =1`$ . We employ the scattering length of $`a=52a_0`$ of Na, where $`a_0`$ is the Bohr radius .

In all calculations reported in this section we used 21 Hermite polynomials each in $`x`$, $`y`$ and $`z`$ directions so that we shall be dealing with a wave function in the form of a cubic array of dimension $`21\times 21\times 21`$. The maxima of $`x`$, $`y`$ and $`z`$ in discretization were chosen consistent with the trap parameters. Typical maxima $`|x|_{\text{max}},|y|_{\text{max}},`$ and $`|z|_{\text{max}}`$ are of the order of 8 for the spherical and anisotropical cases, although a smaller $`|y|_{\text{max}}`$ and $`|z|_{\text{max}}`$ ($`3`$) together with a larger $`|x|_{\text{max}}`$ $`(10)`$ have been used in the axially symmetric case because of large distortion of trap parameters in that case ($`\kappa =\nu =13.585`$).

In table 1 we list the chemical potentials for different symmetries obtained from the PSRK approach as a function of number of atoms $`N_0=2^q`$ in the condensate and compare them with those of a calculation based on a discrete variable representation of the time-independent GP equation . In the present calculation the chemical potentials are extracted from results of time evolution of the GP equation. The results were evaluated at a space point near the center of the BEC and averaged over several samples of calculation. The error (standard deviation) of the time averaged chemical potential is typically of the order of $`0.2\%`$. Considering that the present approach is based on time evolution, the precision is quite satisfactory $``$ about less than a percentage point of discrepancy when compared to results of . It was most difficult to obtain good convergence in the cylindrical case with highly distorted trap ($`\kappa =\nu =13.585`$). A more carefully chosen values of the maxima $`|x|_{\text{max}}`$ $`(10),`$ and $`|y|_{\text{max}}=|z|_{\text{max}}`$ $`(3)`$ were needed in this case.

In the following we study the convergence of the PSRK method in the anisotropical case with $`q=12`$ calculated with $`|x|_{\text{max}}=9`$, $`|y|_{\text{max}}=6.5`$ and $`|z|_{\text{max}}=5`$ which seems to be optimal and were found after some experimentation. These values were taken to be roughly five times the root mean square (rms) sizes $`x,y,z_{\text{rms}}`$ in the respective directions. Here we consider the convergence of these rms sizes. After the proper nonlinearity is introduced in the GP equation, if we continue to integrate using the Runge-Kutta ODE solver routine the wave function and the rms sizes fluctuate a little. This type of fluctuation is common to time stepping methods for a partial differential equation. To quantify this fluctuation we calculate the rms sizes over 250 successive samples generated after 20 time steps of 0.01 each in the Runge-Kutta ODE routine. This corresponds to a total time interval of $`250\times 20\times 0.01`$ or 50 units. Then we calculate the mean rms sizes and the standard deviations which we show in table 2 for different number of expansion points $`N`$ in (3.2). The numerical error in the rms sizes is typically less than one percent for $`N>10`$ and the convergence is quite satisfactory. Considering that we are solving a partial differential equation in four variables this error is small. It will be difficult to obtain similar precision in the FDCN method in this problem even with a significantly larger number of space discretization points ($`N>100`$) in each direction. The FDCN method is no match to the PSRK method in this problem with a smooth potential. However, the FDCN method seems to be very suitable for a rapidly varying potential, such as the optical-lattice potential considered in section 5, which requires a large number of equally distributed spatial discretization points for a faithful reproduction of the potential. Although, the present PSRK method yields satisfactory result for the stationary problem, its main advantage lies in its ability to tackle time-dependent problems as we shall see in the following.

### 4.2 Resonance Dynamics

The appearance of a resonance in the oscillation of a BEC due to a periodic variation of the scattering length has been postulated recently in spherically  and axially symmetric traps . Here we extend this investigation to the more realistic and complicated case of an anisotropic trap.

In the case of a damped classical oscillator under the action of an external periodic force a resonance appears when the frequency of the external force is equal to or a multiple of the natural frequency of oscillation of the system. The natural frequency of oscillation of a trapped three-dimensional BEC along a certain direction ($`x`$, $`y`$, or $`z`$) is twice the trapping frequency in that direction . For a trapped BEC subject to a nonlinear external force due to a periodic variation in the scattering length given by $`n=n_0\mathrm{sin}(\mathrm{\Omega }t)`$, in analogy with the damped classical oscillator a resonance in the oscillation of the BEC along a particular direction is expected when $`\mathrm{\Omega }`$ equals the natural frequency of oscillation in that direction or a multiple of this frequency. In the present numerical study of the anisotropic case we find that this is indeed the case.

To investigate the phenomenon of resonance we solve (2.1) with $`𝒩=8\pi \sqrt{2}n_0\mathrm{sin}(\mathrm{\Omega }t)`$ with $`n_0=0.3`$. Actually, $`n_0`$ has to be less than a critical value $`n_{\text{crit}}`$ ($`0.55`$) in order to avoid collapse for attractive interaction . The actual critical value depends on the asymmetry parameters $`\kappa `$ and $`\nu `$ of the harmonic trap. For the spherically symmetric case $`\kappa =\nu =1`$, $`n_{\text{crit}}=`$ 0.575 . To study the resonance dynamics we use 31 Hermite polynomials in each of $`x`$, $`y`$, and $`z`$ directions. The maximum values of $`|x|`$, $`|y|`$, and $`|z|`$ employed in space discretization are each 15. The resonance is best studied by plotting the rms sizes $`x,y,z_{\text{rms}}`$ vs. time. In this study we employ the trap parameters $`\omega _x\omega _0=354\pi `$ rad/s, $`\kappa =\sqrt{2}`$ and $`\nu =2`$ , so that the unit of time is $`\omega _0^1=0.9`$ ms and of length is $`l/\sqrt{2}=1.115`$ $`\mu `$m. For values of $`\mathrm{\Omega }`$ off the resonance the rms sizes exhibit oscillation of very small amplitude. For resonance frequencies $`\mathrm{\Omega }`$, the rms sizes execute oscillation of large amplitude.

To illustrate the resonance we plot in figures 1 (a), (b), and (c) $`x,y,z_{\text{rms}}`$ vs. time for $`\mathrm{\Omega }=2`$, $`2\sqrt{2}`$, and 4, respectively. These values of $`\mathrm{\Omega }`$ correspond to the natural frequency of oscillation of the condensate along the $`x`$, $`y`$, and $`z`$ directions, respectively . In addition, $`\mathrm{\Omega }=4`$ is also twice the natural frequency of oscillation of the condensate along the $`x`$ direction. So for $`\mathrm{\Omega }=2`$ and $`2\sqrt{2}`$ the rms values of $`x`$ and $`y`$ execute resonance oscillation as shown in figures 1 (a) and (b), respectively. At resonance in a particular direction, the corresponding dimension increases with time in a oscillatory fashion. The rms values of the other components do not show resonance. For $`\mathrm{\Omega }=4`$, rms values of both $`x`$ and $`z`$ exhibit resonance. In case of $`z_{\text{rms}}`$ this corresponds to the lowest harmonic and for $`x_{\text{rms}}`$ this corresponds to the first excited state. However, the behavior of $`x_{\text{rms}}`$ and $`z_{\text{rms}}`$ are different in this case.

## 5 Numerical Results with the FDCN Method

For a smooth trapping potential the PSRK method discussed in the last section yields excellent result with a smaller number of spatial discretization points which are unevenly distributed (at the scaled roots of the Hermite polynomial) compared to the FDCN method employing a relatively large number of evenly distributed spatial discretization points. For smooth potentials the CPU time in the PSRK method could be even an order of magnitude smaller than that in the FDCN method. However, the FDCN method has advantage in the case of a rapidly varying trapping potential requiring a large number of evenly distributed spatial discretization points for a proper description of the trapping potential. One such potential is the optical-lattice trapping potential recently used in BEC experiments in one  and three  dimensions.

The optical-lattice potential created with the standing-wave laser field of wavelength $`\lambda `$ is given by $`V_{\text{opt}}=V_0E_R_{i=1}^3\mathrm{sin}^2(k_Lx_i)`$, with $`E_R=\mathrm{}^2k_L^2/(2m)`$, $`k_L=2\pi /\lambda `$, and $`V_0`$ the dimensionless strength of the optical-lattice potential governed by the intensity of the laser . In terms of the dimensionless laser wave length $`\lambda _0=\sqrt{2}\lambda /l`$ and the dimensionless standing-wave energy parameter $`E_R/(\mathrm{}\omega )=4\pi ^2/\lambda _0^2`$, $`V_{\text{opt}}`$ is given by

$$\frac{V_{\text{opt}}}{\mathrm{}\omega }=V_0\frac{4\pi ^2}{\lambda _0^2}\underset{i=1}{\overset{3}{}}\mathrm{sin}^2\left(\frac{2\pi }{\lambda _0}x_i\right).$$
(5.1)
In the actual experiment  this dimensionless standing-wave optical-lattice potential is superposed on the harmonic trapping potential of (2.1) and we present the solution of (2.1) under the action of spherically symmetric harmonic ($`\kappa =\nu =1`$) as well as the optical-lattice potential (5.1).

To calculate the wave function we discretize the GP equation with time step 0.001 and space step 0.1 along $`x`$, $`y`$, and $`z`$ directions spanning the cubic region $`4<x,y,z<4`$. Consequently, the wave function is defined on a cubic array of dimension $`80\times 80\times 80`$. Starting the time iteration with the known harmonic oscillator solution for nonlinearity $`n=0`$ and $`V_0=0`$, the nonlinearity $`n`$ and the optical-lattice strength $`V_0`$ are slowly increased by equal amount in $`1000n`$ steps of time iteration until the desired values of $`n`$ and $`V_0`$ are obtained. Then without changing any parameters the solution is iterated several thousand times so that a converged solution independent of the initial inputs is obtained. To illustrate the present method we calculate the wave function for the ground state for nonlinearity $`nN_0a/l=10`$, laser wave length $`\lambda _0=1`$ and optical-lattice strength $`V_0=10`$. Two views of the three-dimensional contour plot of the central part of the wave function on a cubic lattice of size $`3\times 3\times 3`$ are shown in figures 2 (a) and (b) as seen from two different angles. The droplets of BEC at each pit of the optical-lattice potential can be identified in figures 2. There are about 10 occupied cites in each of $`x`$, $`y`$, and $`z`$ directions of which the central part is shown in figure 2.

## 6 Conclusion

In this paper we propose and implement a PSRK method  and contrast it with the FDCN method  for the numerical solution of the time-dependent nonlinear GP equation under the action of a three-dimensional trap by time iteration. In the PSRK method the unknown wave function is expanded in a set of known polynomials (Hermite). Consequently, the partial differential GP equation in three space and time variables becomes a set of coupled ODEs in time for the unknown coefficients, which is solved by a fourth-order adaptive step-size control Runge-Kutta method . In the FDCN method the full Hamiltonian is split into the derivative and nonderivative parts. In this fashion the time propagation with the nonderivative parts can be treated very accurately. The time derivative part is treated by the Crank-Nicholson scheme in three independent steps. Both methods lead to stable and accurate results. The final result remains stable for thousands of time iteration of the GP equation.

We applied the PSRK method for the numerical study of certain stationary and time-evolution problems. We solved the GP equation for spherically symmetric, axially symmetric and anisotropic cases and calculated the chemical potential for different nonlinearities. The results compare well with those obtained by a time-independent approach . The two sets of energy values agree to within a fraction of a percentage point. The numerical error in the method is found to be less than one percent with a small number of expansion functions ($`N20`$). The PSRK method was also used to study the resonance dynamics of an anisotropic BEC under the action of periodic sinusoidal variation of the scattering length. When this period of oscillation coincides with the natural frequency of oscillation of the BEC along $`x`$, $`y`$, or $`z`$ directions or a multiple thereof the corresponding rms size executes resonant oscillation . Using the FDCN method we study the ground state of a BEC under the combined action of a harmonic and a periodic optical-lattice trapping potential in three space dimensions.

The domain of applicability of the PSRK and FDCN methods seems to be complementary rather than overlapping. The PSRK method is more efficient and economic for smooth potentials where a relatively small number of expansion functions and unevenly distributed spatial discretization points seems to be adequate. The FDCN method employs a large number of evenly distributed spatial discretization points and is suitable for a rapidly varying potential, such as the optical-lattice potential, favoring such a distribution.

The work is supported in part by the Conselho Nacional de Desenvolvimento Científico e Tecnológico and Fundação de Amparo à Pesquisa do Estado de São Paulo of Brazil. Part of the work of PM is supported by the Department of Science and Technology, Government of India.

## References


text:
# CEEX Exponentiation in QEDWork partly supported by the European Community’s Human Potential Programme under contract HPRN-CT-2000-00149 “Physics at Colliders”, by Polish Government grants KBN 5P03B09320 and 2P03B00122, and by NATO grant PST.CLG.977751.

## Abstract

The aim is to summarize briefly on: (a) Main features of YFS Coherent Exclusive EXponentiation in QED, (b) Examples of recent results relevant for LEP/LC physics program.

This presentation is based mainly on three recent papers . Standard Model calculations for LEP with YFS exponentiation are implemented in many Monte Carlo event generators: for $`e^+e^{}f\overline{f}+n\gamma `$, $`f=\tau ,\mu ,d,u,s,c`$ process: (1a) YFS1 (1987-1989) $`𝒪(\alpha ^1)_{exp}`$ ISR , (1b) YFS2$``$KORALZ (1989-1990) $`𝒪(\alpha ^1+h.o.LL)_{exp}`$ ISR , (1c) YFS3$``$KORALZ (1990-1998) $`𝒪(\alpha ^1+h.o.LL)_{exp}`$ ISR+FSR , (1d) KKMC (1998-02) $`𝒪(\alpha ^2+h.o.LL)_{exp}`$ ISR+FSR+Interf. ($`d\sigma /\sigma =0.2\%`$) ; for the low angle Bhabha process $`e^+e^{}e^+e^{}+n\gamma `$ for $`\theta <6^{}`$ the following ones: (2a) BHLUMI 1.x, (1987-1990) $`𝒪(\alpha ^1)_{exp}`$ , (2b) BHLUMI 2.x, (1990-1996) $`𝒪(\alpha ^1+h.o.LL)_{exp}`$ ($`d\sigma /\sigma =0.07\%`$), ; for the large angle Bhabha $`e^+e^{}e^+e^{}+n\gamma `$ for $`\theta >6^{}`$ (3) BHWIDE (1994-1998) $`𝒪(\alpha ^1+h.o.LL)_{exp}`$ ; and finally for $`WW`$ boson production processes $`e^+e^{}W^+W^{}+n\gamma `$, $`W^\pm f\overline{f}`$ (4a) KORALW (1994-2001) ISR YFS LL ($`d\sigma /\sigma =2\%`$)  and (4b) YFS3WW (1995-2001) YFS expon. + Leading Pole Approx. ($`d\sigma /\sigma =0.4\%`$) .

The detailed description of the Coherent Exclusive Exponentiation (CEEX) can be found in ref. , along with the description of the older Exclusive Exponentiation (EEX), closer to the original Yennie-Frautschi-Suura (YFS) scheme . What is in a nutshell YFS exponentiation? The main steps in YFS exponentiation are:

* Reorganization of the perturbative complete $`𝒪(\alpha ^{\mathrm{}})`$ series such that IR-finite $`\overline{\beta }`$ components are isolated (factorization theorem).
* Truncation of the IR-finite $`\overline{\beta }`$s to finite $`𝒪(\alpha ^n)`$ and calculation of them from Feynman diagrams recursively.

Apart from disentangling properly ultraviolet divergences from the infrared (IR) ones YFS deals consistently with the overlapping and sub-leading IR singularities. To illustrate this point let is show the “map of IR singularities” for $`2\gamma `$ distribution $`D_2(k_1,k_2)`$ in the plane of the energies of two real photons:

The point $`R_{12}`$ is the location of the leading double IR singularity and the lines $`R_1`$ and $`R_2`$ host non-leading single IR singularities. The above topological picture translates into a decomposition of the differential distributions: $`D_0(p_{f_1},p_{f_2},p_{f_3},p_{f_4})=\overline{\beta }_0(p_{f_1},p_{f_2},p_{f_3},p_{f_4}),`$
$`D_1(p_f;𝐤_\mathrm{𝟏})=\overline{\beta }_0(p_f)\stackrel{~}{S}(𝐤_\mathrm{𝟏})+\overline{\beta }_1(p_f;k_1),`$and
$`D_2(𝐤_\mathrm{𝟏},𝐤_\mathrm{𝟐})=\overline{\beta }_0\stackrel{~}{S}(𝐤_\mathrm{𝟏})\stackrel{~}{S}(𝐤_\mathrm{𝟐})`$
$`+\overline{\beta }_1(k_1)\stackrel{~}{S}(𝐤_\mathrm{𝟐})+\overline{\beta }_1(k_2)\stackrel{~}{S}(𝐤_\mathrm{𝟏})+\overline{\beta }_2(k_1,k_2)`$
where typicaly $`p_{f_1}+p_{f_2}=p_{f_3}+p_{f_4}`$ for $`\overline{\beta }_0`$ in $`D_0`$ and $`p_{f_1}+p_{f_2}p_{f_3}+p_{f_4}`$ for $`\overline{\beta }_0`$ in $`D_1`$ and $`D_2`$. As we see the IR-finite part $`\overline{\beta }_1`$ in the decomposition of the 1-photon $`D_1`$ distribution into IR-leading and subleading parts is present as an element in the subleading 2-photon $`D_2`$ distribution. This nontrivial fact allows one to find out the IR-finite $`\overline{\beta }_i`$ functions order by order from standard perturbative calculations. In this particular case the recursive, order-by-order definition of IR-finite $`\overline{\beta }`$s looks as follows:
$`\overline{\beta }_0(p_{f_1},p_{f_2},p_{f_3},p_{f_4})=D_0(p_{f_1},p_{f_2},p_{f_3},p_{f_4}),`$
$`\overline{\beta }_1(p_f;k_1)=D_1(p_f;𝐤_\mathrm{𝟏})\overline{\beta }_0(p_f)\stackrel{~}{S}(𝐤_\mathrm{𝟏})`$
$`\overline{\beta }_2(k_1,k_2)=D_2(𝐤_\mathrm{𝟏},𝐤_\mathrm{𝟐})\overline{\beta }_0\stackrel{~}{S}(𝐤_\mathrm{𝟏})\stackrel{~}{S}(𝐤_\mathrm{𝟐})\overline{\beta }_1(k_1)\stackrel{~}{S}(𝐤_\mathrm{𝟐})+\overline{\beta }_1(k_2)\stackrel{~}{S}(𝐤_\mathrm{𝟏})`$
Having at hand the IR-free $`\overline{\beta }_i`$ the YFS scheme provides immediately for the multi-photon differential distributions all over the phase space and the phase space is integraded without any approximations using MC method.

Below we show an example of the ISR $`𝒪(\alpha ^1)`$ exponentiated multi-photon fully differential distributions for the CEEX variant of YFS exponentiation which is implemented in terms of the IR-free $`\beta _i`$ objects constructed at the amplitude level, rather than using the spin-summed squared $`\overline{\beta }_i`$s. The CEEX we present here very schematically for the process $`e^{}(p_1,\lambda _1)+e^+(p_2,\lambda _2)f(q_1,\lambda _1^{})+\overline{f}(q_2,\lambda _2^{})+\gamma (k_1,\sigma _1)+\mathrm{}\gamma (k_n,\sigma _n)`$:
$`\sigma =\underset{n=0}{\overset{\mathrm{}}{}}\underset{m_\gamma }{}𝑑\mathrm{\Phi }_{n+2}\underset{\lambda ,\sigma _1,\mathrm{},\sigma _n}{}|e^{\alpha B(m_\gamma )}_{n,\sigma _1,\mathrm{},\sigma _n}^\lambda (k_1,\mathrm{},k_n)|^2`$,
$`_0^\lambda =\widehat{\beta }_0^\lambda `$, $`\lambda `$=fermion helicities,
$`_{1,\sigma _1}^\lambda (k_1)=\widehat{\beta }_0^\lambda 𝔰_{\sigma _1}(k_1)+\widehat{\beta }_{1,\sigma _1}^\lambda (k_1)`$
$`_{2,\sigma _1,\sigma _2}^\lambda (k_1,k_2)=\widehat{\beta }_0^\lambda 𝔰_{\sigma _1}(k_1)𝔰_{\sigma _2}(k_2)`$
$`+\widehat{\beta }_{1,\sigma _1}^\lambda (k_1)𝔰_{\sigma _2}(k_2)+\widehat{\beta }_{1,\sigma _2}^\lambda (k_2)𝔰_{\sigma _1}(k_1)`$
$`_{n,\sigma _1,\mathrm{}\sigma _n}^\lambda (k_1,k_2,\mathrm{}k_n)=\widehat{\beta }_0^\lambda 𝔰_{\sigma _1}(k_1)𝔰_{\sigma _2}(k_2)\mathrm{}𝔰_{\sigma _n}(k_n)`$
$`+\widehat{\beta }_{1,\sigma _1}^\lambda (k_1)𝔰_{\sigma _2}(k_2)\mathrm{}𝔰_{\sigma _n}(k_n)`$
$`+𝔰_{\sigma _1}(k_1)\widehat{\beta }_{1,\sigma _2}^\lambda (k_2)\mathrm{}𝔰_{\sigma _n}(k_n)\mathrm{}`$
$`\mathrm{}+𝔰_{\sigma _1}(k_1)𝔰_{\sigma _2}(k_2)\mathrm{}\widehat{\beta }_{1,\sigma _2}^\lambda (k_2)`$ The $`𝒪(\alpha ^1)`$ IR-finite building blocks are:
$`\widehat{\beta }_0^\lambda =\left(e^{\alpha B_4}_\lambda ^{\mathrm{Born}+\mathrm{Virt}.}\right)|_{𝒪(\alpha ^1)},`$
$`\widehat{\beta }_{1,\sigma }^\lambda (k)=_{1,\sigma }^\lambda (k)\widehat{\beta }_0^\lambda 𝔰_\sigma (k)`$
It is to be stressed that everything here was done in terms of the spin $``$-amplitudes! Multiphoton differential distributions are positive by construction (which is not the case for EEX). In KKMC the above CEEX is implemented up to $`𝒪(\alpha ^2)`$ for ISR and FSR111Multiperipheral contributions and non-IR ISR$``$FSR interferences are still incomplete at the $`𝒪\left(\alpha ^2\right)`$ level. They are unimportant for most of LEP observables..

The example of the CEEX aplication is shown in the figure above from ref.  where we present the “state of art” comparison of KK MC vs. Zfitter for ISR+FSR, including IFI=ISR$``$FSR. We conclude that the precision $`d\sigma /\sigma =0.2\%`$ was reached, which is enough for LEP2 data analysis. Similar high precision was obtained for the forward-backward asymmetry .

Recently, the CEEX matrix element was extended in KKMC to the neutrino channel. The contributions from following diagrams are now included:

and also all double photon emission contributions feature the exact matrix element. We summarize below on the improvements in KKMC for the neutrino channel.

* The systematic error is estimated to be 1.3% for $`\nu _e\overline{\nu }_e\gamma `$ and 0.8% for $`\nu _\mu \overline{\nu }_\mu \gamma `$ and $`\nu _\tau \overline{\nu }_\tau \gamma `$.
* For observables with two observed photons we estimate the uncertainty to be about 5%.
* These new improved results were obtained thanks to the inclusion of non-photonic electroweak corrections of the ZFITTER package and due to newly constructed, exact, single and double emission photon amplitudes in the KKMC for the contribution with the $`t`$-channel $`W`$ exchange.
* The virtual corrections for the $`W`$ exchange are at present introduced in the approximated form. The exponentiation scheme CEEX is the same as in the original KKMC program

In the first version of KKMC of ref.  the $`𝒪(\alpha ^2)`$ next-to-leading contributions were incomplete. In ref.  the $`𝒪(\alpha ^2)`$ contribution from the following diagrams were evaluated

and partly included in the development of KKMC – and the remaining contributions were found to be negligible. The contributions in questions  are due to one real and one virtual photon emission and the new numerical results are shown in the following figure

where the new calculation is shown to agree with the older calculations of the refs. . NB. These contributions were also recently (re)calculated in ref. for the purpose of the low energy electron collider experiments.

Conclusions The YFS inspired EEX and CEEX schemes are successful examples of the Monte Carlos based directly on the factorization theorem (albeit for IR soft case for Abelian QED only). They Work well in practice: KORALZ, BHLUMI, YWSWW3, BHWIDE, KK MC. Extension (as far as possible) to all collinear singularities would be very desirable and practically important! The KKMC program is extended to the neutrino channel. Missing fully differential $`2f+1\gamma _{virt}+1\gamma _{real}`$ distributions for $`𝒪\left(\alpha ^2\right)`$ CEEX are now available.


text:
# Three-Neutrino MSW Effect and the Lehmann Mass Matrix

## 1 Introduction

We here review some recent results on analytical solutions of the Mikheyev–Smirnov–Wolfenstein (MSW) effect  for the propagation of three neutrino flavours. Such analytic results have been obtained for both the exponential density  and the linear density .

In the case of an exponential electron density, the three neutrino wave functions can be expressed in terms of generalized hypergeometric functions, $`{}_{2}{}^{}F_{2}^{}`$ . For the linear density, the solutions can be expressed in terms of a Fourier transform . In the case of two flavours, these reduce to parabolic cylinder functions or confluent hypergeometric functions, but of different parameters and arguments .

A particular neutrino mass matrix, originally proposed by Lehmann, Newton and Wu (LNW) for quarks  is also reviewed . Within this model, the current data on atmospheric and solar neutrinos permit a determination of the neutrino masses.

## 2 Exponential density

In a medium where the electron neutrino (described by $`\varphi _1(r)`$) interacts differently from the others, the propagation is given by the equation

$$i\frac{\text{d}}{\text{d}r}\left[\begin{array}{c}\varphi _1\\ \varphi _2\\ \varphi _3\end{array}\right]=\left(\left[\begin{array}{ccc}D(r)& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]+\frac{1}{2p}\left[\begin{array}{ccc}M_{11}^2& M_{12}^2& M_{13}^2\\ M_{21}^2& M_{22}^2& M_{23}^2\\ M_{31}^2& M_{32}^2& M_{33}^2\end{array}\right]\right)\left[\begin{array}{c}\varphi _1\\ \varphi _2\\ \varphi _3\end{array}\right]$$
(1)
The mass matrix is real and symmetric, $`M_{ji}^2=M_{ij}^2(M^2)_{ij}`$, but otherwise arbitrary. (The LNW mass matrix will be discussed later.) Furthermore, $`D(r)=\sqrt{2}G_\mathrm{F}N_e(r)`$, with $`G_\mathrm{F}`$ the Fermi constant and $`N_e(r)`$ the solar electron density.

For the sun, the density  is well approximated by an exponential,

$$N_e(r)=N_e(0)e^{r/r_0},r_00.1\times R_{}.$$
(2)
Introducing a new radial variable: $`u=r/r_0+u_0`$, and performing a rotation on the second and third components, Eq. (1) can be written as

$$i\frac{\text{d}}{\text{d}u}\left[\begin{array}{c}\psi _1(u)\\ \psi _2(u)\\ \psi _3(u)\end{array}\right]=[\begin{array}{ccc}\omega _1+e^u& \chi _2& \chi _3\\ \chi _2& \omega _2& 0\\ \chi _3& 0& \omega _3\end{array}]\left[\begin{array}{c}\psi _1(u)\\ \psi _2(u)\\ \psi _3(u)\end{array}\right].$$
(3)
The eigenvalues of the $`3\times 3`$ matrix

$$\left[\begin{array}{ccc}\omega _1& \chi _2& \chi _3\\ \chi _2& \omega _2& 0\\ \chi _3& 0& \omega _3\end{array}\right]$$
(4)
are denoted $`\mu _1`$, $`\mu _2`$ and $`\mu _3`$, they are the squares of the neutrino masses multiplied by $`r_0/(2p)`$. Together with $`\omega _1`$ and $`\omega _2`$ they control the evolution of the $`\psi _i`$.

Introducing now the variable $`z=ie^u`$, the solutions to Eq. (3) can be expressed in terms of solutions to the third-order ordinary differential equation

$`[(z{\displaystyle \frac{\text{d}}{\text{d}z}}i\mu _1)(z{\displaystyle \frac{\text{d}}{\text{d}z}}i\mu _2)(z{\displaystyle \frac{\text{d}}{\text{d}z}}i\mu _3)`$
$`z(z{\displaystyle \frac{\text{d}}{\text{d}z}}i\omega _2)(z{\displaystyle \frac{\text{d}}{\text{d}z}}i\omega _3)]\psi =0,`$ (5)
namely generalized hypergeometric functions :

$`\psi ^{(1)}=e^{i\mu _1u}{}_{2}{}^{}F_{2}^{}\left[\begin{array}{cc}i(\omega _2\mu _1),& i(\omega _3\mu _1)\\ 1i(\mu _2\mu _1),& 1i(\mu _3\mu _1)\end{array}|ie^u\right]`$ (8)
$`\psi ^{(2)}=e^{i\mu _2u}{}_{2}{}^{}F_{2}^{}\left[\begin{array}{cc}i(\omega _2\mu _2),& i(\omega _3\mu _2)\\ 1i(\mu _1\mu _2),& 1i(\mu _3\mu _2)\end{array}|ie^u\right]`$ (11)
$`\psi ^{(3)}=e^{i\mu _3u}{}_{2}{}^{}F_{2}^{}\left[\begin{array}{cc}i(\omega _2\mu _3),& i(\omega _3\mu _3)\\ 1i(\mu _1\mu _3),& 1i(\mu _2\mu _3)\end{array}|ie^u\right]`$ (14)
The solutions to Eq. (3) are thus

$$\psi _i=C_1\psi _i^{(1)}+C_2\psi _i^{(2)}+C_3\psi _i^{(3)},$$
(15)
where the constants $`C_j`$ are determined by the boundary conditions: $`\psi _1=1`$, $`\psi _2=\psi _3=0`$ at the center of the sun, $`r=0`$.

The $`{}_{2}{}^{}F_{2}^{}`$ functions can be defined in terms of the series expansions

$${}_{2}{}^{}F_{2}^{}\left[\begin{array}{cc}\alpha _1,& \alpha _2\\ \rho _1,& \rho _2\end{array}|z\right]=\underset{k=0}{\overset{\mathrm{}}{}}\frac{(\alpha _1)_k(\alpha _2)_k}{(\rho _1)_k(\rho _2)_k}\frac{z^k}{k!}$$
(16)
where $`(\alpha )_k`$ is a Pochhammer symbol, but (somewhat complicated) asymptotic approximations are actually more useful for numerical work .

## 3 The LNW mass matrix

For quarks, it was found  that a particular, simple texture for the $`d`$ ($`d`$, $`s`$, $`b`$) and $`u`$ ($`u`$, $`c`$, $`t`$) quark mass matrices leads to a good description of the CKM matrix . The mass matrix is assumed to have the form

$$M=\left[\begin{array}{ccc}0& d& 0\\ d& c& b\\ 0& b& a\end{array}\right]$$
(17)
with $`b^2=8c^2`$. This is known as a two-zero texture, but differs from those normally considered (see, e.g., ) in the additional relation between $`b`$ and $`c`$. The eigenvalues are given by $`m_1`$, $`m_2`$, and $`m_3`$, with $`m_1m_3`$.

Actually, in the case of quarks, since there is CP violation, a complex CKM matrix is obtained by replacing the parameter $`d`$ in (17) by $`\pm id`$ for $`u`$ quarks (such that $`M`$ remains Hermitian). The CKM matrix becomes

$$V_{\mathrm{CKM}}=R^\mathrm{T}(u)\mathrm{diag}(i,1,1)R(d),$$
(18)
The Jarlskog determinant  thus obtained is $`J=2.6\times 10^5`$, in good agreement with data.

This same mass matrix has been applied to the case of three neutrinos , and rather good fits to the atmospheric  and solar  neutrino data have been obtained.

The matrix $`M`$ is diagonalized, whereby $`M=RM_{\mathrm{diag}}R^\mathrm{T}`$. For this purpose, the following notation is convenient

$`S_1`$ $``$ $`m_3m_2+m_1,`$
$`=`$ $`a+c`$
$`S_2`$ $``$ $`m_3m_2m_3m_1+m_2m_1,`$
$`=`$ $`8c^2+d^2ac`$
$`S_3`$ $``$ $`m_1m_2m_3`$ (19)
$`=`$ $`ad^2.`$
Eliminating $`d`$ and $`c`$, the resulting cubic equation for the parameter $`a`$ can be written as

$$9a^317S_1a^2+(8S_1^2+S_2)aS_3=0.$$
(20)
A physical solution requires $`a`$ real and positive. This is equivalent to having three real solutions for $`a`$. One of these is negative and two are positive. At any point inside the allowed domain in the $`(m_1/m_3)`$$`(m_2/m_3)`$ plane (See Fig. 1), there are thus two allowed solutions, denoted Solutions 1 and 2.

This diagonalization has to be carried out for both neutrinos and charged leptons, in order to obtain the neutrino mixing matrix

$$U=\left[\begin{array}{ccc}U_{e1}& U_{e2}& U_{e3}\\ U_{\mu 1}& U_{\mu 2}& U_{\mu 3}\\ U_{\tau 1}& U_{\tau 2}& U_{\tau 3}\end{array}\right],$$
(21)
where ($`ϵ=1`$ or $`i`$, depending on whether or not there is CP violation)

$$U=(V_{\mathrm{CKM}}^{\mathrm{}})^{}=R^\mathrm{T}(\mathrm{})\mathrm{diag}(ϵ,1,1)R(\nu )$$
(22)
relates the neutrino mass eigenstates to the flavor states of charged-current interactions:

$$|\nu _e=U_{e1}|\nu _1+U_{e2}|\nu _2+U_{e3}|\nu _3,\mathrm{etc}.$$
(23)
There are four possible solutions to Eq. (20), two for the neutrino sector, and two for the charged lepton sector. Furthermore, the model has two signs (denoted “parities”) to be specified: $`b`$ and $`d`$ can each be either positive or negative (unless $`d`$ is imaginary ). However, only the product of the “$`b`$ parities” in the neutrino and charged lepton sector matters, and similarly for the “$`d`$ parities”, so we may put both “parities” of the charged lepton sector to $`+1`$.

## 4 Fits to data

Let us consider first the atmospheric neutrino data. The Super-Kamiokande results  give $`\mathrm{\Delta }m^2(23)\times 10^3`$ eV, with $`\mathrm{sin}^2(2\theta )1`$. The survival of muon neutrinos is given by

$`P_{\nu _\mu \nu _\mu }(t)`$ $`=`$ $`14[U_{\mu 1}^2U_{\mu 2}^2\mathrm{sin}^2\left({\displaystyle \frac{\mathrm{\Delta }m_{21}^2t}{4p}}\right)+U_{\mu 1}^2U_{\mu 3}^2\mathrm{sin}^2\left({\displaystyle \frac{\mathrm{\Delta }m_{31}^2t}{4p}}\right)`$ (24)
$`+U_{\mu 2}^2U_{\mu 3}^2\mathrm{sin}^2\left({\displaystyle \frac{\mathrm{\Delta }m_{32}^2t}{4p}}\right)],`$
where $`U`$ is the neutrino mixing matrix of Eqs. (21) and (22). In the limit of $`\mathrm{\Delta }m_{21}^2t/4p1`$ this simplifies, and invoking further unitarity, one finds

$$P_{\nu _\mu \nu _\tau }(t)4U_{\mu 3}^2U_{\tau 3}^2\mathrm{sin}^2\left(\frac{\mathrm{\Delta }m_{32}^2t}{4p}\right),$$
(25)
which suggests that one needs $`|U_{\mu 3}U_{\tau 3}|=𝒪(1)`$. This can be achieved within the model (for both Solutions 1 and 2), for $`m_1m_3`$, with also $`m_2`$ small compared with $`m_3`$. Furthermore, the data suggest that the scale $`m_3`$ must be such that $`m_3^2(23)\times 10^3`$ eV, or $`m_3=𝒪(0.05\mathrm{eV})`$.

Fits to atmospheric data confirm this qualitative analysis. Invoking also the solar Cl , Ga , Super-Kamiokande  and SNO  neutrino data, one finds that Solution 2 for the neutrino sector, Solution 1 for the charged lepton sector and “$`b`$ parity”=$`1`$ give good fits for $`m_1m_3`$, with $`m_2`$ also small as compared with $`m_3`$. A $`\chi ^2`$ determined from these different atmospheric and solar survival probabilities leads to good fits (see Fig. 1) with $`m_3`$ about 0.052 eV, $`m_2`$ about 0.01 eV, and $`m_1`$ 0.002–0.004 eV . The “$`d`$ parity” is unimportant.

Figure 1. Fits to the atmospheric and solar neutrino data for Solution 2 (neutrino sector) and 1 (charged leptons) and both “parities” negative. Contours are given at $`\chi ^2=5`$, 10, 15, 20, 25. Sum (right panel) also at 30, 35, 40, 45, 50. The best-fit point (right panel, at the lower left), is marked “low”.

The effect of the CHOOZ data  is to disfavor the region $`m_2=𝒪(m_3)`$, where $`|U_{e3}|`$ is large. But this region is already disfavored by the solar neutrino data, as is seen in the middle panel of Fig. 1. In terms of the more conventional two-flavour analyses for the solar-neutrino sector, these fits roughly correspond to the large-mixing-angle solution.

## 5 Summary

Analytic results for the solutions to the MSW equations for three neutrino flavours are very valuable for a fast scanning over the parameters of some given model for the mass matrix.

The LNW mass matrix is a very constrained model that in the quark sector describes the CKM matrix, and in the neutrino sector gives the mixing in terms of the mass eigenvalues. Applied to the neutrino data, the model gives a very good fit.

The solar neutrino data has also been studied within the same model, using numerical integration methods (no $`{}_{2}{}^{}F_{2}^{}`$’s) . An additional fit was then found at $`m_12.8\times 10^6`$ eV, corresponding to the small-mixing-angle solutions. However, this point is disfavoured by the atmospheric neutrino data, and by the flat electron recoil spectrum .

Acknowledgments

It is a pleasure to thank the organizers of “Beyond the Desert 2002” for creating a very stimulating atmosphere. I should also like to thank R. Gastmans, G. Vigdel and T. T. Wu for very fruitful collaborations. This work is supported in part by the Research Council of Norway.

## References