text:
# SLAC-PUB-8557 hep-ph/0008099 𝒪(𝛼³ln𝛼) corrections to positronium decay rates

## I Introduction

Positronium decays into two and three photons provide an interesting test of bound states Quantum Electrodynamics (QED). While the parapositronium (p-Ps) decay rate is well described by QED, it is known that the decay of orthopositronium (o-Ps) into three photons is still a controversial issue . However, on the theoretical side there has been a major breakthrough recently and both decay rates are now known with $`𝒪(\alpha ^2)`$ accuracy . The next level of precision, i.e. $`𝒪(\alpha ^3)`$ correction, is currently beyond the reach, although parts of the $`𝒪(\alpha ^3)`$ correction that are enhanced by the logarithms of the fine structure constant, can be computed.

The knowledge of these corrections is not of significant phenomenological importance at present since they are much smaller than the experimental accuracy for both p-Ps and o-Ps decays. Nevertheless, from a theoretical viewpoint, it is an interesting problem because at order $`\alpha ^3`$ there are both leading $`𝒪(\alpha ^3\mathrm{ln}^2\alpha )`$ and subleading $`𝒪(\alpha ^3\mathrm{ln}\alpha )`$ corrections. The calculation of the leading $`𝒪(\alpha ^3\mathrm{ln}^2\alpha )`$ corrections is a fairly simple enterprise; it has been done quite some time ago . Here we are interested in subleading logarithmic corrections. In order to compute them, we utilize the nonrelativistic QED in dimensional regularization and we explain how in this framework the logarithmic corrections can be computed using limited amount of information.

Before diving into the description of the calculation, let us summarize our results. For the $`𝒪(\alpha ^3\mathrm{ln}\alpha )`$ corrections to para- and orthopositronium decay rates we find:

$`\mathrm{\Delta }\mathrm{\Gamma }_\mathrm{p}={\displaystyle \frac{\alpha ^3}{\pi }}\mathrm{ln}\alpha \left\{{\displaystyle \frac{367}{90}}+10\mathrm{ln}22A_\mathrm{p}\right\}\mathrm{\Gamma }_\mathrm{p}^{(0)},`$ (1)
$`\mathrm{\Delta }\mathrm{\Gamma }_\mathrm{o}={\displaystyle \frac{\alpha ^3}{\pi }}\mathrm{ln}\alpha \left\{{\displaystyle \frac{229}{30}}+8\mathrm{ln}2+{\displaystyle \frac{A_\mathrm{o}}{3}}\right\}\mathrm{\Gamma }_\mathrm{o}^{(0)}.`$ (2)
The coefficients $`A_{\mathrm{p},\mathrm{o}}`$ describe $`𝒪(\alpha /\pi )`$ corrections to the lowest order decay widths $`\mathrm{\Gamma }_{\mathrm{p},\mathrm{o}}^{(0)}`$. They are

$$A_\mathrm{p}=\frac{\pi ^2}{4}5,A_\mathrm{o}=10.286606(10).$$
(3)
Numerically, the $`𝒪(\alpha ^3\mathrm{ln}\alpha )`$ corrections, Eqs.(1,2), evaluate to $`7.7919\alpha ^3/\pi \mathrm{ln}\alpha \mathrm{\Gamma }_\mathrm{p}^{(0)}`$ for p-Ps and to $`5.517\alpha ^3/\pi \mathrm{ln}\alpha \mathrm{\Gamma }_\mathrm{o}^{(0)}`$ for o-Ps. These corrections are therefore quite comparable with the “leading” $`𝒪(\alpha ^3\mathrm{log}^2\alpha )`$ corrections computed in .

The paper is organized as follows. In the next Section we set up the framework of the calculation. We then continue with detailed discussion of how various contributions to Ps decays at $`𝒪(\alpha ^3\mathrm{ln}\alpha )`$ are computed. In the last Section we present our conclusions.

## II Framework of the calculation

Let us first discuss the framework of our calculation. We work in nonrelativistic QED regularized dimensionally; $`d=32ϵ`$ is the number of spatial dimensions and $`ϵ`$ is the regularization parameter. General features of this technique have been described at length in our previous paper . Here we would like to discuss a new issue which was not considered in  – how logarithmic in $`\alpha `$ corrections can be extracted.

In order to extract logarithms of the fine structure constant in a self-consistent way, we use the fact that the matrix element of any operator in dimensional regularization is a uniform  function of the fine structure constant. This implies that all the dependence on the fine structure constant can be scaled out of any dimensionally regularized matrix element***We stress that this is the feature of dimensional regularization and it is not valid in other regularization schemes.. The scaling, however, should be done in $`d`$ dimensions.

In order to establish the scaling rules, we need to know how different quantities involved in bound state calculations scale with $`\alpha `$. To do that, we rewrite the $`d`$-dimensional Schrödinger equation in “atomic units” familiar from the standard treatment of hydrogen atom in three-dimensional Quantum Mechanics.

Consider the Schrödinger equation for positronium in $`d`$ dimensions:

$$\left(\frac{p^2}{m}\frac{c(d)\alpha }{r^{d2}}\right)\mathrm{\Psi }=E\mathrm{\Psi },$$
(4)
where $`c(d)=\mathrm{\Gamma }(d/21)/\pi ^{d/21}`$ (see ). Let us rescale $`pp\gamma `$ and $`rr/\gamma `$ and choose $`\gamma =(m\alpha /2)^{\frac{1}{1+2ϵ}}`$. Then both the Coulomb Hamiltonian $`H=p^2/mc(d)\alpha /r^{d2}`$ and its eigenvalue $`E`$ scale as $`\gamma ^2/m`$. The wave function of a bound state, being normalized to unity, scales as $`\gamma ^{d/2}`$. Hence, when expressed in atomic units, the energy of a bound state and its properly normalized wave function depend on $`d`$ only. The scaling rules above provide sufficient information to write the matrix element of any operator in atomic units and therefore scale out all the dependence on the fine structure constant.

There is another important point that makes the extraction of the $`\mathrm{ln}\alpha `$ corrections possible using limited amount of information. In order to explain it, we remind the reader that in bound state calculations different contributions to the final result can be distinguished. In particular, there are so-called hard contributions. The $`𝒪(\alpha ^3)`$ hard corrections to the Ps decay rates are described by the three-loop Feynman diagrams for the process $`e^+e^{}2(3)\gamma `$, that have to be computed exactly at the threshold. Schematically, such diagrams generate the correction

$$V_{\mathrm{hard}}=\left(\frac{\alpha }{\pi }\right)^3\left[\frac{t_1}{ϵ^2}+\frac{t_2}{ϵ}+t_3\right]V_{\mathrm{Born}},$$
(5)
( $`t_{13}`$ are some constants) to the annihilation kernel $`V_{\mathrm{Born}}`$ that is responsible for the lowest order decay rate:

$$\mathrm{\Gamma }_{\mathrm{p},\mathrm{o}}^{(0)}=\mathrm{\Psi }|V_{\mathrm{Born}}|\mathrm{\Psi }\mathrm{\Psi }_0^2.$$
(6)
Here $`\mathrm{\Psi }_0`$ stands for the positronium wave function at the origin. The effective potential (5), in turn, generates the correction to the decay rate,

$$\delta _{\mathrm{hard}}\mathrm{\Gamma }=\mathrm{\Psi }|V_{\mathrm{hard}}|\mathrm{\Psi }\left[\frac{t_1}{ϵ^2}+\frac{t_2}{ϵ}+t_3\right]\mathrm{\Psi }_0^2.$$
(7)
If we rewrite $`\mathrm{\Psi }_0`$ in atomic units, Eq.(7) generates logarithms of the fine structure constant. These logarithms are artificial, since we anticipate that other, soft scale, contributions also generate divergences which exactly match and cancel all the divergences in $`V_{\mathrm{hard}}`$. This implies, that the logarithms associated with the rescaling of the wave function at the origin get cancelled as well. Therefore, the easiest way to avoid considering $`V_{\mathrm{hard}}`$ (which is not available at present) is to work with relative, rather than absolute, corrections to the decay width. This automatically discourages Eq.(7) as the source of logarithms of the fine structure constant since in this case there is simply nothing to rescale.

In order to illustrate how these arguments help to compute the $`\mathrm{ln}\alpha `$ corrections, let us consider the matrix element of a nonrelativistic operator $`𝒪`$ that delivers $`𝒪(\alpha ^3)`$ correction to the lowest order annihilation kernel $`V_{\mathrm{Born}}`$. In accordance with the above comment we consider relative correction to the decay rate:

$$\mathrm{\Delta }_𝒪=\frac{\delta _𝒪\mathrm{\Gamma }}{\mathrm{\Gamma }_{\mathrm{p},\mathrm{o}}^{(0)}}=\frac{\mathrm{\Psi }|𝒪|\mathrm{\Psi }}{\mathrm{\Psi }|V_{\mathrm{Born}}|\mathrm{\Psi }}.$$
(8)
The operator $`𝒪`$ is a function of coordinate and momentum operators that act on the positronium wave function. After rescaling of all the quantities on the right hand side of Eq.(8) as described above, we end up with the following equation:

