text:
# MPI–PhT 2003-27 Superheavy dark matter as UHECR source versus the SUGAR data

## 1 Introduction

Protons accelerated by uniformly distributed extragalactic astrophysical sources would be a perfect minimal explanation of the UHECR data above $`10^{19}`$ eV. However, protons with energy $`E>4\times 10^{19}`$ eV loose quickly energy due to pion production on cosmic microwave background photons. Thus the proton spectrum should show the so-called Greisen-Zatsepin-Kuzmin (GZK) cutoff , which is not observed by the AGASA experiment . Moreover, if the small-scale clusters in the arrival directions of UHECR measured by AGASA  are due to point-like sources, one can estimate their number . This number is so small that the nearest source should be located at the distance $`R_{\mathrm{min}}100`$ Mpc . This means that the GZK cutoff is exponentially sharp at $`E6\times 10^{19}`$ eV and even the data of the HiRes experiment  are inconsistent with the expected proton spectrum . This inconsistency becomes even stronger if BL Lacs, which show a statistically significant correlations with the arrival directions of UHECR with energy $`E46\times 10^{19}`$ eV , are sources of UHECR.

A possible solution to this problem would be the existence of superheavy dark matter (SHDM) . Superheavy particles with mass $`M_X10^{1314}`$ GeV can be naturally produced during inflation and would be today the dominant component of dark matter . Such particles will concentrate in galactic halos and the secondaries from their decay could be responsible for the highest energy cosmic rays. It has also been suggested that not decays but annihilations of SHDM particles produce the observed UHECRs .

By construction, this model has two clean signatures: the dominance of photons at the highest energies  and an anisotropy of the arrival directions with an increased flux from the Galactic center . Unfortunately, both signatures are not very clean for the present experiments. Indeed, at 95% C.L., $`30\%`$ of the UHECR above $`E>10^{19}`$ eV can be photons , which means that still most of UHECRs with $`E>4\times 10^{19}`$ eV can be photons without any contradiction to the experimental data. Since experiments in the northern hemisphere do not see the Galactic center, they are not very sensitive to a possible anisotropy of arrival directions of UHECR from SHDM. In contrast, the Galactic center was visible for the old Australian SUGAR experiment.

The anisotropy of the arrival directions using data from the full sky was discussed in Refs. . Reference  compared the flux from the Galactic center to the one from the anti-center and found them to be comparable. Similarly, the full-sky harmonic analysis including AGASA and SUGAR data from Ref.  found no significant anisotropy. In this work, we use a two-component energy spectrum of UHECRs consisting of protons from uniformly distributed, astrophysical sources and the fragmentation products of SHDM calculated in SUSY-QCD. We compare their expected arrival direction distribution to the data of the SUGAR experiment using a Kolmogorov-Smirnov test. Contrary to the harmonic analysis, this test allows to quantify directly the (dis-) agreement of the measured distribution of arrival direction with the expected one in the SHDM model. We consider both decays and annihilations of SHDM.

The paper is organized as follows: In Section II we discuss the status of the SUGAR data. The contribution of SHDM to the UHECR spectrum is discussed in Section III. In Section IV we perform an harmonic analysis of the arrival directions measured by the SUGAR experiment. Then we use a Kolmogorov-Smirnov test in one and two dimensions to check the consistency of the SUGAR data with the probability distribution of arrival directions expected for SHDM in Section V. Finally we conclude in Section VI.

## 2 Assessment of the SUGAR data

In order to use the SUGAR energy spectrum correctly we compare their data given in Ref.  to the energy spectrum measured by the AGASA  and HiRes experiments  in Fig. 1. Rescaling the SUGAR energies calculated with the Hillas prescription by 15% downwards, $`E_{\mathrm{Hillas}}E_{\mathrm{Hillas}}/1.15`$, makes their data consistent with the ones from AGASA. The same is true for the HiRes spectrum if the energy is rescaled up-wards by 25% 111The aperture of HiRes is energy dependent; the rescaling we perform should be seen therefore just as a crude approximation. We have chosen arbitrary the AGASA spectrum as reference, but changing the overall normalization of the spectra affects our results only weakly. As seen from Fig. 1, the SUGAR spectrum has the ankle at the correct place around $`E10^{19}`$ eV and is consistent with the AGASA spectrum in the whole energy range. In particular, the SUGAR spectrum also does not show the GZK cutoff at the highest energies.

The rescaling of the SUGAR data downwards by 15% should be compared to the recent reevaluation of the energy conversion formula used in the Haverah Park experiment . In this reference, the relation between $`\rho (600)`$ and the primary energy has been recalculated using QGSJET  and compared to the original relation suggested by Hillas. The new calibration results in $`30\%`$ lower primary energies.

The SUGAR experiment was a very sparse array of detectors and its energy determination of each single event was therefore rather unprecise. Thus, we shall use as a statistical test later on a method which relies only on the total flux measured by the SUGAR array, but uses not the energy of each single event. Since after the rescaling of the energies measured by SUGAR, $`E_{\mathrm{Hillas}}E_{\mathrm{Hillas}}/1.15`$, its measured flux is consistent with newer experiments like AGASA and HiRes, we conclude that on average the energy determination in the SUGAR experiment was reliable.

The energy conversion formula used in SUGAR to connect the measured muon number $`N_\mu `$ with the primary energy assumes that the primary is a hadron. For photon primaries, predicted to be dominant in the SHDM model, the muon content of the shower is smaller by a factor 5 – 10 . Thus the energies of photon events is expected to be underestimated by the SUGAR experiment. The SUGAR spectrum shown in Fig. 1 would be unchanged at energies $`E<5\times 10^{19}`$ eV, i.e at energies where all three experiments agree after rescaling.