$$\mathrm{\Delta }_𝒪=\alpha ^n\frac{\gamma ^{3n+lϵ}}{m^{3n+kϵ}}\frac{\mathrm{\Psi }|𝒪|\mathrm{\Psi }}{\mathrm{\Psi }|V_{\mathrm{Born}}|\mathrm{\Psi }}|_{\gamma =1},$$
(9)
where $`n=1,2`$ is a power of $`\alpha `$ that explicitly enters $`𝒪`$ and $`l,k`$ are some integer numbers. If the matrix element is finite, we can safely put $`ϵ=0`$ and then the relative correction to the decay width is $`\alpha ^3`$ times the $`\alpha `$-independent ratio and hence no logarithms of $`\alpha `$ appear. Therefore, after the rescaling, the only place where $`\mathrm{ln}\alpha `$ can come from is the expansion of the factor $`\gamma ^{3n+lϵ}`$ in powers of $`ϵ`$; this implies that in order to generate the $`\mathrm{ln}\alpha `$ corrections, the nonrelativistic matrix elements should diverge and only divergent pieces of the matrix elements have to be known to determine logarithmic in $`\alpha `$ corrections to the decay rate. Note also, that because of the relation between $`\gamma `$ and $`\alpha `$,

$$\gamma =\left(\frac{m\alpha }{2}\right)^{\frac{1}{1+2ϵ}},$$

even the operators that scale as integer powers of $`\gamma `$ can generate $`\mathrm{ln}\alpha `$ corrections.

Let us note that the extraction of the $`\mathrm{ln}\alpha `$ corrections in dimensional regularization can lead to some counter-intuitive results; for example, logarithms of the fine structure constant are generated by the operators, that in a more “physical” regularization schemes, such as e.g. the schemes that use either the photon mass or the momentum cut-off $`\lambda `$, can only lead to the logarithms of $`m/\lambda `$ but not to $`\mathrm{ln}\alpha `$. Therefore, it appears that individual contributions to the final result are scheme dependent. Nevertheless, we would like to stress that once dimensional regularization and clear rules for extracting $`\mathrm{ln}\alpha `$ are adopted, there is no other way as to consider all possible contributions; none of them can be disregarded by invoking the fact that in a different regularization scheme a particular operator cannot generate $`𝒪(\mathrm{ln}\alpha )`$ correction.

The calculation of the nonrelativistic contributions that are relevant at order $`𝒪(\alpha ^3\mathrm{ln}\alpha )`$ is described in the following Sections. Some useful integrals, that we need in the calculation, are summarized in Appendix. We give intermediate formulas for the relative logarithmic corrections to the positronium decay rate expressed in units of $`(\alpha ^3/\pi )\mathrm{\Gamma }(1+ϵ)^3(4\pi )^{3ϵ}`$.

## III Irreducible contribution

This particular contribution arises as the average value of a local operator that comes from the Taylor expansion of the one-loop corrections to the annihilation kernel in spatial momenta $`𝒑`$ of electron and positron. For our purpose we need $`𝒪(\alpha p^2)`$ correction to $`V_{\mathrm{Born}}`$. Since this operator is constructed by Taylor expanding in external momentum, the non-analytic dependence on $`p^2`$ cannot appear. Then, from the rescaling argument we know that only divergent piece of the Wilson coefficient of this operator is required. It turns out that this divergent piece can easily be computed using rather general arguments. According to the rules of nonrelativistic QED, we extract Wilson coefficients of various effective operators from the corresponding on-shell scattering amplitudes. The key observation is that the divergence in the Wilson coefficient of the $`𝒪(\alpha p^2)`$ operator causes the divergence in the on-shell annihilation process $`e^+e^{}n\gamma `$, $`n=2,3`$, and that it is in fact a “true” infrared divergence that should be compensated by the real emission of an additional soft photon in the same process.

In order to avoid confusion we stress that the mechanism of cancelling this divergence by real radiation does not apply to the bound state because of C-parity conservation. There is no contradiction, however. The real infrared divergences in the bound state calculations are absent because electron and positron in positronium are off shell. The divergence appears only when we put them on mass shell in order to extract the Wilson coefficient of the relevant operator. The crucial observation is that considering the process in a different kinematic regime (on-shell annihilation), we easily find a divergent piece of the appropriate Wilson coefficient from the known amplitude of the real soft photon emission.

Requiring that virtual and real corrections to the on-shell annihilation process $`e^+e^{}n\gamma `$ add up to a finite quantity, we get the following correction to the annihilation kernel $`V_{\mathrm{Born}}`$:

$$V_{\mathrm{irr}}=\frac{2\alpha }{3\pi ϵ}\frac{p^2+p^2}{m^2}V_{\mathrm{Born}}.$$
(10)
Applying the rescaling arguments, we end up with the following correction induced by the irreducible operator (10):

$$\mathrm{\Delta }_{\mathrm{irr}}=\frac{4}{3}\mathrm{ln}\alpha .$$
(11)
## IV “Hard loop” contributions

These contributions arise in the second order of the nonrelativistic perturbation theory. This means that the corresponding nonrelativistic operators $`𝒪`$ are of the form $`VGV^{}`$, where $`G`$ is the reduced Green function of the Coulomb Hamiltonian from Eq.(4) and $`V,V^{}`$ are some local operators with one of the two originating from a hard one-loop correction.

There are two sources of the “hard loop” contributions. The first one is the one-loop renormalization of the annihilation kernel:

$$V_{\mathrm{p},\mathrm{o}}=\frac{\alpha }{\pi }A_{\mathrm{p},\mathrm{o}}V_{\mathrm{Born}}.$$
(12)
The corresponding $`𝒪(\alpha ^3\mathrm{ln}\alpha )`$ correction to the decay rate is then easily related to the $`𝒪(\alpha ^2\mathrm{ln}\alpha )`$ correction computed in . We find:

$$\mathrm{\Delta }_{\mathrm{p},\mathrm{o}}=\left(\frac{7}{6}𝑺^22\right)A_{\mathrm{p},\mathrm{o}}\mathrm{ln}\alpha .$$
(13)
The second “hard loop” contribution corresponds to the “hard” piece in $`𝒪(m\alpha ^5)`$ effective potential. It reads :

$$V_{\mathrm{hl}}=\frac{\alpha ^2}{3m^2}\frac{\mathrm{\Gamma }(1+ϵ)}{(4\pi )^ϵ}\left[\frac{1}{ϵ}+\frac{39}{5}12\mathrm{ln}2+𝑺^2\left(\frac{32}{3}+6\mathrm{ln}2\right)2\mathrm{ln}m\right]\delta (𝒓),$$
(14)
where $`𝑺`$ is the operator of the total spin. There is no problem with defining total spin operator here, since it multiplies explicitly finite quantity.

The $`𝒪(\alpha ^3)`$ correction to the decay rate generated by the potential from Eq.(14) then reads:

$$\delta _{\mathrm{hl}}\mathrm{\Gamma }=2\mathrm{\Psi }|V_{\mathrm{Born}}GV_{\mathrm{hl}}|\mathrm{\Psi },$$
(15)
and is proportional to the Green function at the origin, $`G(0,0)`$. All necessary results for this Green function can be found in . Finally, we obtain:

$$\mathrm{\Delta }_{\mathrm{hl}}=\frac{\mathrm{ln}^2\alpha }{3}+\left(\frac{1}{6ϵ}+2\mathrm{ln}2\frac{59}{30}𝑺^2\left(\frac{16}{9}+\mathrm{ln}2\right)+\mathrm{ln}m\right)\mathrm{ln}\alpha .$$
(16)
Let us note that since there is an explicitly divergent term in Eq.(16), one may wonder whether or not the difference in the energy $`E`$ of the bound state in $`d`$ and three dimensions should be taken into account. We have checked that the cancellation of all divergent terms in the final result for the $`𝒪(\alpha ^3\mathrm{ln}\alpha )`$ correction to the decay rate occurs for arbitrary $`E`$ and for this reason we use three dimensional expression for this (rescaled) quantity to present individual contributions as well.

## V Seagull contribution

The correction to the decay rate,

$$\delta _\mathrm{s}\mathrm{\Gamma }=2\mathrm{\Psi }|V_{\mathrm{Born}}GV_\mathrm{s}|\mathrm{\Psi },$$
(17)
is induced by the double seagull effective potential $`V_\mathrm{s}`$,

$$V_\mathrm{s}(𝒒)=\frac{\pi ^2\alpha ^2}{m^2}\frac{d^dk}{(2\pi )^d}\frac{𝒫_{ij}(𝒌)𝒫_{ij}(𝒌^{})}{kk^{}(k+k^{})},$$
(18)
where $`𝒌^{}=𝒒𝒌`$ and $`𝒫_{ij}(𝒌)=\delta _{ij}k_ik_j/k^2`$.

To facilitate the calculation of this potential, we introduce auxiliary integration variable $`k_0`$ and rewrite $`V_\mathrm{s}(𝒒)`$ as follows:

$$V_\mathrm{s}(𝒒)=\frac{\alpha ^2}{m^2}_0^{\mathrm{}}\frac{dk_0}{2\pi }\frac{d^dk}{(2\pi )^d}T_{ij}(𝒌)T_{ij}(𝒒𝒌),$$
(19)
with

$$T_{ij}(𝒌)=\frac{4\pi 𝒫_{ij}(𝒌)}{k^2+k_0^2}.$$
(20)
It is convenient to consider Fourier transform of $`V_\mathrm{s}(𝒒)`$:

$$V_\mathrm{s}(𝒓)=\frac{\alpha ^2}{m^2}_0^{\mathrm{}}\frac{dk_0}{2\pi }T_{ij}(𝒓)T_{ij}(𝒓),$$
(21)
where $`T_{ij}(𝒓)`$ stands for

$$T_{ij}(𝒓)=\delta _{ij}\frac{d^dk}{(2\pi )^d}e^{i𝒌𝒓}\frac{4\pi }{k^2+k_0^2}+\frac{1}{k_0^2}_i_j\frac{d^dk}{(2\pi )^d}e^{i𝒌𝒓}\left(\frac{4\pi }{k^2}\frac{4\pi }{k^2+k_0^2}\right).$$
(22)
To compute $`T_{ij}(𝒓)`$ we use

$$\frac{d^dk}{(2\pi )^d}e^{i𝒌𝒓}\frac{4\pi }{k^2+k_0^2}=2\left(\frac{k_0}{2\pi r}\right)^{\frac{d}{2}1}K_{\frac{d}{2}1}(k_0r),$$
(23)
where $`K_\nu (x)`$ is the modified Bessel function of the second kind.

Inserting (23) to (22), and integrating over $`k_0`$ in Eq.(21), we obtain:

$$V_\mathrm{s}(𝒓)=\frac{\alpha ^2r^{3+4ϵ}}{2\pi m^2}\frac{\mathrm{\Gamma }(1ϵ)^2}{(4\pi )^{2ϵ}}\left[1ϵ\frac{178\mathrm{ln}2}{3}+𝒪(ϵ^2)\right].$$
(24)
We then rescale the relative correction to the decay rate and obtain:

$$\frac{\delta _\mathrm{s}\mathrm{\Gamma }}{\mathrm{\Gamma }_{\mathrm{p},\mathrm{o}}^{(0)}}=\frac{\alpha ^2\gamma ^{14ϵ}}{2\pi m}\frac{\mathrm{\Gamma }(1ϵ)^2}{(4\pi )^{2ϵ}}\left[1ϵ\frac{178\mathrm{ln}2}{3}+𝒪(ϵ^2)\right]d^drr^{3+4ϵ}G(𝒓,0)\frac{\mathrm{\Psi }(𝒓)}{\mathrm{\Psi }_0}.$$
(25)
As we explained previously, we need only a divergent part of the integral in Eq.(25). Since in three dimensions the Green function $`G(𝒓,0)`$ behaves as $`𝒪(r^1)`$ for small values of $`r`$, we may expand the wave function in series around $`r=0`$ and keep only two first terms in such an expansion:

$$\frac{\mathrm{\Psi }(𝒓)}{\mathrm{\Psi }_0}=1\frac{c(d)r^{4d}}{4d}+𝒪(r^{82d}).$$
(26)
The Green function $`G`$ can be written as a sum of three pieces $`G=G_0+G_1+G_{\mathrm{multi}}`$ according to the number of Coulomb interactions$`G_0`$ is the free Green function, $`G_1`$ is the single Coulomb correction to $`G_0`$, and $`G_{\mathrm{multi}}`$ accounts for two and more Coulomb interactions.. In three dimensions, $`G_0r^1`$, $`G_1\mathrm{ln}r`$ and $`G_{\mathrm{multi}}r^0`$ as $`r0`$. Therefore, the first term from the expansion Eq.(26) is sufficient to extract the singularities caused by the contributions of $`G_1`$ and $`G_{\mathrm{multi}}`$. For $`G_{\mathrm{multi}}`$ we derive:

$$d^drr^{3+4ϵ}G_{\mathrm{multi}}(𝒓,0)=4\pi _0^1𝑑rr^{1+2ϵ}G_{\mathrm{multi}}(0,0)+𝒪(1)=\frac{3}{ϵ}+𝒪(1).$$
(27)
In order to analyze $`G_{0,1}`$ contributions, it is convenient to switch to the momentum space. We obtain:

$`{\displaystyle d^drr^{3+4ϵ}G_1(𝒓,0)}={\displaystyle \frac{4^ϵ\pi ^{3/2ϵ}\mathrm{\Gamma }(ϵ)}{\mathrm{\Gamma }(3/22ϵ)}}{\displaystyle \frac{d^dp}{(2\pi )^d}\frac{G_1(𝒑)}{p^{2ϵ}}},`$ (28)
$`{\displaystyle d^drr^{3+4ϵ}G_0(𝒓,0)\left(1\frac{c(d)r^{1+2ϵ}}{1+2ϵ}\right)}={\displaystyle \frac{2\pi }{ϵ}}{\displaystyle \frac{d^dp}{(2\pi )^d}\frac{G_0(𝒑)}{p^{2ϵ}}\left(1\frac{\pi ϵ}{p^{1+2ϵ}}\right)}+𝒪(1).`$ (29)
The integrals in Eqs.(28,29) can be expressed through the integrals listed in Appendix:

$`{\displaystyle \frac{d^dp}{(2\pi )^d}\frac{G_1(𝒑)}{p^{2ϵ}}}`$ $`=`$ $`16\pi \left[I_2(1+ϵ,1,1)+I_1(1,1,1)I_2(1,1,1)\right]+𝒪(ϵ),`$ (30)
$`{\displaystyle \frac{d^dp}{(2\pi )^d}\frac{G_0(𝒑)}{p^{2z}}}`$ $`=`$ $`2I_0(1,z).`$ (31)
Substituting all the relevant expressions into Eq.(25) and expanding in $`ϵ`$, we arrive at the final result for the seagull contribution:

$$\mathrm{\Delta }_\mathrm{s}=\frac{3\mathrm{ln}^2\alpha }{2}\left(\frac{1}{2ϵ}+\frac{4\mathrm{ln}2}{3}+\frac{5}{3}3\mathrm{ln}m\right)\mathrm{ln}\alpha .$$
(32)
## VI Retardation

In this Section we discuss the retardation effect, caused by the exchange of a photon with a typical momentum of the order of the inverse Bohr radius $`km\alpha `$ between electron and positron in the bound state. The corresponding formula reads:

$`\delta _{\mathrm{ret}}\mathrm{\Gamma }`$ $`=`$ $`2\mathrm{\Psi }\left|V_{\mathrm{Born}}G{\displaystyle \frac{d^dk}{(2\pi )^d}j_i^ee^{i𝒌𝒓_e}\frac{4\pi \alpha }{2k}\frac{𝒫_{ij}(𝒌)}{k+HE}j_j^pe^{i𝒌𝒓_p}}\right|\mathrm{\Psi }`$ (34)
$`\mathrm{\Psi }\left|V_{\mathrm{Born}}\right|\mathrm{\Psi }\mathrm{\Psi }\left|{\displaystyle \frac{d^dk}{(2\pi )^d}j_i^ee^{i𝒌𝒓_e}\frac{4\pi \alpha }{2k}\frac{𝒫_{ij}(𝒌)}{(k+HE)^2}j_j^pe^{i𝒌𝒓_p}}\right|\mathrm{\Psi }+\mathrm{H}.\mathrm{c}.,`$
where $`𝒋^e=𝒑^e/m[𝝈𝒌,𝝈]/(4m)`$ and $`𝒋^p=𝒑^p/m+[𝝈^{}𝒌,𝝈^{}]/(4m)`$ are electron and positron currents, respectively.

The appearance of two matrix elements in Eq.(34) is related to the fact that we need the second order correction to the wave function of a bound state. In contrast to the first order, at second order of perturbation theory one should be careful to maintain the normalization of the wave function. This is the reason why the second term in Eq.(34) appears. The general formula for the second order correction to the wave function can be found in .

When writing Eq.(34), we have assumed that a magnetic photon with the momentum $`𝒌`$ is emitted by the electron at the point $`𝒓_e`$ and absorbed by the positron at the point $`𝒓_p`$. Using the fact that $`HEm\alpha ^2`$ and $`km\alpha `$ for the retardation effects, we expand the r.h.s. of Eq.(34) in powers of $`(HE)/k`$ and obtain:

$$\delta _{\mathrm{ret}}\mathrm{\Gamma }=2\mathrm{\Psi }|V_{\mathrm{Born}}GV_{\mathrm{ret}}|\mathrm{\Psi }\mathrm{\Psi }|V_{\mathrm{Born}}|\mathrm{\Psi }\mathrm{\Psi }|V_{\mathrm{ret}}^{}|\mathrm{\Psi },$$
(35)
where the “retardation potential” $`V_{\mathrm{ret}}`$ reads:

$$V_{\mathrm{ret}}=\frac{d^dk}{(2\pi )^d}\frac{4\pi \alpha 𝒫_{ij}(𝒌)}{2k^3}j_i^ee^{i𝒌𝒓_e}(HE)j_j^pe^{i𝒌𝒓_p}+\mathrm{H}.\mathrm{c}.,$$
(36)
and $`V_{\mathrm{ret}}^{}=V_{\mathrm{ret}}/E`$. Let us illustrate how one deals with such expressions using the spin-dependent parts of the currents as an example. The corresponding contribution to $`V_{\mathrm{ret}}`$ then reads:

$$V_{\mathrm{ret}}^{\mathrm{spin}}=\frac{\pi \alpha }{16m^2}\frac{d^dk}{(2\pi )^d}\frac{[𝝈𝒌,\sigma _i][𝝈^{}𝒌,\sigma _i^{}]}{k^3}\left\{(HE)e^{i𝒌𝒓}+e^{i𝒌𝒓}(HE)\right\}+\mathrm{H}.\mathrm{c}.,$$
(37)
where the relative coordinate $`𝒓=𝒓_e𝒓_p`$ has been introduced. It is now easy to see that if we insert this result into Eq.(35) and use the Schrödinger equation both for the wave function $`\mathrm{\Psi }`$ and for the reduced Green function $`G`$, we obtain a finite correction to the decay rate. As we explained in the Introduction, finite contributions cannot generate logarithms of the fine structure constant.