The angular acceptance $`\eta (\delta )`$ as function of declination $`\delta `$ averaged over time of an experiment at geographical latitude $`b`$ ($`b=30.5^{}`$ for SUGAR) observing showers with maximal zenith angle $`\theta _{\mathrm{max}}`$ is

$$\eta (\delta )_0^{\alpha _{\mathrm{max}}}𝑑\alpha \mathrm{cos}(\theta )\left[\mathrm{cos}(b)\mathrm{cos}(\delta )\mathrm{sin}(\alpha _{\mathrm{max}})+\alpha _{\mathrm{max}}\mathrm{sin}(b)\mathrm{sin}(\delta )\right]$$
(1)
where

$$\xi =\frac{\mathrm{cos}(\theta _{\mathrm{max}})\mathrm{sin}(b)\mathrm{sin}(\delta )}{\mathrm{cos}(b)\mathrm{cos}(\delta )}$$
(2)
and

$$\alpha _{\mathrm{max}}=\{\begin{array}{c}\mathrm{arccos}(\xi )\mathrm{for}1\xi 1,\hfill \\ \pi \mathrm{for}\xi <1,\hfill \\ 0\mathrm{for}\xi >1.\hfill \end{array}$$
(3)
We have checked that the zenith angle distribution of the SUGAR events agrees with the theoretical predicted one, $`dN_{\mathrm{th}}d\theta \mathrm{sin}(\theta )\mathrm{cos}(\theta )`$, above $`E>4\times 10^{19}`$ eV. At lower energies, the acceptance of the experiment becomes energy dependent and deviations from $`dN_{\mathrm{th}}`$ start to grow.

## 3 Superheavy dark matter contribution to UHECR spectrum

We fix the contribution of SHDM to the total UHECR flux following the assumptions of Ref. : we assume that no galactic astrophysical sources contribute to the cosmic ray flux above $`10^{19}`$ eV and that the extragalactic cosmic ray flux can be characterized by an injection spectra of protons with a single power law, $`j_{\mathrm{ex}}(E)E^\alpha `$. For the choice of $`\alpha =2.7`$, this energy spectra modified by redshift, $`e^+e^{}`$ and pion production fits very well the measured spectra below $`E<68\times 10^{18}`$ eV . The only difference with  is that we take into account that total number of sources is small , if the small scale clusters measured by the AGASA experiment are due to point-like sources. The AGASA data favor as minimal distance to the nearest source $`D_{\mathrm{min}}100`$ Mpc . Therefore, the contribution of protons from extragalactic sources has a sharp cutoff. For the calculations of the proton spectrum we used the code . We use then the SUSY QCD fragmentation functions $`D(x,M_x)`$ of superheavy particles with mass $`M_X`$ calculated in Ref.  to model the flux $`j_{\mathrm{DM}}(E)D_\gamma (x,M_X)`$. The total UHECR flux is thus

$$j(E)=(1ϵ)j_{\mathrm{ex}}(E)+ϵj_{\mathrm{DM}}(E).$$
(4)
We fix the constant $`ϵ`$ determining the relative contribution of SHDM to the UHECR flux by a fit of $`j(E)`$ to the AGASA data . In Fig. 2, we show our fits for the case of a harder $`1/E^{2.3}`$ (left) and a softer $`1/E^{2.7}`$ (right) injection spectrum of extragalactic protons. In the first case, the contribution of SHDM to the UHECR spectrum below the GZK cutoff is minimal and starts to dominate only at highest energies $`E>6\times 10^{19}`$ eV. For this choice of injection spectrum, the contribution of galactic sources dominate for $`E<10^{19}`$ eV. In the second case, SHDM gives a larger contribution at lower energies $`E<6\times 10^{19}`$ than before, and again starts to dominate at $`E>6\times 10^{19}`$ eV. Note, that most UHECRs above $`E>6\times 10^{19}`$ eV should be photons in this model, but this does not contradict the rather weak existing bound of 30% of photons at $`E>10^{19}`$ eV .

In the following, we shall use conservatively the case of the harder $`1/E^{2.3}`$ injection spectrum if not otherwise stated. Then the contribution of SHDM to the UHECR spectrum is fixed by Eq. (4).

## 4 Harmonic analysis

In order to compare SUGAR data with an uniform distribution typical for extragalactic sources we have performed an one-dimensional harmonic analysis. As usual we sum

$$a_k=\frac{2}{n}\underset{a=1}{\overset{n}{}}\mathrm{cos}(k\varphi )\text{and}b_k=\frac{2}{n}\underset{a=1}{\overset{n}{}}\mathrm{sin}(k\varphi )$$
(5)
over the $`n`$ data points.

The amplitude $`r_k`$ and and phase $`\varphi _k`$ of the $`k`$.th harmonic are given by

$$r_k=\sqrt{a_k^2+b_k^2}\text{and}\varphi _k=\mathrm{arctan}(b_k/a_k)$$
(6)
with chance probability

$$p_{\mathrm{ch}}=\mathrm{exp}\left(nr_k^2/4\right).$$
(7)
The direction to the signal is $`\varphi =k\varphi _k`$.

Results of a harmonic analysis in right ascension $`\alpha `$ depending on $`E_{\mathrm{min}}`$ and $`\theta _{\mathrm{max}}`$ are given in Table 1. The results for all harmonics show generally good agreement with an isotropic distribution for any cutoff energy we have used. Only the third harmonics shows some anisotropy, in particular at the highest energies, $`E>8\times 10^{19}`$ eV; however its phase does not points towards the galactic center (lying at $`\alpha =266^{}`$). Generally, all harmonics point instead towards $`\alpha 130^{}`$. Reference  derived the probability distribution (pdf) that a data set with phase $`\varphi _1`$ and amplitude $`r_1`$ is drawn from an arbitrary pdf. However, we are not aware of a generalization to higher harmonics and, in particular, of a method to combine the information content of several harmonics. In the next section, we use therefore a Kolmogorov-Smirnov test to quantify the (dis-) agreement between the expected distribution of arrival direction and the SUGAR data.

## 5 Kolmogorov-Smirnov tests

The pdf to detect an event with energy $`E`$ and arrival direction $`\alpha ,\delta `$ is a combination of the isotropic extragalactic and the SHDM flux,

$$p(E,\alpha ,\delta )\eta (\delta )\left[j_{\mathrm{ex}}(E)+j_{DM}(E)_0^{s_{\mathrm{max}}}𝑑sn_{\mathrm{DM}}(r(\alpha ,\delta ))\right]$$
(8)
where $`s_{\mathrm{max}}=R_E\mathrm{cos}\theta +\sqrt{R_h^2R_E^2\mathrm{sin}^2\theta }`$ is given by the extension $`R_h=100`$ kpc of the DM halo and $`\theta `$ is the angle relative to the direction to the GC. As explained, the relative size of the two contributions is fixed by the fit to Eq. (4).

For the two-dimensional test, we have integrated Eq. (8) over energy,

$$P_{2d}(\alpha ,\delta )=_{E_{\mathrm{min}}}^{E_{\mathrm{max}}}𝑑Ep(E,\alpha ,\delta ),$$
(9)
where $`E_{\mathrm{min}}`$ and $`E_{\mathrm{max}}`$ are the minimal and maximal energy considered in the UHECR spectrum. We have used $`E_{\mathrm{max}}=10^{2122}`$ eV, but the results do depend on very weakly on the exact value. By contrast, the value of $`E_{\mathrm{min}}`$ has a strong influence on the results obtained.

For the one-dimensional test, we have integrated Eq. (9) over the declination,

$$P_{\mathrm{CR}}(\alpha )=_{\pi /2}^{\pi /2}𝑑\delta \mathrm{cos}\delta P_{2d}(\alpha ,\delta ).$$
(10)
In the standard one-dimensional Kolmogorov-Smirnov (KS) test, the maximal difference $`D`$ between the cumulative probability distribution function $`P(x)=𝑑x^{}p(x^{})`$ and the cumulative distribution function of the data,

$$S(x)=\frac{1}{n}\underset{i}{}\theta (x_ix),$$
(11)
is used as estimator for the belief that the data are drawn from the distribution $`p`$. A variant of this test which is equally sensitive on differences for all $`x`$ and is especially well suited for data on $`S^1`$ uses instead of $`D`$ the symmetric estimator

$$V=D_++D_{}=\mathrm{max}[S(x)P(x)]+\mathrm{max}[P(x)S(x)].$$
(12)
The significance of a certain value of $`V`$ is calculated with the formula given in . Since the exposure of a ground-array experiment is uniform in right ascension $`\alpha `$, we use $`\alpha `$ as variable in the one-dimensional KS test. More exactly, we use as pdf Eq. (8) integrated over $`dE`$ and $`d\delta \mathrm{cos}\delta `$.

As simplest test, we assume that all SUGAR events above $`E_{\mathrm{min}}`$ are produced by SHDM. Thus we compare $`S(\alpha )`$ with $`P(\alpha )=_0^\alpha 𝑑\alpha ^{}P_{\mathrm{CR}}(\alpha ^{})`$. The result is shown for $`E>8\times 10^{19}`$ eV in Fig. 3a for two different values of the maximal zenith angle, $`\theta _{\mathrm{max}}=45^{}`$ and $`\theta _{\mathrm{max}}=55^{}`$. While for $`\theta _{\mathrm{max}}=45^{}`$ SHDM is disfavoured at the two sigma level for realistic values of the core radius, $`R_c20`$ kpc, the SUGAR data have for the choice of $`\theta _{\mathrm{max}}=55^{}`$ a rather large probability $`p`$ to be consistent with the SHDM hypothesis, $`p20\%`$. In Fig. 3b, we compare the two commulative distributions $`S(\alpha )`$ and $`P(\alpha )`$ for $`R_s=15`$ kpc, $`E_{\mathrm{min}}=8\times 10^{19}`$ eV and $`\theta _{\mathrm{max}}=45^{}`$. Inspecting $`S(\alpha )`$ makes it clear that the data in this case are not uniformly distributed but clustered around $`\alpha 130^{}`$ and $`\alpha 350^{}`$. Since none of these two directions coincide with the position of the GC, this data set disfavours the SHDM hypothesis more strongly than one would expect for uniformly distributed events from extragalactic sources. However, one should use a rather low value of $`E_{\mathrm{min}}`$ to minimize the uncertainties in the SUGAR energy determination and we will therefore not rely on these results.

We consider therefore next as more realistic test the case that both SHDM and extragalactic sources contribute to the UHECR spectrum. Then the dependence of $`p`$ on $`E_{\mathrm{min}}`$ should be diminished. More exactly, one would expect in the case that the SHDM hypothesis is disfavoured by the data that decreasing $`E_{\mathrm{min}}`$ first decreases $`p`$. This decrease of $`p`$ should continue down until $`E(34)\times 10^{19}`$ eV, i.e. until a point where the signal-to-background ratio becomes considerably smaller than one. Decreasing $`E_{\mathrm{min}}`$ even further should result in a increase of $`p`$ because now practically all new events are from extragalactic sources.

In Fig. 4a, we show the dependence of the probability on the energy cutoff for $`R_c=15`$ kpc and $`\theta _{\mathrm{max}}=55^{}`$. The two thick solid lines show $`p`$ for a combination of SHDM and uniform sources according Eq. (4); the upper one corresponds to an injection spectrum $`1/E^{2.3}`$, the lower one to an injection spectrum $`1/E^{2.7}`$. The behaviour of $`p`$ suggest that the minimum for $`E_{\mathrm{min}}7\times 10^{19}`$ eV is a fluctuation similar to the maximum around $`E_{\mathrm{min}}5\times 10^{19}`$ eV. In the range $`E_{\mathrm{min}}(34)\times 10^{19}`$ eV, the fluctuations adding an additional event are relatively small. Therefore, we consider the probability in this range as more reliable indicator for the consistency of SHDM with the SUGAR arrival directions; we conclude that the SUGAR data have the probability $`p=520\%`$ to be consistent with the SHDM depending on the injection spectrum of the extragalactic protons.

The thin solid line shows how consistent the SUGAR data are with an isotropic distribution. This distribution has also a minimum around $`E_{\mathrm{min}}7\times 10^{19}`$ eV where the events cluster around two arrival directions. After including more low-energy data, the SUGAR arrival directions are consistent with an isotropic distribution. Finally, the dashed line shows the consistency of the SUGAR data with the assumption that all UHECR events above $`E_{\mathrm{min}}`$ are from SHDM. It is clear that only values of $`E_{\mathrm{min}}`$ above $`E_{\mathrm{min}}4\times 10^{19}`$ eV are compatible with the SUGAR data. Similar, the spectral shape of the flux in the SHDM model allows a dominance of SHDM in the UHECR spectrum only above $`E>6\times 10^{19}`$ eV . In Fig. 4b we compare results from one- and two-dimensional KS tests (of $`\alpha `$ and $`\delta `$) as function of the energy cutoff and find that they give rather similar results.

Finally, we consider the model where not decays but annihilations of SHDM particles produce the observed UHECRs . In the original version of this model it was suggested that the flux of the clumpy component dominates over the one from the smooth SHDM profile by 3 orders of magnitude. On the other side, it was shown in a recent paper  that the contribution of the clumpy component can be just a factor few larger than the one of the smooth component. Moreover, the newest numerical calculations show that the contribution of clumps is even subdominant and that it is very unlikely that a nearby clump will outshine the Galactic center . Because of the arguments above, we assume that the clumpy part of the SHDM gives a subdominant contribution to the UHECR flux. In the opposite case our results for the SHDM model with annihilations will be less significant, depending on the relative contribution of the two components.

Since the flux is now $`n_{\mathrm{DM}}^2`$, the anisotropy in this model is much stronger than for decaying SHDM. This can be clearly seen in Fig. 5b where we show the dependence of the probability of annihilating SHDM on the core radius $`R_c`$ assuming that all events above $`E`$ are from SHDM. Even for core radii as large as 30 kpc, annihilating SHDM is disfavoured by two sigma. Figure 5a shows similar to Fig. 4a the dependence of the probability on the energy cutoff for $`R_c=15`$ kpc and $`\theta _{\mathrm{max}}=55^{}`$. The two thick solid lines show $`p`$ for a combination of SHDM and uniform sources according Eq. (4); the upper one corresponds to an injection spectrum $`1/E^{2.3}`$, the lower one to an injection spectrum $`1/E^{2.7}`$. Depending on the injection spectrum of extragalactic protons, annihilations of SHDM are disfavoured by the SUGAR data between 3 and 4$`\sigma `$.

## 6 Conclusions

In this paper we have tested the consistency of the SHDM model with the SUGAR data. In order to use the SUGAR data, we have compared its energy spectrum to the one of AGASA and found that they are compatible after rescaling down the SUGAR energies by 15%. We have assumed that the energy spectrum in the region $`10^{19}\mathrm{eV}<E<6\times 10^{19}`$, i.e. between ankle and GZK cutoff, is dominated by protons coming from uniformly distributed extragalactic sources. After fitting the relative contributions of SHDM decay products and extragalactic protons to the AGASA data, we have performed Kolmogorov-Smirnov tests of the SUGAR data. As result we have found that SUGAR data are able to disfavour strongly extreme case like annihilations of SHDM without clumps ($`5\sigma `$) or decaying SHDM (99% CL) assuming their contribution to the UHECR flux dominates down to $`E=4\times 10^{19}`$ eV. The phenomenologically most interesting case, decaying SHDM dominating the UHECR spectrum only above $`E>6\times 10^{19}`$ eV, is consistent with the SUGAR date with 5–20% probability. Thus the SUGAR data do not disfavour strongly this model but they neither support it. A statistically significant test of this model can be done by the Pierre Auger Observatory.

## Acknowledgements

We would like to thank L. Anchordoqui, R. Dick, J. Primack, K. Shinozaki, M. Teshima, S. Troitsky and S. White for discussions and comments. This work was supported by the Deutsche Forschungsgemeinschaft (DFG) within the Emmy Noether program.

Note added: After completing this work we have learnt about a preprint of H. B. Kim and P. Tinyakov  discussing also the consequences of the SUGAR data for the SHDM hypothesis. These authors come to similar conclusion as ours. We thank H. B. Kim and P. Tinyakov for sending us the preprint before publication.


format: text

text:
# How to derive and compute 1,648 diagrams

## Abstract

We present the first calculation for many-electron atoms complete through fourth order of many-body perturbation theory. Owing to an overwhelmingly large number of underlying diagrams, we developed a suite of symbolic algebra tools to automate derivation and coding. We augment all-order single-double excitation method with 1,648 omitted fourth-order diagrams and compute amplitudes of principal transitions in Na. The resulting ab initio relativistic electric-dipole amplitudes are in an excellent agreement with 0.05%-accurate experimental values. Analysis of previously unmanageable classes of diagrams provides a useful guide to a design of even more accurate, yet practical many-body methods.

preprint: UNR Apr 2003-

Many-body perturbation theory (MBPT) has proven to be a powerful tool in physics Fetter and Walecka (1971) and quantum chemistry Szabo and Ostlund (1982). Although MBPT provides a systematic approach to solving many-body quantum-mechanical problem, the number and complexity of analytical expressions and thus challenges of implementation grow rapidly with increasing order of MBPT (see Fig.1) Indeed, because of this complexity it has proven to be difficult to go beyond the complete third order in calculations for many-electron atoms (see, e.g., review Sapirstein (1998)). At the same time, studies of higher orders are desirable for improving accuracy of ab initio atomic-structure methods. Such an improved accuracy is required, for example, in interpretation of atomic parity violation PNC  and unfolding cosmological evolution of the fine-structure constant $`\alpha `$ Webb et al. (2001).

Here we report the first calculation of transition amplitudes for alkali-metal atoms complete through the fourth order of MBPT. We explicitly computed 1,648 topologically distinct Brueckner-Goldstone diagrams. To overcome such an overwhelming complexity we developed a symbolic problem-solving environment that automates highly repetitive but error-prone derivation and coding of many-body diagrams. Our work illuminates a crucial role of symbolic tools in problems unsurmountable in the traditional “pencil-and-paper” approach. Indeed, third-order calculation Blundell et al. (1987); Johnson et al. (1996) was a major research project most likely to have required a year to accomplish. As one progresses from the third to the present fourth order (see Fig. 1) there is a 50-fold increase in the number of diagrams. Simple scaling shows that present calculations require half a century to complete. With our tools derivation and coding take just a few minutes. Similar symbolic tools were developed by Perger et al. (2001), however their package is so far limited to well-studiedBlundell et al. (1990a); Johnson et al. (1996) third order of MBPT. In contrast, we explore a wide range of new, previously unmanageable, classes of diagrams.

As an example of application of our symbolic technology, we compute electric-dipole amplitudes of the principal $`3p_{3/2}3s_{1/2}`$ and $`3p_{1/2}3s_{1/2}`$ transitions in Na. We augment all-order single-double excitations method Safronova et al. (1998) with 1,648 diagrams so that the formalism is complete through the fourth order ( see Fig.1 ). The results are in excellent agreement with 0.05%-accurate experimental values Jones et al. (1996). Thus our computational method not only enables exploration of a wide range of previously unmanageable classes of diagrams but also defines new level of accuracy in ab initio relativistic atomic many-body calculations. Atomic units $`|e|=\mathrm{}=m_e=4\pi \epsilon _01`$ are used throughout this paper.

Method— A practitioner of atomic MBPT typically follows these steps: (i) derivation of diagrams, (ii) angular reduction, and (iii) coding and numerical evaluation. Below we highlight components of our problem-solving environment designed to assist a theorist in these technically-involved tasks. First we briefly reiterate MBPT formalism Derevianko and Emmons (2002) for atoms with a single valence electron outside a closed-shell core. For these systems a convenient point of departure is a single-particle basis generated in frozen-core ($`V^{N1}`$) Dirac-Hartree-Fock (DHF) approximation Kelly (1969). In this approximation the number of MBPT diagrams is substantially reduced Lindgren and Morrison (1986); Blundell et al. (1987). The lowest-order valence wavefunction $`|\mathrm{\Psi }_v^{(0)}`$ is simply a Slater determinant constructed from core orbitals and proper valence state $`v`$. The perturbation expansion is built in powers of residual interaction $`V_I`$ defined as a difference between the full Coulomb interaction between the electrons and the DHF potential. The $`n^{\mathrm{th}}`$-order correction to the valence wavefunction may be expressed as

$$|\mathrm{\Psi }_v^{(n)}=R_v\left\{QV_I|\mathrm{\Psi }_v^{(n1)}\right\}_{\mathrm{linked}},$$
(1)
where $`R_v`$ is a resolvent operator modified to include so-called “folded” diagrams Derevianko and Emmons (2002), projection operator $`Q=1|\mathrm{\Psi }_v^{(0)}\mathrm{\Psi }_v^{(0)}|`$, and only linked diagrams Lindgren and Morrison (1986) are to be kept. From this recursion relation we may generate corrections to wave functions at any given order of perturbation theory. With such calculated corrections to wavefunctions of two valence states $`w`$ and $`v`$, $`n^{\mathrm{th}}`$-order contributions to matrix elements of an operator $`\widehat{Z}`$ are

$$Z_{wv}^{(n)}=\underset{k=0}{\overset{n1}{}}\mathrm{\Psi }_w^{(nk1)}|Z|\mathrm{\Psi }_v^{(k)}_{\mathrm{val},\mathrm{conn}}+Z_{wv,\mathrm{norm}}^{(n)}.$$
(2)
Here $`Z_{wv,\mathrm{norm}}^{(n)}`$ is a normalization correction arising due to an intermediate normalization scheme employed in derivation of Eq. (1). Subscript “$`\mathrm{val},\mathrm{conn}`$” indicates that only connected diagrams involving excitations from valence orbitals are included in the expansion.

Equations (1) and (2) completely define a set of many-body diagrams at any given order of MBPT. In practice the derivations are carried out in the second quantization and the Wick’s theorem Lindgren and Morrison (1986) is used to simplify products of creation and annihilation operators. Although the application of the Wick’s theorem is straightforward, as order of MBPT increases, the sheer length of expressions and number of operations becomes quickly unmanageable. We developed a symbolic-algebra package written in Mathematica Wolfram (1999) to carry out this task. The employed algorithm relies on decision trees and pattern matching, i.e., programming elements typical to artificial intelligence applications. With the developed package we fully derived fourth-order corrections to matrix elements of univalent systems Derevianko and Emmons (2002).

This is one of the fourth-order terms from Ref. Derevianko and Emmons (2002)

$$\underset{abc}{}\underset{mnr}{}\frac{z_{bv}\stackrel{~}{g}_{canr}g_{nrcm}\stackrel{~}{g}_{mwab}}{(\epsilon _w\epsilon _b)(\epsilon _{mw}\epsilon _{ab})(\epsilon _{nrw}\epsilon _{abc})}.$$
(3)
There are 524 such contributions in the fourth order nmb . Here abbreviation $`\epsilon _{xy\mathrm{}z}`$ stands for $`\epsilon _x+\epsilon _y+\mathrm{}\epsilon _z`$, with $`\epsilon _x`$ being single-particle DHF energies. Further, $`g_{ijlk}`$ are matrix elements of electron-electron interaction in the basis of DHF orbitals. The quantities $`\stackrel{~}{g}_{ijlk}`$ are antisymmetric combinations $`\stackrel{~}{g}_{ijlk}=g_{ijlk}g_{ijkl}`$. The summation is over single-particle DHF states. Core orbitals are enumerated by letters $`a,b,c`$ and complementary excited states are labelled by $`m,n,r`$. Finally matrix elements of the operator $`\widehat{Z}`$ in the DHF basis are denoted as $`z_{ij}`$.

The summations over magnetic quantum numbers are usually carried out analytically. This “angular reduction” is the next major technically-involved step. We also automate this task. The details are provided in Ref. Derevianko . Briefly, the angular reduction is based on application of the Wigner-Eckart (WE) theorem Edmonds (1985) to matrix elements $`z_{ij}`$ and $`g_{ijkl}`$. The WE theorem allows one to “peel off” dependence of the matrix elements on magnetic quantum numbers in the form of 3j-symbols and reduced matrix elements. In the particular case of fourth-order terms, such as Eq. (3), application of the WE theorem results in a product of seven 3j-symbols. To automate simplification of the products of 3j-symbols we employed a symbolic program Kentaro developed by Takada (1992).

The result of angular reduction of our sample term (3) is

$`{\displaystyle \underset{abcmnr}{}}{\displaystyle \underset{L}{}}{\displaystyle \frac{\delta _{j_aj_m}\delta _{j_bj_w}\left(1\right)^{j_a+j_c+j_n+j_r}}{(2L+1)\sqrt{(2j_a+1)(2j_w+1)}}}`$
$`{\displaystyle \frac{bzvZ_L(canr)X_L(nrcm)Z_0(mwab)}{(\epsilon _w\epsilon _b)(\epsilon _{mw}\epsilon _{ab})(\epsilon _{nrw}\epsilon _{abc})}}.`$
Here the reduced quantities $`izj`$, $`X_L(ijkl)`$, and $`Z_L(ijkl)`$ depend only on total angular momenta and principal quantum numbers of single-particle orbitals.

As a result of angular reduction we generate analytical expressions suitable for coding. We also automated the tedious coding process by developing custom parsers based on Perl and Mathematica. These parsers translate analytical expressions into Fortran90 code. The resulting code is very large - it is about 20,000 lines long and were it be programmed manually, it would have required several years to develop. For numerical evaluation we employed a B-spline libraryJohnson et al. (1988). All the fourth-order results were computed with a sufficiently large basis of 25 out of 30 lowest-energy ($`E>mc^2`$) spline functions for each partial wave through $`h_{11/2}`$.

At this point we have demonstrated feasibility of working with thousands of diagrams in atomic MBPT. Now we apply our computational technique to high-accuracy calculation of transition amplitudes in Na.

Fourth-order diagrams complementary to single-double excitations method. One of the mainstays of practical applications of MBPT is an assumption of convergence of series in powers of the perturbing interaction. Sometimes the convergence is poor and then one sums certain classes of diagrams to “all orders” using iterative techniques. In fact, the most accurate many-body calculations of parity violation in Cs by Dzuba et al. (1989) and Blundell et al. (1990b) are of this kind. These techniques, although summing certain classes of MBPT diagrams to all orders, still do not account for an infinite number of residual diagrams (see Fig. 1). In Ref. Derevianko and Emmons (2002) we proposed to augment a given all-order technique with some of the omitted diagrams so that the formalism is complete through a certain order of MBPT. As in that work, here we consider an improvement of all-order single-double (SD) excitation method employed in Ref. Blundell et al. (1990b). Here a certain level $`n`$ of excitations from lowest-order wavefunction refers to an all-order grouping of contributions in which $`n`$ core and valence electrons are promoted to excited single-particle orbitals. The SD method is a simplified version of the coupled-cluster expansion truncated at single and double excitations.

The next step in improving the SD method would be an inclusion of triple excitations. However, considering present state of available computational power, the complete incorporation of triples seems as yet impractical for heavy atoms. Here we investigate an alternative illustrated in Fig. 1 : we compute the dominant contribution of triples in a direct fourth-order MBPT for transition amplitudes. We also account for contribution of disconnected quadruple excitations in the fourth order. In Ref. Derevianko and Emmons (2002), we separated these complementary diagrams into three major categories by noting that triples and disconnected quadruples enter the fourth order matrix element $`Z_{wv}^{(4)}`$ via (i) an indirect effect of triples and disconnected quadruples on single and double excitations in the third-order wavefunction — we denote this class as $`Z_{0\times 3}`$; (ii) direct contribution to matrix elements labelled as $`Z_{1\times 2}`$; (iii) correction to normalization denoted as $`Z_{\mathrm{norm}}`$. Further these categories were broken into subclasses based on the nature of triples, so that

$`\left(Z_{wv}^{(4)}\right)_{\mathrm{non}\mathrm{SD}}=Z_{1\times 2}(T_v)+Z_{1\times 2}(T_c)+`$
$`Z_{0\times 3}(S_v[T_v])+Z_{0\times 3}(D_v[T_v])+`$
$`Z_{0\times 3}(S_c[T_c])+Z_{0\times 3}(D_v[T_c])+`$
$`Z_{1\times 2}\left(D_{\mathrm{n}l}\right)+Z_{0\times 3}\left(D_{\mathrm{n}l}\right)+Z_{\mathrm{norm}}(T_v).`$
Here we distinguished between valence ($`T_v`$) and core ($`T_c`$) triples and introduced a similar notation for singles ($`S`$) and doubles ($`D`$). Notation like $`S_v[T_c]`$ stands for effect of second-order core triples ($`T_c`$) on third-order valence singles $`S_v`$. Diagrams $`D_{\mathrm{nl}}`$ are contributions of disconnected quadruples (non-linear contributions from double excitations). The reader is referred to Ref. Derevianko and Emmons (2002) for further details and discussion.

Transition amplitudes in Na. Using our problem-solving environment we derived the 1,648 complementary diagrams Derevianko and Emmons (2002), carried out angular reduction Derevianko , and generated Fortran 90 code suitable for any univalent system. As an example we evaluate reduced electric-dipole matrix elements of $`3s_{1/2}3p_{1/2,3/2}`$ transitions in Na (eleven electrons)com (a). Our numerical results are presented in Table 1. Analyzing this table we see that leading contributions come from valence triples $`T_v`$. Similar conclusion can be drawn from our preliminary calculations for heavier Cs atom. Dominance of valence triples ($`T_v`$) over core triples ($`T_c`$) may be explained by smaller energy denominators for $`T_v`$ terms. Representative diagrams for these relatively large contributions are shown in Fig. 2. Based on this observation we propose to fully incorporate valence triples into a hierarchy of coupled-cluster equations and add a perturbative contributions of core triples. Such an all-order scheme would be a more accurate and yet practical extension of the present calculations.

Another point we would like to discuss is a sensitivity of our results to higher-order corrections. In Table 1, all large contributions add up coherently, possibly indicating a good convergence pattern of MBPT. However, we found large, factor of 100, cancellations of terms inside the $`Z_{0\times 3}(S_v[T_v])`$ class. In principle higher-order MBPT corrections may offset a balance between cancelling terms and an all-order treatment is desired. Fortunately, the $`Z_{0\times 3}(S_v[T_v])`$ class of diagrams (combined with parts of $`Z_{1\times 2}(T_v)`$) have been taken into account in all-order SDpT (SD + partial triples) method Blundell et al. (1990b); Safronova et al. (1999). We correct our results for the difference between all-order com (b) and our fourth-order values for these diagrams (last row of Table 1). These all-order corrections modify our final values of complementary diagrams by 15%.

In Table 2 we add our complementary diagrams to SD matrix elements Safronova et al. (1998) and compare with experimental values. Several high-accuracy experiments have been carried out for Na, resolving an apparent disagreement between an earlier measurement and calculated lifetimes (see review Volz and Schmoranzer, 1996, and references therein). In Table 2 we compare with the results of the two most accurate experimentsJones et al. (1996); Volz et al. (1996). The SD method Safronova et al. (1998) overestimates these experimental values by 2.5 $`\sigma `$ and 2.8 $`\sigma `$ respectively ($`\sigma `$ is experimental uncertainty). With our fourth-order corrections taken into consideration the comparison significantly improves. The resulting ab initio matrix elements for both $`3p_{1/2}3s_{1/2}`$ and $`3p_{3/2}3s_{1/2}`$ transitions are in an excellent agreement with 0.05%-accurate values from Ref. Jones et al. (1996) and differ by 1.2$`\sigma `$ from less-accurate results of Ref. Volz et al. (1996). Considering this agreement it would be desirable to have experimental data accurate to 0.01%.

We demonstrated that symbolic tools can replace multi-year detailed development efforts in atomic MBPT with an interactive work at a conceptual level, thus enabling an exploration of ever more sophisticated techniques. As an example, we presented the first calculations for many-electron atoms complete through fourth order, a task otherwise requiring half a century to complete. Although even at this level the computed transition amplitudes for Na indicate a record-setting ab initio accuracy of a few 0.01%, the calculations allowed us to gain insights into relative importance of various contributions and to propose even more accurate yet practical many-body method. With an all-order generalization Derevianko and Emmons (2002) of the derived diagrams we plan to address a long-standing problem Dzuba et al. (1989); Blundell et al. (1990b) of improving theoretical accuracy of interpretation of parity violation in Cs atom Wood et al. (1997).

We would like to thank W.R. Johnson, V.A. Dzuba, W.F. Perger, and K. Takada for useful discussions and M.S. Safronova for providing detailed breakdown of SDpT and SD results. This work was supported in part by the National Science Foundation.


text:
# References

Schrödinger Equation for Nanoscience

Janina Marciak-Kozlowska
Institute of Electron Technology, Al. Lotników 32/46 02–668 Warsaw Poland

Miroslaw Kozlowski
Physics Department, Science Teachers College and Institute of Experimental Physics,
Warsaw University, Hoża 69, 00-681 Warsaw Poland
email: mirkoz@fuw.edu.pl

## Abstract

The second order (in time) Schrödinger equation is proposed. The additional term (in comparison to Schrödinger equation) describes the interaction of particles with vacuum filled with virtual particle – antiparticle pairs (zitterbewegung).

Key words: Schrödinger equation; Nanoscience; Zitterbewegung.

Quantum mechanics has been remarkably successful in all realms of atoms molecular and solids. But even more remarkable is the fact that quantum theory still continues to fascinate researches. Interest in quantum mechanics both theoretical and experimental is probably greater now that it ever has been.

In this article we develop the modified Schrödinger equation which describes the structure of matter on the subatomic level i.e. for characteristic dimension $`r_n<d<r_a`$ where $`r_n(\text{nucleus radius})\text{fm}`$, $`r_a(\text{atom radius})\text{nm}`$. To that aim we use the analogy between the Schrödinger equation and diffusion equation (Fourier equation). The quantum Fourier equation which describes the heat (mass) diffusion on the atomic level has the form :

$$\frac{T}{t}=\frac{\mathrm{}}{m}^2T.$$
(1)
When the real time $`tit/2`$ and $`T\mathrm{\Psi }`$, Eq. (1) has the form of the free Schrödinger equation:

$$i\mathrm{}\frac{\mathrm{\Psi }}{t}=\frac{\mathrm{}^2}{2m}^2\mathrm{\Psi }.$$
(2)
The complete Schrödinger equation has the form:

$$i\mathrm{}\frac{\mathrm{\Psi }}{t}=\frac{\mathrm{}^2}{2m}^2\mathrm{\Psi }+V\mathrm{\Psi }$$
(3)
where $`V`$ denotes the potential energy. When we go back to real time $`t2it`$, $`\mathrm{\Psi }T`$ the new parabolic quantum heat transport equation (quantum Fokker-Planck equation) is obtained:

$$\frac{T}{t}=\frac{\mathrm{}}{m}^2T\frac{2V}{\mathrm{}}T.$$
(4)
Equation (4) describes the quantum heat transport for $`\mathrm{}t>\tau `$. For ultrashort time processes, $`\mathrm{}t<\tau `$ one obtains the generalized quantum hyperbolic heat transport equation:

$$\tau \frac{^2T}{t^2}+\frac{T}{t}=\frac{\mathrm{}}{m}^2T\frac{2V}{\mathrm{}}T.$$
(5)
The structure and the solutions of Eq. (5) for ultrashort thermal processes was investigated in the monograph: M. Kozlowski, J. Marciak-Kozlowska: From quarks to bulk matter, Hadronic Press, USA, 2001. The generalized heat transport equation (5) leads to modified Schrödinger equation (MSE). After substitution $`tit/2,T\mathrm{\Psi }`$ in Eq. (5) one obtains:

$$i\mathrm{}\frac{\mathrm{\Psi }}{t}=\frac{\mathrm{}^2}{2m}^2\mathrm{\Psi }+V\mathrm{\Psi }2\tau \mathrm{}\frac{^2\mathrm{\Psi }}{t^2}.$$
(6)
The additional term (in comparison to Schrödinger equation) describes the interaction of electrons with surrounding space-time filled with virtual positron-electron pairs, i.e. zitterbewegung.

One can conclude that for time period $`\mathrm{}t<\tau `$ and distance $`\mathrm{}r<c\tau `$ the description of quantum phenomena needs some revision. The numerical values for $`\mathrm{}t`$ and $`\mathrm{}r`$ can be calculated as follows. Considering that $`\tau `$

$`\tau `$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{m\alpha ^2c^2}}10^{17}\mathrm{s}`$ (7)
$`c\tau `$ $`=`$ $`{\displaystyle \frac{\mathrm{}c}{m\alpha ^2c^2}}1\mathrm{nm},`$
we conclude that with the help of MSE we can visit the inner atomic environement. On the other hand for $`\mathrm{}t>10^{17}`$ s, $`\mathrm{}r>`$ nm, in MSE the second derrivative term can be omitted and as result the SE is obtained, i.e.

$$i\mathrm{}\frac{\mathrm{\Psi }}{t}=\frac{\mathrm{}^2}{2m}^2\mathrm{\Psi }+V\mathrm{\Psi }.$$
(8)
The visit of the inner structure of the atom can be quite interesting, for the fact that atom radius remains strictly constant during the universe expansion .