The analysis of the spin-independent contribution to $`V_{\mathrm{ret}}`$ is only slightly more cumbersome. We again pull out $`HE`$ both to the right and to the left by commuting it with either electron or positron current and obtain the following expression for the retardation potential:

$$V_{\mathrm{ret}}=\frac{\pi \alpha }{m^2}\frac{d^dk}{(2\pi )^d}\frac{𝒫_{ij}(𝒌)}{k^3}\left(\{HE,e^{i𝒌𝒓}p_ip_j\}+[p_i,[H,p_j]]e^{i𝒌𝒓}\right)+\mathrm{H}.\mathrm{c}.,$$

where $`\{HE,e^{i𝒌𝒓}p_ip_j\}`$ denotes the anticommutator of the two operators. Using this expression for the retardation potential in Eq.(35) and applying equations of motion, we arrive at the following correction to the decay rate:

$$\delta _{\mathrm{ret}}\mathrm{\Gamma }=\frac{4\pi \alpha }{m^2}\mathrm{\Psi }|V_{\mathrm{Born}}G[U_{ij}[C,p_i],p_j]|\mathrm{\Psi },$$
(38)
where $`C=c(d)\alpha /r^{d2}`$ is the Coulomb potential, and $`U_{ij}(𝒓)`$ is defined as

$$U_{ij}(𝒓)=\frac{d^dk}{(2\pi )^d}\frac{𝒫_{ij}(𝒌)}{k^3}e^{i𝒌𝒓}=\frac{\mathrm{\Gamma }(ϵ)r^{2ϵ}}{6\pi ^{2ϵ}}\left(\delta _{ij}ϵn_in_j\right).$$
(39)
Performing rescaling and computing the double commutator in Eq.(38), we get the result for the relative correction to the decay rate,

$$\frac{\delta _{\mathrm{ret}}\mathrm{\Gamma }}{\mathrm{\Gamma }_{\mathrm{p},\mathrm{o}}^{(0)}}=\frac{8\alpha \gamma ^{22ϵ}}{3m^2}\frac{\mathrm{\Gamma }(2ϵ)\mathrm{\Gamma }(3/2ϵ)}{\pi ^{3/2ϵ}}d^drr^{3+4ϵ}G(𝒓,0)\frac{\mathrm{\Psi }(𝒓)}{\mathrm{\Psi }_0},$$
(40)
which is very similar to the correction induced by the seagull potential, Eq.(25). We can, therefore, borrow much of the analysis from the previous Section. We finally obtain:

$$\mathrm{\Delta }_{\mathrm{ret}}=2\mathrm{ln}^2\alpha \left(\frac{2}{3ϵ}\frac{4\mathrm{ln}2}{3}+44\mathrm{ln}m\right)\mathrm{ln}\alpha .$$
(41)
## VII Ultrasoft contribution

By definition, the ultrasoft contribution is due to the photons with energy and momentum of order $`m\alpha ^2`$. Such soft photons cannot resolve the structure of the bound state and for this reason they interact directly with the positronium as a whole. Since the positronium is chargeless, the interaction is necessarily of the dipole nature.

A general formula for the correction to the decay rate caused by ultrasoft contributions is:

$`\delta _{\mathrm{us}}\mathrm{\Gamma }`$ $`=`$ $`2\mathrm{\Psi }\left|V_{\mathrm{Born}}G{\displaystyle \frac{2p_i}{m}}{\displaystyle \frac{d^dk}{(2\pi )^d}\frac{4\pi \alpha }{2k}\frac{𝒫_{ij}(𝒌)}{EHk}\frac{2p_j}{m}}\right|\mathrm{\Psi }`$ (43)
$`\mathrm{\Psi }\left|V_{\mathrm{Born}}\right|\mathrm{\Psi }\mathrm{\Psi }\left|{\displaystyle \frac{2p_i}{m}}{\displaystyle \frac{d^dk}{(2\pi )^d}\frac{4\pi \alpha }{2k}\frac{𝒫_{ij}(𝒌)}{(EHk)^2}\frac{2p_j}{m}}\right|\mathrm{\Psi }.`$
After integration over directions of $`𝒌`$ and performing the rescaling we get:

$$\frac{\delta _{\mathrm{us}}\mathrm{\Gamma }}{\mathrm{\Gamma }_{\mathrm{p},\mathrm{o}}^{(0)}}=\frac{4^{2ϵ}\pi \alpha \gamma ^{24ϵ}}{m^{22ϵ}}\frac{\mathrm{\Omega }_d}{(2\pi )^d}\frac{d1}{d}\left(0\left|G𝒑_0^{\mathrm{}}\frac{dkk^{2ϵ}}{k+HE}[H,𝒑]\right|\stackrel{~}{\mathrm{\Psi }}\mathrm{\Psi }\left|𝒑_0^{\mathrm{}}\frac{dkk^{12ϵ}}{(k+HE)^2}𝒑\right|\mathrm{\Psi }\right).$$
(44)
Both matrix elements in Eq.(44) must be computed in atomic units; also $`\stackrel{~}{\mathrm{\Psi }}=\mathrm{\Psi }/\mathrm{\Psi }_0`$, and $`\mathrm{\Omega }_d=2\pi ^{3/2ϵ}/\mathrm{\Gamma }(3/2ϵ)`$ is the $`d`$-dimensional angular volume.

The second matrix element in Eq.(44) is easy to computeRecall that only the terms that are singular for $`ϵ0`$ are needed.:

$$\mathrm{\Psi }\left|𝒑_0^{\mathrm{}}\frac{dkk^{12ϵ}}{(k+HE)^2}𝒑\right|\mathrm{\Psi }=\frac{\mathrm{\Psi }\left|𝒑(HE)^{2ϵ}𝒑\right|\mathrm{\Psi }}{2ϵ}=\frac{1}{2ϵ}.$$
(45)
Then, consider the first matrix element from Eq.(44). It is convenient to write it as a sum:

$$0\left|G𝒑_0^{\mathrm{}}\frac{dkk^{2ϵ}}{k+HE}[H,𝒑]\right|\stackrel{~}{\mathrm{\Psi }}=M_0+M_1,$$
(46)
where the two terms correspond to the number of Coulomb interactions between the moments of emission and absorption of the ultrasoft photon§§§If two or more Coulomb photons are exchanged, the resulting matrix element becomes finite and, in accordance with the argument given in Section II, it cannot generate $`𝒪(\mathrm{ln}\alpha )`$ correction..

Let us consider $`M_0`$. In this case, integrating over $`k`$ we obtain:

$$M_0=\frac{1}{2ϵ}0\left|G\left(H_0E\right)^{2ϵ}𝒑[H,𝒑]\right|\stackrel{~}{\mathrm{\Psi }},$$
(47)
where $`H_0=p^2/2`$ is the free Hamiltonian in atomic units. Using

$$𝒑[H,𝒑]=[𝒑,[H,𝒑]]+[H,𝒑]𝒑=4\pi \delta (𝒓)+[H,𝒑]𝒑,$$
(48)
we represent $`M_0`$ as the sum of two terms,

$$M_0=M_\delta +M_\psi ,$$
(49)
where

$`M_\delta `$ $`=`$ $`{\displaystyle \frac{4\pi }{2ϵ}}0\left|G\left(H_0E\right)^{2ϵ}\right|0,`$ (50)
$`M_\psi `$ $`=`$ $`{\displaystyle \frac{1}{2ϵ}}0\left|G\left(H_0E\right)^{2ϵ}[H,𝒑]𝒑\right|\stackrel{~}{\mathrm{\Psi }}.`$ (51)
We now analyze the two terms in Eq.(49) separately. In order to compute $`M_\delta `$ we use the fact that only two first terms in the expansion of the Green function in the number of Coulomb exchanges, $`G_0`$ and $`G_1`$, diverge at zero spatial separation. Since $`G_{\mathrm{multi}}(0,0)`$ is finite, we can compute it for $`d=3`$. On the other hand, both $`G_0`$ and $`G_1`$ are known explicitly (see e.g. ) and hence the corresponding integrals can easily be computed. Using expressions for integrals summarized in Appendix, we arrive at the following expression for $`M_\delta `$:

$$M_\delta =\frac{2^{2+2ϵ}\pi }{ϵ}I_0(1+2ϵ,0)\frac{2^{5+2ϵ}\pi ^2}{ϵ}I_1(1+2ϵ,1,1)\frac{3}{ϵ}.$$
(52)
Consider the second matrix element $`M_\psi `$. In this case, it is sufficient to separate $`G=G_0+(GG_0)`$. Since there is an overall factor $`ϵ^1`$ in Eq.(51), the finite matrix element that contains $`GG_0`$ can be calculated in three dimensions. Using the expression for $`GG_0`$ from , we derive:

$$M_\psi ^\mathrm{f}=\frac{1}{2ϵ}0\left|(GG_0)[H,𝒑]𝒑\right|\stackrel{~}{\mathrm{\Psi }}=\frac{2}{ϵ}.$$
(53)
The matrix element containing $`G_0`$ can be rewritten in momentum space as

$$M_\psi ^\mathrm{i}=\frac{1}{2ϵ}0\left|G_0\left(H_0E\right)^{2ϵ}[H,𝒑]𝒑\right|\stackrel{~}{\mathrm{\Psi }}=\frac{1}{2ϵ}\frac{d^dp^{}d^dp}{(2\pi )^{2d}}\left(\frac{2}{p^2+1}\right)^{1+2ϵ}\frac{4\pi (𝒑^{}𝒑)𝒑}{(𝒑^{}𝒑)^2}\frac{\mathrm{\Psi }(𝒑)}{\mathrm{\Psi }_0},$$
(54)
where $`\mathrm{\Psi }(𝒑)`$ is the Fourier transform of $`\mathrm{\Psi }(𝒓)`$. Since $`M_\psi ^\mathrm{i}`$ has an overall divergence and hence we need the integral up to a constant only, we can use the three dimensional expression for the wave function $`\mathrm{\Psi }(𝒑)`$:

$$M_\psi ^\mathrm{i}=\frac{4\pi ^2}{ϵ}\frac{d^dp^{}d^dp}{(2\pi )^{2d}}\left(\frac{2}{p^2+1}\right)^{1+2ϵ}\frac{(𝒑^{}𝒑)𝒑}{(𝒑^{}𝒑)^2}\left(\frac{2}{p^2+1}\right)^2.$$
(55)
We now rewrite the scalar product in the numerator as a linear combination of denominators and obtain our final expression for $`M_\psi =M_\psi ^\mathrm{f}+M_\psi ^\mathrm{i}`$:

$$M_\psi =\frac{4^{2+ϵ}\pi ^2}{ϵ}\left[I_1(2ϵ,1,2)I_1(1+2ϵ,1,1)I_1(1+2ϵ,0,2)\right]+\frac{2}{ϵ}.$$
(56)
Consider next the matrix element $`M_1`$. We start with the following expression:

$$M_1=0\left|G𝒑_0^{\mathrm{}}𝑑kk^{2ϵ}\frac{1}{k+H_0E}C\frac{1}{k+H_0E}[H,𝒑]\right|\stackrel{~}{\mathrm{\Psi }}.$$
(57)
Simple power counting shows that we can safely take $`\stackrel{~}{\mathrm{\Psi }}`$ at the origin, $`\stackrel{~}{\mathrm{\Psi }}1`$, and also replace the Green function $`G`$ by its high-momentum asymptotics, $`G(p)2/p^2`$. We then obtain symmetric and uniform expression

$$M_1=8\pi _0^{\mathrm{}}𝑑kk^{2ϵ}\frac{d^dp^{}d^dp}{(2\pi )^{2d}}\frac{2}{p^2(p^2+1+2k)}\frac{4\pi 𝒑^{}𝒑}{(𝒑^{}𝒑)^2}\frac{2}{p^2(p^2+1+2k)}.$$
(58)
If we now rescale both $`𝒑`$ and $`𝒑^{}`$ as $`𝒑\sqrt{2k+1}𝒑`$, the integration over $`k`$ factorizes and we obtain:

$$M_1=\frac{8\pi ^2}{ϵ}\left[2I_2(1,1,1)2I_1(1,1,1)I_0(1,1)^2\right].$$
(59)
Finally, using explicit expressions for the integrals from Appendix, we arrive at our final result for the ultrasoft correction to the decay rate:

$$\mathrm{\Delta }_{\mathrm{us}}=\frac{16\mathrm{ln}\alpha ^2}{3}+\left(\frac{4}{3ϵ}+8\mathrm{ln}2+\frac{20}{9}8\mathrm{ln}m\right)\mathrm{ln}\alpha .$$
(60)
## VIII Conclusions

The sum of all the contributions from Eqs.(11,13,16,32,41,60) gives the final result for the $`𝒪(\alpha ^3\mathrm{ln}\alpha )`$ corrections to Ps decay rate:

$`{\displaystyle \frac{\mathrm{\Delta }\mathrm{\Gamma }_\mathrm{p}}{\mathrm{\Gamma }_\mathrm{p}^{(0)}}}={\displaystyle \frac{\alpha ^3}{\pi }}\left[{\displaystyle \frac{3}{2}}\mathrm{ln}^2\alpha +\mathrm{ln}\alpha \left\{{\displaystyle \frac{367}{90}}+10\mathrm{ln}22A_\mathrm{p}\right\}\right]={\displaystyle \frac{\alpha ^3}{\pi }}\left[{\displaystyle \frac{3}{2}}\mathrm{ln}^2\alpha +7.919\mathrm{ln}\alpha \right],`$ (61)
$`{\displaystyle \frac{\mathrm{\Delta }\mathrm{\Gamma }_\mathrm{o}}{\mathrm{\Gamma }_\mathrm{o}^{(0)}}}={\displaystyle \frac{\alpha ^3}{\pi }}\left[{\displaystyle \frac{3}{2}}\mathrm{ln}^2\alpha +\mathrm{ln}\alpha \left\{{\displaystyle \frac{229}{30}}+8\mathrm{ln}2+{\displaystyle \frac{A_\mathrm{o}}{3}}\right\}\right]={\displaystyle \frac{\alpha ^3}{\pi }}\left[{\displaystyle \frac{3}{2}}\mathrm{ln}^2\alpha 5.517\mathrm{ln}\alpha \right].`$ (62)
Numerically, these corrections cause a negligible change in the theoretical prediction for p-Ps and o-Ps lifetimes at the current level of precision. It is interesting to note, however, that the magnitude of the leading $`𝒪(\alpha ^3\mathrm{ln}^2\alpha )`$ and the subleading $`𝒪(\alpha ^3\mathrm{ln}\alpha )`$ corrections is comparable; in case of o-Ps they almost cancel each other.

Our results Eqs.(61,62) are in agreement with two recent calculations of $`𝒪(\alpha ^3\mathrm{ln}\alpha )`$ corrections . In Ref. the result for $`𝒪(\alpha ^3\mathrm{ln}\alpha )`$ correction to o-Ps decay rate has been obtained numerically, where as in Ref. analytical methods similar to ours have been employed. We believe that the achieved agreement between three independent calculations ensures that the results, Eqs.(61,62), are correct.

As we mentioned, the $`𝒪(\alpha ^3\mathrm{ln}\alpha )`$ correction to Ps decay rates at present is not very interesting phenomenologically. A more important question, which we think we fully addressed in this paper, is how the logarithms of the fine structure constant can be efficiently extracted in the bound state calculation when the dimensional regularization is used to regulate the nonrelativistic dynamics. It is true that dimensional regularization offers many technical advantages in the calculation. This does not go without a price, however, since one has to be extremely careful in defining basic objects of the nonrelativistic theory, e.g. the wave functions and energies. If this is not done, one is left guessing whether or not is the calculation correct.

Our key observation, which we think cures such problems and makes our calculation unambiguous, is the fact that the matrix elements in $`d`$ dimensions are the uniform functions of the fine structure constant, and that the corresponding power of $`\alpha `$ can be determined by expressing the matrix elements in “$`d`$-dimensional” atomic units. We think that these arguments have not been spelled out before in the literature on one hand, and that they are necessary to make a convincing case, on the other.

Finally, let us note that the technique discussed in this paper can obviously be used in other bound state QED problems, as well as for the heavy quarkonium states in QCD.

## IX Acknowledgments

This work was supported in part by DOE under grant number DE-AC03-76SF00515, and by the Russian Foundation for Basic Research grant 00-02-17646. We are grateful to the BNL High Energy Theory Group, where part of this work has been done, for the hospitality extended to us.

## Appendix

Definition of the integrals that were used in the derivation.

$$I_0(a,b)=\frac{d^dp}{(2\pi )^d}\left(\frac{1}{p^2+1}\right)^a\left(\frac{1}{p^2}\right)^b=(4\pi )^{d/2}\frac{\mathrm{\Gamma }(a+bd/2)\mathrm{\Gamma }(d/2b)}{\mathrm{\Gamma }(a)\mathrm{\Gamma }(d/2)},$$
(63)
$`I_1(a,b,c)`$ $`={\displaystyle \frac{d^dp^{}d^dp}{(2\pi )^{2d}}\left(\frac{1}{p^2+1}\right)^a\left(\frac{1}{(𝒑^{}𝒑)^2}\right)^b\left(\frac{1}{p^2+1}\right)^c}=(4\pi )^d`$ (65)
$`\times {\displaystyle \frac{\mathrm{\Gamma }(a+b+cd)\mathrm{\Gamma }(a+bd/2)\mathrm{\Gamma }(b+cd/2)\mathrm{\Gamma }(d/2b)}{\mathrm{\Gamma }(a)\mathrm{\Gamma }(c)\mathrm{\Gamma }(d/2)\mathrm{\Gamma }(a+2b+cd)}},`$
$`I_2(a,b,c)`$ $`={\displaystyle \frac{d^dp^{}d^dp}{(2\pi )^{2d}}\left(\frac{1}{p^2}\right)^a\left(\frac{1}{(𝒑^{}𝒑)^2}\right)^b\left(\frac{1}{p^2+1}\right)^c}=(4\pi )^d`$ (67)
$`\times {\displaystyle \frac{\mathrm{\Gamma }(a+b+cd)\mathrm{\Gamma }(a+bd/2)\mathrm{\Gamma }(d/2a)\mathrm{\Gamma }(d/2b)}{\mathrm{\Gamma }(a)\mathrm{\Gamma }(b)\mathrm{\Gamma }(c)\mathrm{\Gamma }(d/2)}}`$


format: text

text:
# LARGE ORBITAL MAGNETIC MOMENT IN FeBr2

## Abstract

Magnetic moment of the Fe2+ ion in FeBr2 has been calculated to be 4.32 $`\mu _B`$ in the magnetically-ordered state at 0 K. It is composed from the spin moment of 3.52 $`\mu _B`$ and the orbital moment of +0.80 $`\mu _B.`$ These calculations show that for the meaningful analysis of e-m properties of FeBr2 the spin-orbit coupling is essentially important and that the orbital moment is largely unquenched.

PACS 75.10.Dg 71.70.Ch 71.28+d 75.30.Mb

Keywords: crystal-field interactions, spin-orbit coupling, Fe2+ ion, FeBr2

FeBr2 exhibits the strong metamagnetic properties . At 4.2 K an external field of 3.15 T causes a jump of the magnetization from the almost zero value to a very big value of 110 emu/g . This latter value corresponds to 4.4 $`\mu _B`$ per the iron-ion moment. This value exceeds a theoretical value of 4.0 $`\mu _B`$ expected for the spin-only moment. The explanation of this large magnetic moment is a subject of long-time debate within the 3d-magnetism theoretisians \[3-9\].

We have calculated the atomic-like magnetic moment of the Fe2+ ion in the ionic FeBr2 compound to be 4.32 $`\mu _B`$ in the magnetically-ordered state at 0 K. This moment is composed from the spin moment of 3.52 $`\mu _B`$ and the orbital moment of +0.80 $`\mu _B`$. In the paramagnetic state the total moment amounts to 4.26 $`\mu _B`$ (3.48 and 0.78 $`\mu _B`$).

The calculations have been performed within the quantum atomistic solid-state theory, that points out the existence of the discrete energy spectrum in the 3d-ion containing compounds associated with quasi-atomic 3d states. We have considered 6 d electrons of the Fe2+ ion to form the highly-correlated atomic like electronic system 3d6. Such the system is characterized by L=2 and S=2 \[10-12\]. We have taken into account the crystal field (CEF) interactions (cubic CEF parameter B4 = +200 K, off-cubic trigonal distortion parameter B$`{}_{}{}^{0}{}_{2}{}^{}`$= -30 K) and the intra-atomic spin-orbit coupling ($`\lambda _{so}`$ = -150 K).

The crystal-field interactions in the hexagonal structure of FeBr2 cause the Fe-ion moment to lie along the hexagonal axis. The hexagonal axis is the trigonal axis of the local Br1- octahedron. As a consequence the antiferromagnet FeBr2 is very good example of the Ising-type magnetic compound. The metamagnetic transition is associated with the reversal of the AF moment. The Fe2+-ion moment strongly depends on the local symmetry. The local symmetry is reflected in the symmetry of the CEF. The low-temperature magnetic moment can vary from the zero value to almost 6 $`\mu _B`$ as the function of the local surrounding of the Fe2+ ion. The zero value is the manifestation of the Kramers and Jahn-Teller theorems for the highly-correlated electronic system with the even number of involved electrons.

In conclusion, the Fe2+ ion in FeBr2 exhibits large atomic-like magnetic moment with the very substantial orbital magnetic moment.The derived electronic structure exhibits the discrete energy spectrum, like in Refs 6 and 13 (the Co3+ ion is equivalent to the Fe2+ ion), in contrary to the presently-in-fashion results with the continuous energy spectrum for the 3d states \[7-9\].

Figure captions:Fig.1. The calculated temperature dependence of the local magnetic moment of the Fe2+ ion in FeBr2 together with the spin and orbital contributions.


text:
# Science User Scenarios for a VO Design Reference Mission

## 1 Science Requirements for Data Mining

One of the major functions of a Virtual Observatory (VO) is to facilitate data mining and knowledge discovery within the very large astronomical databases that are now coming on-line (or soon will be). A similarly important function of the VO is to facilitate linkages and cross-archive investigations utilizing these new data in conjunction with the rich legacy data archives that preceded them. The scientific teams that generate large (multi-Terabyte) databases cannot begin to tap their full scientific potential. Thus a significant portion of the astronomical research community and a comprehensive suite of research tools should be brought to bear on extracting the maximum scientific return for the huge investment in these large astronomical facilities, large surveys, and large scientific data systems. One approach to this problem can be identified as “data mining”.

What is data mining and why is applicable to scientific research? Data mining is defined as an information extraction activity whose goal is to discover hidden facts contained in databases. Data mining has taken the business community by storm and the phrase has become a bit overworked to describe some fairly routine functions of marketing. Even so, there are consequently now a vast array of resources and research techniques available for exploitation by the scientific communities. It is useful therefore to examine a further categorization of data mining thrusts and their sub-components, since these are likewise applicable to the scientific exploration of large astronomical databases.

In the marketing community, data mining is used to find patterns and relationships in data by using sophisticated techniques to build models — abstract representations of reality. A good model is a useful guide to understanding that reality and to making decisions. There are two main types of data mining models: descriptive and predictive. Descriptive models describe patterns in data and are generally used to create meaningful subgroups or clusters. Predictive models are used to forecast explicit values, based upon patterns determined from known results. These models are applicable to scientific inquiry as well.

There is another differentiation of data mining into two categories that we find particularly appropriate to knowledge discovery in large astronomical databases: event-based mining and relationship-based mining. At the risk of trivializing some fairly sophisticated techniques, we classify event-based mining scenarios into four orthogonal categories:

* Known events / known algorithms — use existing physical models (descriptive models) to locate known phenomena of interest either spatially or temporally within a large database.
* Known events / unknown algorithms — use pattern recognition and clustering properties of data to discover new observational (in our case, astrophysical) relationships among known phenomena.
* Unknown events / known algorithms — use expected physical relationships (predictive models) among observational parameters of astrophysical phenomena to predict the presence of previously unseen events within a large complex database.
* Unknown events / unknown algorithms — use thresholds to identify transient or otherwise unique (“one-of-a-kind”) events and therefore to discover new phenomena.

Similarly, for relationship-based mining, we identify three classes of association-driven scenarios that would find application in astronomical research:

* Spatial associations — identify events (astronomical objects) at the same location in the sky.
* Temporal associations — identify events occurring during the same or related periods of time.
* Coincidence associations — use clustering techniques to identify events that are reasonably co-located within a multi-dimensional parameter space.

Therefore, from this discussion, we can derive a reduced set of science requirements for data mining. These requirements correspond to the following set of exploratory approaches to mining large databases : Object Cross-Identification, Object Cross-Correlation, Nearest-Neighbor Identification, and Systematic Data Exploration. (a) “Object cross-identification” refers to the classical problem of connecting the source list in one catalog (or observation database) to the source list in another, in order to derive new astrophysical understanding of the cross-identified objects (e.g., gamma-ray burst counterparts). (b) “Object cross-correlation” refers to the application of “what if” scenarios to the full suite of parameters in a database (e.g., identify distant galaxies as $`U`$-band dropouts in a color-color scatter plot from the HDF survey). (c)“Nearest-neighbor identification” refers to the general application of clustering algorithms in multi-dimensional parameter space (e.g., finding the closest known population of young stars – in the TW Hydrae association – through their similar kinematics, X-ray emission, $`H\alpha `$, and Li abundance). (d) “Systematic data exploration” refers to the application of the broad range of event-based and relationship-based queries to a database in the hope of making a serendipitous discovery of new objects or a new class of objects (e.g., finding new types of variable stars, such as “bumpers”, in the MACHO database).

## 2 User Scenario #1: Estimating the Galaxy Interaction Rate

It is well established that a significant fraction of all galaxies have been involved in a galaxy-galaxy interaction and perhaps a merger at some time(s) in their past. The rate of these interactions is not yet well determined empirically: either the current rate for galaxies in the nearby Universe, or the cosmologically evolving rate in the distant Universe. Numerical simulations of the galaxy population and of the evolving hierarchical structure within various cosmological scenarios give a handle on the interaction and merger rates, which naturally depend on the choice of cosmological model. In general, the simulations confirm the importance and relatively high frequency of occurrence of interactions and mergers. Given the cosmological significance of interactions to galaxy formation and evolution, it is important to derive a firm value for the galaxy interaction rate observationally, for comparison with the numerical models, which will in turn help to narrow the plausible range of cosmological models, galaxy formation models, and galaxy evolution models.

We attempted an initial exploration of several on-line databases in order to estimate the galaxy interaction rate. We began by exploring an on-line catalog of galaxies (available through NASA’s ADC = Astronomical Data Center): the Updated Zwicky Catalog of Falco et al. (1999). This catalog identifies multiple-galaxy groupings, which we used to reduce the full list of 19,000 galaxies to the set of 1800 multiples. We then selected a very small sub-sample from this list to conduct a proof-of-concept investigation. We used existing catalog visualization tools and archive linkage tools at the ADC to find all possible NASA mission data and most of the all-sky survey data for these selected objects (Kargatis et al. 1999). We then identified characteristics in the optical images or in the IRAS fluxes or in the X-ray emissions to verify that the associated multiple galaxy systems are in fact (to high probability) bound groups (pairs, triples, quartets, etc.). The expectation that these small galaxy-galaxy separations and other evidences for physical association do in fact imply an on-going interaction was often confirmed through inspection of the DSS (Digital Sky Survey) imagery, which showed signs of interaction in many cases (e.g., distorted morphologies). Thus, by applying knowledge of astrophysical signatures of interactions, we were able to explore multiple databases (ADC catalogs, NASA mission archives, and ground-based sky surveys) in a coherent organized manner. We estimate that the galaxy interaction rate in the local Universe is approximately 8%.

## 3 User Scenario #2: In Search of the CIB

Among the exciting results of the COBE mission was the discovery by the DIRBE team of an extragalactic CIB (Cosmic Infrared Background; Hauser et al. 1988). There has been a storm of activity in the research community to identify the sources of the CIB and to understand their power sources (i.e., what powers the strong IR emissions? are they dust-enshrouded quasars? or dusty starbursts? or both?). Some possible counterparts to the CIB include: (a) ultraluminous IR galaxies (ULIRGs; see Sanders & Mirabel 1996 for a review); (b) SCUBA submm sources (Barger et al. 1999; Blain et al. 1999); (c) IR-selected AGN (a new population of AGN identified through the 2MASS survey; Beichman et al. 1998); or (d) Extremely Red Objects (EROs), which have been reported in several recent deep surveys (Smail et al. 1999; Thompson et al. 1999).

We initiated a proof-of-concept search scenario for identifying potential candidate contributors to the CIB. Our approach is similar to that of Haasrma & Partridge (1998), except that we are applying the full power of on-line databases and linkages between these databases, archives, and published literature. Our search scenario involved finding object cross-identifications among the IRAS Faint Source Catalog and FIRST survey catalog, and then attempting to find those commonly identified objects also within other databases, such as the HST observation log. In a very limited sample of targets that we investigated to test our “ADC as a mini-NVO” approach to the problem, we did find one object in common among the HST-IRAS-FIRST databases: a known hyperluminous infrared galaxy (HyLIRG) at $`z`$=0.780 harboring an AGN, which was specifically imaged by HST because of its known HyLIRG characteristics. In this extremely limited test scenario, we did in fact find what we were searching for: a distant IR-luminous galaxy that is a likely contributor to the CIB, or else similar in characteristics to the more distant objects that likely comprise the CIB.

The preliminary results of the investigations described above for User Scenarios #1 and #2 (including screen shots of the user interfaces employed in the studies) are presented at the following web site (under the category “How to use the ADC for scientific research projects”):

http://adc.gsfc.nasa.gov/adc/how\_to.html


text:
# 1 Muon-neutrino charged-current events per kiloton year, assuming no oscillation, for a 4–MW upgraded FNAL wide-band conventional source, and for a 20-GeV muon storage ring giving 1020 decays per year. The numbers are scaled from Ref. [1].

SLAC–PUB–8587
August 2000

Conventional Beams or Neutrino
Factories: The Next Generation of
Accelerator-Based Neutrino Experiments ***Work supported by Department of Energy contract DE–AC03–76SF00515.

Burton Richter

Stanford Linear Accelerator Center

Stanford, California 94309-0450

Abstract

The purpose of this paper is to provoke a discussion about the right next step in accelerator-based neutrino physics. In the next five years many experiments will be done to determine the neutrino mixing parameters. However, the small parameters $`\theta _{13},\mathrm{\Delta }m_{21}^2`$, and the CP violating phase are unlikely to be well determined. Here, I look at the potential of high-intensity, low-energy, narrow-band conventional neutrino beams to determine these parameters. I find, after roughly estimating the possible intensity and purity of conventional neutrino and anti-neutrino beam, that $`\mathrm{sin}^2\theta _{13}`$ can be measured if greater than a few parts in ten thousand, $`\mathrm{\Delta }m_{21}^2`$ can be measured if it is greater than $`4\times 10^5`$ (eV)2, and the CP violating phase can be measured if it is greater than 20 and the other parameters are not at their lower bounds. If these conclusions stand up to more detailed analysis, these experiments can be done long before a muon storage ring source could be built, and at much less cost.

Neutrino physics has become one of the central studies of current high-energy physics. Experiments on neutrinos from the sun have shown that not enough of them reach the earth to be consistent with the sun s energy output. Experiments on neutrinos produced in the atmosphere have shown that muon-type neutrinos seem to be depleted while passing through the earth compared to those incident on detectors directly from above. Neutrinos no longer seem to be the zero-mass passive participants in the sub-nuclear world that they were thought to be only a few years ago. Instead, the picture that is emerging shows a system of three, or possibly four, neutrino-mass eigenstates that mix in different amounts to form the familiar flavor eigenstates, the electron, muon and tau neutrino types. What is produced in a reaction is a flavor eigenstate which evolves as it travels because of a change in relative phases of the different mass eigenstates arising from their relative mass differences.

Each year brings new information that is beginning to fill in a still incomplete picture. The continuation of several current experiments and a host of soon to be begun experiments will further fill in the picture. However, enough is known now to indicate that certain of the interesting parameters, in particular the mixing angle known as $`\theta _{13}`$, the mass difference $`\mathrm{\Delta }m_{21}^2`$, the magnitude if any of CP violation in the neutrino sector, and the ordering of the neutrino masses cannot be well determined, if determined at all, from the current crop of running and planned experiments. Thus attention has begun to turn toward high-intensity accelerator sources to get further information.

Much of this attention has been focused on the potential of high-energy (20 GeV and above) muon storage rings as candidate sources.  These devices require considerable R&D to determine whether they are indeed feasible, and, if feasible, will be very expensive. In this paper I will discuss the potential of high-intensity, conventional neutrino beams to get at these small parameters. My conclusion is that low-energy, conventional muon-neutrino beams of attainable intensity are serious candidates for the next generation systems.

Before analyzing what can be done in the real world, it is useful to look at a primitive, two-neutrino, model to illustrate why low energy may be better than high energy in untangling the puzzle. Consider two species with a small mixing amplitude between them. The “signal-to-noise” ratio in an experiment looking for the appearance of neutrino type two in a beam of neutrino type one is given by

$$\frac{P(\nu _1\nu _2)}{P(\nu _1\nu _1)}=\frac{A^2\mathrm{sin}^2(\mathrm{\Delta }m^2L/4E)}{1A^2\mathrm{sin}^2(\mathrm{\Delta }m^2L/4E)}$$
(1)
where $`A`$ is the mixing amplitude, $`\mathrm{\Delta }m^2`$ is the difference of the squares of the masses, $`L`$ is the distance from the source to the detector, and $`E`$ is the beam energy. The optimum signal-to-noise ratio comes when the sine term is equal to one, i.e., $`\mathrm{\Delta }m^2L/4E`$ is equal to an odd integer multiple of $`\pi /2`$. However, all of the muon storage ring designs have high energy, making this factor small with the known mass difference. Hence, a very sophisticated background rejection is required in the detector.

It is also sometimes argued that, because neutrino cross sections increase with energy and the flux of neutrinos increases with the square of the energy of the parent particle, high-energy beams are better than low energy ones. However, in looking for the small terms that are unlikely to be determined by the present round of experiments, the appearance of a neutrino species different from the primary species gives the most sensitive tests. This probability is proportional to $`E^2`$, leaving only a single power of the energy as a potential advantage for the high-energy beams. Thus, this rationale for the choice of high energy is not as overwhelming as it might appear.

In addition, in the real three-neutrino world, there is a further complication from high-energy beams, tau-lepton production. In the scenarios using a storage-ring source, tau-lepton production is a potentially serious complication in determining the small parameters, since its leptonic and semi-leptonic decays can make a tau event look like an electron or a muon event. It is proposed to solve this problem by determining, for example, the sign of electrons produced by muon neutrinos  which would be effective but is very difficult in a large detector. It is a simplification if one can stay below the tau-lepton production threshold.

Now I turn to the standard three-neutrino formulation (recent results from Super Kamiokande appear to give 95% confidence in this formulation). Equation (3) includes the small terms that are often dropped in oscillation analyses. In Eq. (3), $`C_{ij}`$ means $`\mathrm{cos}\theta _{ij},S_{ij}`$ means $`\mathrm{sin}\theta _{ij},\mathrm{\Delta }m_{ij}^2`$ means $`m_i^2m_j^2`$, $`\delta `$ is the CP violating phase, $`L`$ is the distance from the source to the detector, and $`E`$ is the neutrino energy. The matter effect is given by

$$a=2\sqrt{2}G_Fn_eE=7.6\times 10^5\rho (gm/cm^3)E(GeV)(eV)^2$$
(2)
where $`G_F`$ is the Fermi weak-coupling constant, and $`n_e`$ is the electron density. Equation (3) is exact when the matter term of Eq. (2) is zero; otherwise $`a/\mathrm{\Delta }m_{31}^2`$ is assumed to be small.

$`P(\nu _\mu \nu _e)=4C_{13}^2S_{13}^2S_{23}^2\mathrm{sin}^2{\displaystyle \frac{\mathrm{\Delta }m_{31}^2L}{4E}}\times \left(1+{\displaystyle \frac{2a}{\mathrm{\Delta }m_{31}^2}}\left(12S_{13}^2\right)\right)`$
$`+8C_{13}^2S_{12}S_{13}S_{23}(C_{12}C_{23}\mathrm{cos}\delta S_{12}S_{13}S_{23})\mathrm{cos}{\displaystyle \frac{\mathrm{\Delta }m_{32}^2L}{4E}}\mathrm{sin}{\displaystyle \frac{\mathrm{\Delta }m_{31}^2L}{4E}}\mathrm{sin}{\displaystyle \frac{\mathrm{\Delta }m_{21}^2L}{4E}}`$
$`8C_{13}^2S_{13}^2S_{23}^2\mathrm{cos}{\displaystyle \frac{\mathrm{\Delta }m_{32}^2L}{4E}}\mathrm{sin}{\displaystyle \frac{\mathrm{\Delta }m_{31}^2L}{4E}}{\displaystyle \frac{aL}{4E}}\left(12S_{13}^2\right)`$ (3)
$`8C_{13}^2C_{12}C_{23}S_{12}S_{13}S_{23}\mathrm{sin}\delta \mathrm{sin}{\displaystyle \frac{\mathrm{\Delta }m_{32}^2L}{4E}}\mathrm{sin}{\displaystyle \frac{\mathrm{\Delta }m_{31}^2L}{4E}}\mathrm{sin}{\displaystyle \frac{\mathrm{\Delta }m_{21}^2L}{4E}}`$
$`+4S_{12}^2C_{13}^2\left\{C_{12}^2C_{23}^2+S_{12}^2S_{23}^2S_{13}^22C_{12}C_{23}S_{12}S_{23}S_{13}\mathrm{cos}\delta \right\}\mathrm{sin}^2{\displaystyle \frac{\mathrm{\Delta }m_{21}^2L}{4E}}.`$
It is not my purpose here to do an exhaustive analysis of all of the possibilities, but to show some examples of the potential, and in what follows I will assume that $`\theta _{23}=\theta _{12}=\pi /4`$. The Super Kamiokande atmospheric data favors this for $`\theta _{23}`$ and their solar data favor the large mixing-angle solution for the solar neutrino deficit. I will also take $`\mathrm{\Delta }m_{32}^2=3\times 10^3(eV)^2`$ which is the central value of the latest Super-K data,  $`\mathrm{\Delta }m_{31}^2\mathrm{\Delta }m_{32}^2`$, $`S_{13}^21`$, and the matter term $`(a)=0`$.

Now, look at the results of four experiments; $`\nu _\mu \nu _e`$ and $`\overline{\nu }_\mu \overline{\nu }_e`$ at two values of $`L/E`$ such that $`\mathrm{\Delta }m_{32}^2L/4E`$ equals $`\pi /2`$ and $`\pi `$. The result when the $`\mathrm{\Delta }m_{32}^2`$ term is $`\pi /2`$ is,

$$P_+(\pi /2)=P(\nu _\mu \nu _e)+P(\overline{\nu }_\mu \overline{\nu }_e)=4S_{13}^2+\mathrm{sin}^2\frac{\mathrm{\Delta }m_{21}^2L}{4E}$$
(4)
$$P_{}(\pi /2)=P(\nu _\mu \nu _e)P(\overline{\nu }_\mu \overline{\nu }_e)=4S_{13}\mathrm{sin}\delta \mathrm{sin}\frac{\mathrm{\Delta }m_{21}^2L}{4E}$$
(5)
and, when the mass term is $`\pi `$, is

$$P_+(\pi )=P(\nu _\mu \nu _e)+P(\overline{\nu }_\mu \overline{\nu }_e)=\mathrm{sin}^2\frac{\mathrm{\Delta }m_{21}^2L}{4E}$$
(6)
$$P_{}(\pi )=P(\nu _\mu \nu _e)P(\overline{\nu }_\mu \overline{\nu }_e)=0.$$
(7)
The next question is, are there enough neutrinos to do useful experiments. It will require a detailed design of a beam to get a precise answer, but we can get a rough idea by extrapolating from work already done at FNAL on the MINOS beam design and from the FNAL storage ring source study. The present FNAL proton beam power for the MINOS beam is about 400 kW. The beam intensity is limited by space-charge at injection into the 8-GeV Booster and at injection into the 100–GeV Main Injector. If the FNAL linac s energy was increased to about 1 GeV and the Booster s energy increased to about 16 GeV, the FNAL beam on target could be increased to about 4 MW.  The improvements to the injector chain are straightforward. The proton target is more difficult than the present system and will certainly require more shielding. Target R&D is being intensively pursued for the Spallation Neutron Source and should present no fundamental problems for neutrino production. The target for a conventional beam is somewhat easier per unit beam power than that for a muon storage ring source. The conventional beam is designed to produce pions of multi-GeV energy, while the muon source uses pions of low energy. For maximum production, the conventional source is thus thinner and absorbs less power than the storage ring source.

The comparative yields from a high-intensity, conventional, wide-band beam and a muon storage ring can be estimated by scaling the numbers given in the FNAL report. Table 1 shows the numbers. Using those numbers, a 20 kiloton detector at 732 km would see about $`9\times 10^4`$ events per year from the low-energy, wide-band, conventional beam.

At 732 km, the neutrino beam energies required to have $`(\mathrm{\Delta }m_{32}^2L/4E)`$ equal to $`\pi /2`$ and $`\pi `$ would be about 2 GeV and 1 GeV respectively ($`\pi `$ mesons of about 9 and 4.5 GeV). There is no design available for a narrow band beam, and so I use the spectrum of the MINOS low-energy beam to make a rough estimate. The spectrum has a low-energy peak and a high-energy tail, each generating about equal number of events  Much improved focusing can be achieved in a narrow band beam, which increases the flux. The peak in the MINOS spectrum has a width of about $`\pm `$ 30%, while it is desirable for background rejection to have a narrower width, thus reducing the flux. I will assume that the spectrum has a width of ($`\pm `$ 10–20% and that the focusing and width effects balance out in determining the beam intensity. Finally, I will give up another factor of two because this estimate may be too optimistic, or the FNAL Main Injector upgrade may not be capable of reaching the 4 MW power level. The result for neutrinos in the 20-kiloton detector at 732 km from the source is

$$Y_{cc}(2GeV)=2.5\times 10^4\mathrm{events}\mathrm{per}\mathrm{year}.$$
(8)
A 1–GeV beam will have a lower rate. There will be more pions in the decay channel, but the neutrino beam divergence will increase by a factor of two and the charged-current neutrino cross-section will decrease by a factor of two resulting in a reduction of about a factor of four in event rate

$$Y_{cc}(1GeV)=6\times 10^3\mathrm{events}\mathrm{per}\mathrm{year}.$$
(9)
Note that antineutrino cross-sections are about half of those for neutrinos.

Electron neutrinos in the beam come from the chain $`\pi \mu e`$ in the decay channel. They amount to about 0.1% total at the far detector, and amount to $`(12)\times 10^4`$ in the muon neutrino energy range. However, a more serious background will be neutral-current interactions initiated by muon and tau neutrinos in the detector. These can give $`\pi ^0`$ mesons, and highly asymmetric conversions of the $`\pi ^0`$ decay $`y`$’s can be confused with electron-neutrino-charged current events. These have to be controlled by detector design. Here I will assume that they give a background of 0.1% of the “unoscillated” muon-neutrino flux.

The optimum data-taking strategy depends on “a priori” knowledge. If, for example, KamLAND determines $`\mathrm{\Delta }m_{21}^2`$, there is no reason to take data at 1 GeV. Here, I will simply assume that data equivalent to 104 charged-current events (in the absence of oscillation) is collected for $`\nu _\mu `$ and $`\overline{\nu }_\mu `$ at both 2 GeV and 1 GeV. I also ignore the effect of the finite energy spread in the beam.

From Eq. (6), the minimum detectable value of $`\mathrm{\Delta }m_{21}^2`$ different from zero at the three standard-deviation level is

$$\left(\mathrm{\Delta }m_{21}^2\right)_{\mathrm{min}}=4.2\times 10^5(eV)^2.$$
(10)
This is independent of $`\theta _{13}`$ as long as $`C_{13}`$ is near one. It is near the lower end of the solar Large Mixing Angle solution allowed region and, if this large, might be measurable by the KamLAND experiment.

Equation (4) couples $`\mathrm{\Delta }m_{21}^2`$ and $`S_{13}^2`$. Table 2 gives the 3$`\sigma `$ lower bound for $`S_{13}^2`$ for various values of $`\mathrm{\Delta }m_{21}^2`$ ranging from the minimum detectable in this proposed experiment to the maximum in the solar LMA solution. What is happening, as the usually ignored $`\mathrm{\Delta }m_{21}^2`$ term increases, is that more events have to come from the $`S_{13}^2`$ term to meet the three standard deviation requirement.

The last item to look at is the sensitivity to the CP violating phase $`\delta `$. This is given in Eq. (5) and depends on both $`S_{13}`$ and $`\mathrm{\Delta }m_{21}^2`$ (note also that the sign of $`\mathrm{\Delta }m_{21}^2`$ comes in). If either of these is very small, CP violation becomes impossible to measure in this or any other experiment. To give an idea of the sensitivity, I will take $`\mathrm{\Delta }m_{21}^2=2\times 10^4`$, the center of its range, and take $`S_{13}^2`$ at the 3$`\sigma `$ limit for this $`\mathrm{\Delta }m_{21}^2`$, i.e. $`9\times 10^4`$. The minimum detectable value for the phase that differs from zero by three $`\sigma `$ is,

$$|\mathrm{sin}\delta |_{\mathrm{min}}=0.35$$
(11)
where Eq. (4) has been used to determine the error on the non-CP violating yield.

The determination of $`\delta `$ is most sensitive to the matter term. With Sato s  first order analysis, the matter term is 30% of the CP violating term for the parameters used above. I believe the matter effect needs to be taken to higher order. It is interesting to note that as the neutrino beam energy is increased, the matter term gets bigger while the CP violating term gets smaller, another argument for low-energy beams.

In summary, I have analyzed here the potential of high-intensity, low-energy, narrow-band conventional neutrino beams. The “gedanken” experiments outlined here give interesting limits on the measurement of the small terms among the neutrino-mixing parameters. These limits are better than those from the entry level storage ring neutrino factory.  It is well worth the time of the experts to see if my assumptions on potential beam intensity and purity, and background rejection are reasonable. If they are, these experiments can be carried out sooner, and at less cost than those with a muon storage ring source.