text:
# Scaling issues in ensemble implementations of the Deutsch-Jozsa algorithm

thanks: Permanent Address: Department of Physics, Guru Nanak Dev University, Amritsar 143005, INDIAthanks: Present Address: Department of Physics, Bucknell University, Lewisburg, PA 17837, USA

## Abstract

We discuss the ensemble version of the Deutsch-Jozsa (DJ) algorithm which attempts to provide a “scalable” implementation on an expectation-value NMR quantum computer. We show that this ensemble implementation of the DJ algorithm is at best as efficient as the classical random algorithm. As soon as any attempt is made to classify all possible functions with certainty, the implementation requires an exponentially large number of molecules. The discrepancies arise out of the interpretation of mixed state density matrices.

Conventional NMR implementations of quantum computing algorithms require the preparation of pseudopure states Knill et al. (1997); Chuang et al. (1998); Cory et al. (1998); Schulman and Vazirani (1999); Boykin et al. (2002). There have been a few proposals to efficiently implement the Deutsch Jozsa (DJ) algorithm on an NMR quantum information processor using highly mixed states Woodward and Brüschweiler (2000); Myers et al. (2002). The basic idea in these ensemble schemes is to avoid pseudopure state preparation which would require exponential resources and instead work with highly mixed states close to thermal equilibrium. These schemes need to be carefuly examined for their ‘quantum character’ and their efficiencies compared to classical random algorithms. In this paper we show that for the DJ problem, a parallel can be drawn between these ensemble implementations and classical random algorithms.

We begin with a brief recapitulation of the DJ problem. Consider the set of functions $`f:\{0,1\}^n\{0,1\}`$. If all $`2^n`$ inputs map to the same output then the function is ‘constant’ and if half the outputs map to $`0`$ while the other half to $`1`$ then the function is ‘balanced’. Functions which are neither ‘constant’ nor ‘balanced’ are not considered here. The task here is to determine the constant or balanced nature of a given function. Given an oracle which evaluates $`f(x)`$ at an input $`x\{0,1\}^n`$, no deterministic classical algorithm can carry out such a classification with certainty without using at least $`2^{n1}+1`$ invocations of the oracle. The quantum DJ algorithm on the other hand accomplishes the classification task by invoking the oracle only once Cleve et al. (1998); Deutsch and Jozsa (1992). This it does by using a quantum oracle defined through the unitary transformation on an $`n`$ qubit argument $`|x`$ and a one qubit target $`|y`$,

$$|x|y\stackrel{\widehat{U}_f}{}|x|yf(x).$$
(1)
If a query to the oracle is assumed to come at a unit cost then the quantum DJ algorithm provides an exponential speedup over its classical counterpart. The $`n+1`$ qubits need to be in a pure quantum state for the algorithm to work. In liquid state NMR at room temperatures, which is the most successful implementation of quantum information processing to date, the quantum states of the spins are far from pure. Therefore, to emulate the standard version of the DJ algorithm, one has to prepare the system in a special ‘pseudopure’ state where the ensemble is divided into two parts: a small subset in a given pure state, and the rest acting as a uniform background with no contribution to the signal. However, the present preparation schemes for such states lead to an exponential loss of signal because the subset of spins which can be prepared in a ‘pure’ state decreases exponentially with the number of qubits Warren (1997); Boykin et al. (2002).

Recently, alternative schemes have been proposed to circumvent this difficulty Woodward and Brüschweiler (2000); Myers et al. (2002). These effectively use a computer with $`n+1`$ qubits, where the first $`n`$ qubits are represented by a density matrix $`\frac{1}{2}I`$ (a fully mixed state) and the last qubit is in the pure state $`|0`$. The initial state thus is

$$\rho _{\mathrm{in}}=\frac{1}{2^n}II\mathrm{}I|00|$$
(2)
and rewriting it in the computational basis $`\{|x|0|x=0\mathrm{}x=2^n1\}`$ gives

$$\rho _{\mathrm{in}}=\frac{1}{2^n}\underset{x=0}{\overset{2^n1}{}}|xx||00|$$
(3)
This preparation is followed by the standard quantum oracle query described in Eq. (1), yielding

$$\rho _{\mathrm{out}}=\frac{1}{2^n}\underset{x=0}{\overset{2^n1}{}}|xx||f(x)f(x)|$$
(4)
Before we extract the constant or balanced nature of the function we note that this is not an entangled state. The entanglement is missing because of the special choice of the initial state. As a matter of fact the oracle is capable of generating entanglement and the standard pure state version of the quantum algorithm relies on entanglement Cleve et al. (1998); Deutsch and Jozsa (1992).

The information about the function is contained entirely in the target qubit whose reduced density matrix is

$$\rho _{\mathrm{target}}=\frac{1}{2^n}\underset{x=0}{\overset{2^n1}{}}|f(x)f(x)|.$$
(5)
The expectation value of $`\sigma _z`$ in this state will immediately reveal the ‘constant’ or ‘balanced’ nature of the function $`f`$, with the result:

$$\sigma _z=\{\begin{array}{cc}\pm 1\hfill & \mathrm{Constant}f\hfill \\ 0\hfill & \mathrm{Balanced}f\hfill \end{array}$$
(6)
By actually carrying out such a measurement the function can be classified with a single invocation of the quantum oracle without the associated problems of preparing pure or pseudo-pure states. This scheme is particularly suitable for implementation on an NMR quantum information processor where the thermal equilibrium state can be easily transformed into the maximally mixed state of Eq. (2) and expectation value measurements are natural.

For comparison we describe a classical scenario which in essence mimics the ‘quantum’ scheme described above. Instead of the NMR qubits, consider classical bit strings of length $`n+1`$. Further, assume that we have $`2^n`$ such strings and each string is in a different state for the first $`n`$ bits, thereby providing representation to all possible states of the first $`n`$-bits. The $`(n+1)`$th bit in each string is set to ‘0’ and acts as the target bit. Now the application of the classical oracle ($`xf(x)`$ with $`f(x)`$ appearing on the target) to all the copies will yield the function values at all the $`2^n`$ input points and this value will be stored in the target bit in each copy. The constant or balanced nature of the function can then be obtained by adding these values by appropriate gates. The values add to $`2^n`$ or $`0`$ for a constant function and to $`2^n/2`$ for a balanced function. This is analogous to the case of the expectation value quantum algorithm using maximally mixed states. In our view this scheme is thus fully classical, using separable states at all stages and camouflaged in the language of quantum mechanics. The exponential resource is explicit in the classical situation while it is hidden in the definition of the input density matrix of Eq. (2) for the ensemble quantum case. The deceptively simple fact that one has effectively prepared the state of the first $`n`$ qubits in the density matrix $`\frac{1}{2^n}\left(III\mathrm{}I\right)`$ requires that we have at least $`2^n`$ molecules!

Every density operator can be viewed as an ensemble of pure states occurring with certain probabilities. The existence of a decomposition

$$\rho =\underset{i}{}p_i|\psi _i\psi _i|\mathrm{with}p_i0,p_i=1$$
(7)
implies an ensemble $`\{p_i,|\psi _i\}`$ for $`\rho `$, where the state $`|\psi _i`$ occurs with probability $`p_i`$. It is to be noted that the states $`|\psi _i`$ need not be orthogonal and the decomposition given above is not unique for mixed states. However, whatever can be determined from $`\rho `$ can be consistently thought of as deriving from any one of the ensembles. The ensemble scheme that culminates in Eq. (6) is illustrated by considering the situation in which each ensemble member is taken to be in one of the computational basis (pure) states, $`\left\{|x|0|x=0\mathrm{}2^n1\right\}`$, each of which occurs with probability $`p_x=1/2^n.`$ This corresponds to the input density operator of Eq. (3). Now view the scheme as it is applied on each member of the ensemble, which means that each member of the ensemble independently computes the function on the input state of that member and the output appears on the $`(n+1)`$th qubit of each member molecule. The measurement of the average of $`\sigma _z`$ then ostensibly reveals the constant or balanced nature of the function as described in Eq. (6).

Now imagine that one chooses to work with a fixed number of molecules $`M`$. As the number of qubits $`n`$ increases, soon one will reach a stage when $`M<2^n/2`$. In this case more than half the computational basis states cannot find representation in the ensemble because there are simply not enough molecules! Therefore there are always balanced functions which will have the same value over these $`M`$ states and will get classified as constant, despite assuming the best situation, where all the molecules are assumed to be in different states. We will see later that the actual scheme is even more inefficient because all the molecules cannot be in different computational basis states. Thus, for the scheme to work for all functions one needs at least $`2^n/2`$ molecules in the ensemble, a number which grows exponentially with $`n`$. In other words, when the number of molecules is smaller than $`2^n`$, there is no way one can prepare the input density matrix of Eq. (2). We will return to this point later.

How many balanced functions escape classification for a given $`M`$? Assume $`N=2^n`$ is the input set size. The number of constant functions is $`2`$ which is independent of $`n`$ while the number of balanced function is $`{}_{}{}^{N}C_{N/2}^{}`$. If the scheme is used with a number of molecules $`1M\frac{N}{2}`$, the balanced functions which escape classification are the ones which have same value ($`0`$ or $`1`$) for the first $`M`$ inputs. The number of functions which have the value $`1`$ ($`0`$) for the first $`M`$ inputs is the same as the number of ways one can distribute the remaining $`N/2M`$ $`1`$’s ($`0`$’s) on $`NM`$ inputs, giving:

$$\begin{array}{c}\hfill \mathrm{No}.\mathrm{of}\mathrm{Balanced}\mathrm{Fns}.\\ \hfill \mathrm{Classified}\mathrm{as}\mathrm{Constant}\end{array}=\mathrm{\hspace{0.17em}2}\left(\genfrac{}{}{0pt}{}{NM}{N/2M}\right)$$
(8)
Dividing this by the total number of balanced functions gives the fraction of balanced functions for which the schemes fails

$$\mathrm{Failure}\mathrm{Fraction}=2\left(\genfrac{}{}{0pt}{}{NM}{N/2M}\right)/\left(\genfrac{}{}{0pt}{}{N}{N/2}\right)$$
(9)
This fraction diminishes quite fast as $`M`$ increases from $`0`$ toward $`\frac{N}{2}`$, increasing the efficiency of the algorithm. The fact that for most cases the scheme will work with a relatively small number of molecules has nothing to do with quantum mechanics. The classical randomized algorithm too will work to the same extent. In fact the above counting is valid for the classical algorithm as well. It is well known that there is an efficient randomized classical algorithm for the DJ problem and we conclude that the expectation value ensemble scheme in the best case is equivalent to it Brazier and Plenio (2003).

Even if we use a molecule number $`M>2^n`$ how safely can we say that the initial maximally mixed state has been realized and we are able to classify all functions with certainty? It is possible that, in a given experimental implementation of the algorithm, all ensemble members are in the same initial state $`|x^{}|0.`$ In this event, the algorithm only evaluates $`f(x^{})`$ and the measurement outcome will be that for a constant function regardless of the nature of $`f`$. Thus, in contrast to the conventional Deutsch-Jozsa algorithm, there can be no way of determining the function type with certainty. In this sense, this ensemble algorithm for solving the Deutsch-Jozsa problem is not deterministic and must be compared to probabilistic classical algorithms. We shall consider the probability with which each correctly determines the function type and show that, regardless of the ensemble size, the standard probabilistic classical algorithm is superior to this ensemble quantum algorithm.

Suppose that the ensemble consists of $`M`$ identical, independent $`n+1`$ qubit molecules. Each member of the ensemble will be subject to the unitary of Eq. (1), which can be re-expressed as

$$\widehat{U}_f=\widehat{P}_0(f)I+\widehat{P}_1(f)\sigma _x$$
(10)
where

$`\widehat{P}_0(f)`$ $`:=`$ $`{\displaystyle \underset{x:f(x)=0}{}}|xx|`$
$`\widehat{P}_1(f)`$ $`:=`$ $`{\displaystyle \underset{x:f(x)=1}{}}|xx|`$ (11)
project onto subspaces of $`n`$ qubit argument while $`I`$ and $`\sigma _x`$ act on the target qubit. For the target qubit, $`\sigma _z`$ is approximated by

$$\overline{z}:=\frac{1}{M}\underset{j=1}{\overset{M}{}}z_j$$
(12)
where $`z_j=\pm 1`$ are the outcomes of projective measurement ($`z_j=1`$ corresponding to $`|00|`$ and $`z_j=1`$ to $`|11|`$) on the target qubit for individual ensemble members.

The only assumption that we make about the ensemble members’ initial states is that they occur with probabilities described by the density operator of Eq. (3). Then, for any ensemble member $`j`$,

$`\mathrm{Prob}(z_j=+1|f)`$ $`=`$ $`\text{Tr}_{\mathrm{arg}}\left(\widehat{P}_0(f)\rho _{\mathrm{in}}\right)`$ (13a)
$`\text{Prob}(z_j=1|f)`$ $`=`$ $`\text{Tr}_{\mathrm{arg}}\left(\widehat{P}_1(f)\rho _{\mathrm{in}}\right)`$ (13b)
where measurements are performed immediately after algorithm unitaries and the traces are taken over the argument register only. Note that each constant function yields one measurement outcome with certainty: $`z_j=+1`$ for $`f=0`$ and $`z_j=1`$ for $`f=1.`$ Thus $`\overline{z}=+1`$ for $`f=0`$ and $`\overline{z}=1`$ for $`f=0`$. Whenever $`\overline{z}`$ departs from $`\pm 1`$ it is clear that $`f`$ is balanced. However, the extent to which such a departure is noticeable depends on the available measurement resolution, which can be expressed in terms of outcomes of function register measurements on individual ensemble members. Suppose that it is possible to distinguish two ensemble averages only when they differ in $`R`$ (out of $`M`$) or more individual measurement outcomes. Then we regard two ensemble averages as distinct provided that $`\left|\overline{z}\overline{z}^{}\right|R/M.`$ This motivates the following protocol for deciding the algorithm outcome:

$`\overline{z}1R/M`$ $``$ $`f=0`$ (14a)
$`\overline{z}1+R/M`$ $``$ $`f=1`$ (14b)
$`1R/M>\overline{z}>1+R/M`$ $``$ $`f\text{balanced}.`$ (14c)
The issue is to determine the probability with which this protocol will correctly identify the function type. Constant functions will always be identified correctly and we need only to find the probability that a balanced function will give $`\overline{z}1R/M`$ or $`\overline{z}1+R/M`$. These are the probabilities that a balanced function will return $`z_j=1`$ or $`z_j=+1`$ at most $`R1`$ times respectively. For a balanced function $`\text{Prob}(z_j=+1)=\text{Prob}(z_j=1)=1/2.`$ The probability that we incorrectly declare a balanced function to be constant is

$`p_{\mathrm{fail}}`$ $`=`$ $`2{\displaystyle \underset{k=0}{\overset{R1}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{M}{k}}\right)\left({\displaystyle \frac{1}{2}}\right)^k\left({\displaystyle \frac{1}{2}}\right)^{Mk}`$ (15)
$`=`$ $`{\displaystyle \frac{1}{2^{M1}}}{\displaystyle \underset{k=0}{\overset{R1}{}}}\left({\displaystyle \genfrac{}{}{0pt}{}{M}{k}}\right).`$
In the best conceivable case $`R=1,`$ giving $`p_{\mathrm{fail}}=1/2^{M1}.`$

To account for the spatial resources offered by the ensemble we consider the application of $`\widehat{U}_f`$ to the ensemble containing $`M`$ members as equivalent to $`M`$ oracle calls. We must then compare this ensemble algorithm to a classical random algorithm that uses $`M`$ oracle calls. In the classical random algorithm one begins by choosing $`x_1`$ randomly and evaluating $`f(x_1)`$. The next step is to choose $`x_2x_1`$, evaluate $`f(x_2)`$ and compare the result to $`f(x_1).`$ If the two differ then $`f`$ is balanced. If not pick $`x_3`$ which differs from both $`x_1`$ and $`x_2,`$ and compare $`f(x_3)`$ to $`f(x_2)`$ and $`f(x_1)`$, etc…. The algorithm terminates when $`f`$ returns different outcomes or has been evaluated on $`2^n/2+1`$ different inputs. The classical random algorithm never misidentifies a constant function and identifies a balanced function $`f`$ as constant only when $`f(x_1)=f(x_2)=\mathrm{}=f(x_M)`$. The probability of failure is the probability with which this occurs. The outcome $`f(x_1)=0`$ occurs with probability $`(N/2)/N`$. Given $`f(x_1)=0`$, $`f(x_2)=0`$ occurs with probability $`(N/21)/(N1)`$. Continuing, the probability that $`x_1,\mathrm{},x_M`$ are all such that $`f(x_k)=0`$ is:

$`p_{\mathrm{classical}}^0`$ $`=`$ $`{\displaystyle \frac{N/2}{N}}{\displaystyle \frac{N/21}{N1}}\mathrm{}{\displaystyle \frac{N/2M+1}{NM+1}}`$ (16)
$`=`$ $`\left({\displaystyle \genfrac{}{}{0pt}{}{N/2}{M}}\right)/\left({\displaystyle \genfrac{}{}{0pt}{}{N}{M}}\right)`$
Similarly the probability $`p_{\mathrm{classical}}^1`$ that $`x_1,\mathrm{},x_M`$ are all such that $`f(x_k)=1`$ can be computed and it turns out to be same as $`p_{\mathrm{classical}}^0`$. Thus the probability of failure is:

$$p_{\mathrm{fail}}^{\mathrm{classical}}=2\left(\genfrac{}{}{0pt}{}{N/2}{M}\right)/\left(\genfrac{}{}{0pt}{}{N}{M}\right).$$
(17)
As expected, this turns out to be same as the failure fraction given in Eq. (9) for the quantum ensemble version with the unrealistic assumption that every molecule is in a different computational basis state. However, note that for $`N/2>k>0`$, $`(N/2k)/(Nk)<1/2,`$ which implies that

$$p_{\mathrm{fail}}^{\mathrm{classical}}<2\left(\frac{1}{2}\right)^Mp_{\mathrm{fail}}.$$
(18)
Thus the probability of failure for the classical random algorithm is strictly less than that of the ensemble quantum version discussed here.

The central issue is therefore one of interpreting density matrices. What is relevant here are the inferences that can be drawn from measurement outcomes on quantum systems whose states are described by density matrices. In all cases the density matrix merely provides the probability distribution for outcomes of various measurements. The accuracy with which such a distribution is realized improves with an increasingly large ensemble. Imagine a single quantum system which is handed over to us with no information about it. What quantum state or density matrix will we be able to assign to it? To express our complete lack of information about this system we have to assign equal weightage to all possible outcomes in all bases and therefore a density matrix proportional to identity is the best choice. In this extreme case, measurement yields one of all possible outcomes and one cannot reliably infer anything from this. The density operator merely reflects our lack of knowledge in the state of the system. It is only when measurements are performed on many copies described by the same density operator that outcomes or more precisely, the average outcome, carry any meaningful information. This is no more than standard statistical sampling and for an $`n`$ qubit density matrix proportional to the identity, the variance scales as $`2^{2n}`$, indicating that one typically needs $`O(2^n)`$ samples (measurement outcomes) to make sensible inferences from measurements. Accordingly the ensemble size would have to scale as $`O(2^n)`$ before we can consider this density operator to have been realized accurately, at least in terms of measurement outcome averages. Here the ensemble begins to appear as a collection of quantum systems with states described according to the density matrix. It should be noted that, for this version of the DJ algorithm, the situation is less dire since inferences are made from measurements on the target qubit alone. Hence $`n`$ does not appear in the failure probability in Eqn. (15). However, as clear from Eqn. (18) a classical random algorithm does the task better.

These ensemble computing ideas might work for other algorithms and give a genuine speed up over classical or classical random algorithms. One possibility is efficient simulation of quantum systems Knill and Laflamme (1998) and our result does not pertain to this. It is worthwhile to explore the exact implications of this model which will be taken up elsewhere.

###### Acknowledgements.

Numerous stimulating discussions with R. B. Griffiths and Kavita Dorai are acknowledged. The research effort is funded by the NSF Grant No. 0139974.


format: text

text:
# Medium influence of the nucleon nucleon cross section on the fragmentation

## 1 ACKNOWLEDGMENT

We thank Prof.Bao-An Li for helpful discussions.
This work is supported by the Major State Basic Research Development Program in China Under Contract No.G2000077400, ”100-person project” of the Chinese Academy of Sciences, the National Natural Science Foundation of China under Grants Nos.10175080, 10004012,10175082 and The CAS Knowledge Innovation Project No. KJCX2-SW-N02

## Figure captions

An impact parameter average value of the multiplicity of intermediate mass fragments,$`<N_{imf}>_b`$ as a function of the beam energy for $`U_0^{sym}`$, $`U_1^{sym}`$ and $`U_2^{sym}`$ with $`\sigma ^{iso}`$ and $`\sigma ^{noiso}`$ for the mass symmetry system 76Kr + 76Kr and asymmetry mass system 112Sn + 40Ca(see text).

The time evolution of the $`N_{imf}`$ with MDI and NOMDI for four cases (see text) and systems $`{}_{}{}^{76}Kr+^{76}Kr`$ at E=150 MeV/nucleon and b= 4.0 fm.

The time evolution of the ratio of nuclear density to normal density $`\frac{\rho (\sigma ^{iso},\alpha )}{\rho _0}`$ for different cases(see text).


text:
# Enhanced thermoelectric properties of NaxCoO2 whisker crystals

thanks: Corresponding author: Prof. Hisao Yamauchi
Materials and Structures Laboratory, Tokyo Institute of Technology
4259 Nagatsuta, Midori-ku, Yokohama 226-8503, Japan
Phone: +81-45-924-5315, Fax: +81-45-924-5365
Electronic mail: yamauchi@msl.titech.ac.jp

## Abstract

Single-crystalline whiskers of thermoelectric cobalt oxide, NaxCoO2, have been grown by an unconventional method from potassium-containing compositions, and their transport properties, and chemical compositions were determined. The growth mechanism was analyzed and interpreted by means of thermogravimetric analysis. At 300 K, electrical resistivity $`\rho `$ and thermoelectric power $`S`$ of the whisker are 102 $`\mu \mathrm{\Omega }`$ cm and 83 $`\mu `$V/K, respectively. The power-factor, $`S^2/\rho `$, thus is enhanced up to $``$ 68 $`\mu `$W/K2 cm.

preprint: to appear in Appl. Phys. Lett.

Thermoelectric materials that can convert heat flow directly to electric flow/current and vice versa are considered to play a key role in “clean” energy utility of the next generation Mahan97a . However, requirements for high-efficiency thermoelectric materials are not easily satisfied. The thermoelectric efficiency is evaluated by the following figure-of-merit: $`Z=S^2/\rho \kappa `$, where $`S`$, $`\rho `$, and $`\kappa `$ denote thermoelectric power, electrical resistivity, and thermal conductivity, respectively. Therefore, large $`S`$ values and small $`\rho `$ and $`\kappa `$ values are simultaneously required to realize high thermoelectric efficiency. For the last few years thermoelectric oxides have been widely investigated owing to their high chemical stability at temperatures higher than $``$800 K. The sodium-cobalt oxide, NaxCoO2, is one of the most promising candidates among the thermoelectric oxides. Terasaki et al. Terasaki97a  found that this compound (originally expressed as NaCo2O4) exhibits simultaneously high thermoelectric power and low resistivity ($`S`$ $``$ 100 $`\mu `$V/K and $`\rho `$ $``$ 200 $`\mu \mathrm{\Omega }`$ cm at 300 K). However, due to its two-dimensional crystal structure, the compound has large electromagnetic anisotropy Terasaki97a , and thus polycrystalline samples always show higher resistivity values ($`\rho `$ $``$ 3 m$`\mathrm{\Omega }`$ cm at 300 K Terasaki00a ). For reducing the resistivity and thereby enhancing the $`Z`$ value of NaxCoO2, it is indispensable to develop a method to prepare textured or single-crystalline materials.

Whisker-type single crystals are of significant importance due to the following reasons: (1) the one-dimensional shape of whiskers should be favorable for some synthetic textures, and (2) whiskers may be convenient for some tiny-dimensional applications such as natural “thermoelectric wires” taking advantage of their characteristic shape. It was shown that whiskers of Bi2Sr2Co2Oy Funahashi01a ; Funahashi02a  exhibit significantly enhanced thermoelectric properties as compared with polycrystalline samples. Moreover, whiskers are known to be high-quality single crystals with low defect densities. Therefore, precise investigations on the crystal structure and physical properties of NaxCoO2 can be made using whisker crystals. In this letter we report thermoelectric properties of NaxCoO2 whiskers prepared using an original growth technique.

It is known that NaxCoO2 forms plate-like crystals when NaCl is used as a flux. To obtain NaxCoO2 whiskers we used potassium-containing compositions. Dry powders of Na2CO3, K2CO3, and Co3O4 were mixed with cation ratios, Na : K : Co = 0.525 : 0.225 : 1.00, and calcined at 750C in air employing a “rapid heat-up” technique Motohashi01a  to avoid evaporation of the alkali metals. The calcined powder was ground, pressed into pellets, and fired at 900C for 12 h in air. During the firing, partial melting apparently occurred, and whiskers were found to form at the surface of the pellets (Fig. 4).

Figure 4 shows an x-ray powder diffraction pattern (measured with Cu$`K_{\alpha 1}`$ radiation) for pulverized whiskers. All the peaks in the pattern can be indexed according to the $`\gamma `$-NaxCoO2 phase Fouassier73a . Crystal structure analysis using a four-circle diffractometer and synchrotron radiation revealed that the whiskers possess a hexagonal structure with space group, $`P6_3/mmc`$ (No. 194) Hanashima02a , being in agreement with the neutron diffraction data for a polycrystalline sample of Na0.74CoO2 Balsys96a . The grown whiskers are up to 1.6 mm in length, 15 - 40 $`\mu `$m in width, and 1.5 - 4.0 $`\mu `$m in thickness: the shape is strongly one-dimensional in noticeable contrast to plate-like crystals grown by a NaCl-flux method Terasaki97a ; Fujita01a . The growth direction of the whiskers was found to be along the ab plane, while the smallest dimension is along the $`c`$ axis.

Electrical resistivity ($`\rho `$) and thermoelectric power ($`S`$) of the whiskers were measured: $`\rho `$ was measured using a four-point-probe apparatus (Quantum Design; PPMS) with a dc current of 5 mA along the growth (i.e. in-plane) direction, and $`S`$ using a steady-state technique with a typical temperature gradient of 0.5 K/cm. Figure refr-and-S(a) shows the dependence of $`\rho `$ on temperature for a NaxCoO2 whisker. The $`\rho `$ value decreases with decreasing temperature, exhibiting typical metallic behavior (i.e. $`d\rho /dT>0`$) in the whole temperature range. This is in accordance with the results reported for single crystals in previous works Terasaki97a ; Fujita01a . An anomaly in resistivity is observed at $`T_t`$ $``$ 290 K, which is considered to be related to a first-order transition Tojo02a , while such an anomaly has not been observed for single-crystalline samples of NaxCoO2 Terasaki97a ; Fujita01a . The absolute value of resistivity for the present whiskers is by a factor of two smaller than that for single crystals of NaxCoO2 grown with the NaCl flux method (cf. $`\rho `$ $``$ 200 $`\mu \mathrm{\Omega }`$ cm at 300 K) Terasaki97a . Figure 4(b) shows the dependence of $`S`$ on temperature for a NaxCoO2 whisker. Both the absolute value and the dependence on temperature of $`S`$ are in good agreement with those previously reported for single-crystalline and polycrystalline samples Terasaki97a ; Terasaki00a . The values of $`\rho `$ and $`S`$ at 300 K for the present whisker crystal are 102 $`\mu \mathrm{\Omega }`$ cm and 83 $`\mu `$V/K, respectively. The power-factor ($`S^2/\rho `$) is then calculated at $``$ 68 $`\mu `$W/K2 cm, being as large as that for a crystal grown by the NaCl-flux method Terasaki97a  and by a factor of $``$ 30 larger than that for polycrystalline samples Terasaki00a .

In order to determine the chemical composition of grown whiskers, we performed ICP-AES (Seiko Instruments: SPS-1500VR) and SEM-EDS (Hitachi: S-4500) analyses. In Table 1, the actual cation composition of the whiskers as determined by the two independent analysis methods is given together with those of the calcined powder and the fired pellets. Neither the whiskers nor the pellets contained potassium, whereas the actual composition of the calcined powder was nearly the same as the nominal one. This result indicates that fast evaporation of potassium occurs during firing (at $``$ 900C), which is apparently essential for the growth of whiskers. The chemical formula determined for the whisker (in terms of cation composition) is Na0.5CoO2, being highly sodium deficient and free from potassium.

To investigate the whisker formation mechanism, thermogravimetric (TG) measurements were performed. A powder mixture of raw materials of 4.295 mg was heated with a rate of 3.2C/min in a thermobalance (Perkin Elmer: Pyris 1) with an isothermal heating period of 12 h at 900C. The TG curve is shown in Fig. 4. Since K2CO3 is highly hygroscopic, the first weight loss about 100C is attributed to the evaporation of water. The second weight loss is seen starting about 500C. This is most likely due to the decarbonation process, though the total weight loss in this region ($``$ 6.4 %) is a little too large to be explained with decarbonation (5.8 %) only. Therefore it is considered that a small amount of alkali metal(s) evaporates as well. The two weight drops about 100C and 500C were also seen for the potassium-free sample of Na0.55CoO2 Motohashi01a . After the large weight drop due to decarbonation, a temperature window was found to appear about the calcination temperature (650-830C) where the TG curve is relatively flat. This should correspond to the stability region of the NaxKyCoO2 phase. In fact, x-ray diffraction measurements showed that the powder product after calcination at 750C is of single phase of hexagonal NaxKyCoO2 phase, whereas chemical analyses for the cation composition revealed essentially no deviation from the nominal potassium content (see Table 1). Upon increasing the temperature up to 900C an additional weight loss is seen indicating that the fast evaporation of potassium from the NaxKyCoO2 phase starts to occur in this temperature range. In the inset of Fig. 4 the normalized weight is plotted against time. During the isothermal heating period of 12 hours at 900C the weight loss continues until potassium completely gets depleted from the sample: the evaporation rate of potassium from NaxKyCoO2 was estimated at $``$ 6.5 mol%/h. The TG result is in good agreement with the ICP analysis data (Table 1), according to which no potassium was found to remain in both the pellets and the whiskers. Taking into account the results of TG and ICP-AES analyses, it is most likely that a large amount of potassium evaporation induces partial melting at the surface of the pellet, and the partially-melted surface plays an essential role in the whisker growth.

In summary, we have developed an unconventional method to grow whisker single crystals of NaxCoO2. Using potassium-containing compositions, whisker crystals of NaxCoO2 are found to form after firing at 900C. At 300 K, $`\rho `$ and $`S`$ of the whiskers are 102 $`\mu \mathrm{\Omega }`$ cm and 83 $`\mu `$V/K, respectively. Calculated power factor ($`S^2/\rho `$) at this temperature shows a record-high value of $``$ 68 $`\mu `$W/K2 cm for the sodium-cobalt-oxide system. From ICP-AES and SEM-EDS analyses, it has been revealed that the grown whiskers are highly sodium deficient ($`x`$ $``$ 0.50) and free of potassium. Increasing sodium content up to 0.75 can probably lead to even higher thermoelectric power ($`S^2/\rho `$) values, therefore preparation of whiskers with higher sodium contents are on the way in our laboratory.

We are grateful to Prof. I. Terasaki of Waseda University for his contribution in thermoelectric power measurements. The present work was supported by Grants-in-aid for Scientific Research (Contract Nos. 11305002 and 14038219) from the Ministry of Education, Culture, Sports, Science and Technology of Japan.


text:
# The tidal perturbation of the low-mass globular cluster NGC 5466

## 1. Introduction

NGC 5466 is a metal-poor halo cluster located at $`d=16.3`$ kpc from the Sun (Cassisi et al. 2001) and $`z=15.6`$ kpc above the Galactic plane (see Fig.1). It is one of the less luminous clusters in the globular cluster system of our Galaxy. It has a low central velocity dispersion of only 1.5 km s-1 (Pryor & Meylan 1993) and a relatively large core radius ($`r_c=6.4`$ pc, Lehmann & Scholz 1997). This implies a dynamical mass of about $`5\times 10^4M_{}`$ (Pryor et al. 1991). The internal properties of NGC 5466 are thus similar to - though less extreme than - those of the low-mass halo cluster Palomar 5, which was recently found to have long tails of tidal debris (Odenkirchen et al. 2001, 2003) that reveal ongoing tidal disruption.
Existing measurements of the proper motion of NGC 5466 (see Odenkirchen et al. 1997) allow us to assess its local Galactic orbit. We took the available kinematic data for NGC 5466 and a model of the Galactic potential developed by Allen & Santillan (1991) and integrated the equations of motion to $`\pm 100`$ Myr from present. The resulting orbit is plotted in Fig.1. It demonstrates that NGC 5466 is on a highly eccentric orbit with $`R_{max}>40`$ kpc, but has crossed the Galactic disk about 50 Myr ago at a distance of about 8 kpc from the Galactic center. This suggests that NGC 5466 has recently undergone a tidal shock, which - due to the low mass and low concentration of this cluster - might lead to a significant loss of stars.

## 2. Analysis of APM data

We analysed the spatial distribution of stars in the outer parts of NGC 5466 using data from the APM catalogue111data available at www.ast.cam.ac.uk/~apmcat/. This catalogue is based on scans of the first Palomar Sky Survey (POSS-I) O and E plates. We extracted a $`3^{}\times 3^{}`$ field centered on the cluster and selected all sources classified as “starlike”. To weed out stars that are most likely not associated with the cluster we restricted the sample to the color-magnitude box defined by $`18.3R20.0`$ and $`0.7BR1.9`$. This selects stars on the upper main sequence and near the main-sequence turn-off of NGC 5466. Since the photometry obtained from the POSS-I plates is of moderate accuracy, the selection necessarily is rather coarse.
Fig.2 shows the distribution of the selected cluster member candidates on the sky. While the inner part of the cluster appears spherical, there is a stellar halo around it, which has a non-spherical shape. This halo appears elongated towards the south-west and the north-east, i.e., stars are seen preferentially near the directions towards the Galactic center (see dashed line in Fig.1) and anticenter. The orientation of this cluster halo is further analysed in Fig.3, where we plot star counts in the range $`9^{}r18^{}`$ as a function of position angle. The position angle of the major-axis of the cluster’s halo is about +30 while the position angle of the Galactic center-anticenter line is +43 . This provides first evidence for a possible tidal perturbation of the cluster and eventual mass loss because the Galactic tidal field will drag stars away from the cluster along the center-anticenter line (see, e.g., the example of Pal 5 in Odenkirchen et al. 2003). However, the APM data do not allow to trace cluster members to more than $`20^{}`$ from the center of NGC 5466. Thus deeper and more accurate photometry is needed to explore the distribution of cluster stars in more detail.

## 3. Comparison to simulated tidal debris

Fig.4 shows a prediction of the spread of tidal debris from NGC 5466. We simulated tidal mass loss from the cluster by releasing particles from the Lagrange points located on the Galactic centre - anticenter line in time steps of 10 Myr. The particles were assumed to have peculiar velocities of size 1.5 km s-1 (equal to the measured velocity dispersion of the cluster). The dots show the present location of particles released $`200`$ Myrs ago. Stars lost during the last 50 Myrs, i.e., after the cluster’s recent disk crossing, are highlighted as fat dots. The distribution of the fat dots is similar to the observed halo of cluster member candidates around NGC 5466. The orientation of the observed halo is compatible with a transition from the center-anticenter line towards alignment with the cluster’s orbit. Deeper and cleaner data might allow to resolve the sharp bending near the Lagrange points suggested by the simulation.

## References

Allen C., & Santillan A. 1991, Rev. Mex. Astron. Astrofis., 22, 255

Cassisi S., De Santis R., & Piersimoni A. M. 2001, MNRAS, 326, 342

Lehmann I., & Scholz R.D. 1997, A&A320, 776

Odenkirchen M., Brosche P., Geffert M., & Tucholke H.-J. 1997, NewA 2, 477

Odenkirchen M., Grebel E.K., Rockosi C.M., et al. 2001, ApJ, 548, L165

Odenkirchen M., Grebel E.K., Dehnen W., Rix H.W., et al. 2003, AJ, accepted

Pryor C., Mc Clure R.D., Fletcher J.M., & Hesser J.E. 1991, AJ, 102, 1026

Pryor C., & Meylan G. 1993, in ASP Conf. Ser. Vol. 50, Structure and Dynamics of Globular Clusters, eds. S.G. Djorgovski & G. Meylan, p.357


text:
# Lattice QCD Production on Commodity Clusters at Fermilab

## I Introduction

Large scale QCD Monte Carlo calculations have typically been performed on either commercial supercomputers or specially built massively parallel computers. Commodity clusters equipped with high performance networking equipment present an attractive alternative, achieving superior performance to price ratios and offering clear upgrade paths. The U.S. Department of Energy, through the SciDAC (Scientific Discovery through Advanced Computing) program, is supporting the investigation of commodity clusters as well as purpose built machines for lattice gauge computations.

Lattice QCD codes are memory bandwidth, floating point, and network intensive. Successful design of clusters to run these codes requires knowledge of the bottlenecks which control performance. In this paper we examine a number of aspects of performance, observed while running the MILC lattice QCD code. For single systems, we discuss the effects of the various memory architectures. We examine optimizations such as data layout and SSE-assisted matrix algebra, and consider SMP behavior. We discuss scaling of the MILC code on Myrinet-connected clusters. Using a modified version of GM which allows control of bandwidth and latency, we examine the sensitivity of the code to network performance.

## II The SciDAC Lattice Gauge Computing Initiative

In 2001, a collaboration which includes most U.S. lattice theorists, applied for funding to the DOE through the SciDAC initiative scidac . Support for fiscal years FY2002 through FY2004 was awarded, with funding primarily for software development, but also for prototype clusters. The SciDAC Lattice Gauge Computing Project scidac\_lqcd  pursues hardware investigations along two paths: special purpose computers (QCDOC qcdoc ) under development by Columbia University, and commodity clusters. A primary goal of the project is to design and implement the software infrastructure necessary to allow legacy and new lattice gauge codes to run on both types of hardware. At the present time it is anticipated that purpose built machines like the QCDOC will achieve the best performance to price ratio for several years, with clusters taking the lead in the following years. Ideally the U.S. community will develop and have access to several facilities which will each provide sustained computing performance on large problems of order one to ten teraflops.

SciDAC software scidac\_soft  will include low level application program interfaces for communications (QMP) and for linear algebra (QLA). QMP will offer features similar to a small subset of MPI, with less overhead. It will run transparently over the QCDOC mesh, Myrinet GM, gigabit ethernet meshes, and MPI. QLA will perform linear algebra math kernels which will run on a single node; it will be lattice aware, offering operations at single sites and across the entire sublattice held by a single computational node. QLA will have low level optimizations for some architectures, for example, through the use of the SSE (Streaming SIMD Extension) instructions on x86 CPUs. Lattice wide computations, which involve communications, will be addressed by the QDP API. Lattice input and output, including parallel file I/O when allowed by hardware, will be provided by the QIO API.

## III The Fermilab Clusters

Shown in Table 1 are the parameters of the lattice QCD clusters at Fermilab. All of these systems were purchased with supplemental grants from the Department of Energy, the latter three under the SciDAC program. On all clusters we use the OpenPBS batch software with the Maui scheduler to control user jobs. The MPICH implementation of MPI with the Myricom-supplied driver is used for the communications API; currently we use version 1.2.4..8a, built against GM version 1.6.3.

As is traditional with MPICH, user jobs are launched on the set of nodes assigned by the batch queue system by means of a user script started on the rank 0 node. The user script invokes mpirun, which uses rsh to start the user binary on each of the nodes. Unfortunately, there is a 256 process barrier encountered for large jobs. Because standard rsh uses a pair of privileged TCP ports to connect to the target node, and privileged ports are limited in number, large jobs fail because of exhaustion of these TCP ports. Instead of rsh, use instead use ssh, which can be configured to use non-privileged ports. To avoid burdening users with the task of managing ssh keys, we also configure the ssh daemons to accept RSA host authentication.

## IV Single Node Performance

For several years we have used the performance of the Dirac inverter in the MILC improved staggered code (“su3\_rmd\_symzk1\_asqtad”) to benchmark the performance of systems. Figure 1 is a typical performance graph, which shows the sustained floating point operations delivered by different processors as a function of the lattice size, when the MILC code is run as a single process on a single computer. The performance for lattices smaller than 44 is dominated by the floating point power of the processor. Near 44 on the processors shown, the lattice begins to surpass the L2 cache size (512K bytes), and for larger lattices the memory bandwidth of the computer determines performance. The processor with the fastest front side bus in the graph has the highest performance in main memory, even though it does not have the highest clock speed.

To help assess floating point performance, we have implemented a modified version of McCalpin’s “Streams” benchmark streams . The new benchmark, qcdstream, uses several SU3 math kernels with a variety of memory access patterns qcdstream . Table 2 shows the performance of two of the kernels on systems using 2.0 GHz Xeon processors based on the E7500 chipset (interleaved DDR memory). Also shown is the measured performance of 900 MHz Itanium 2 processors on the same kernel. The Xeon processors ran standard “C” language code, optimized by the GCC compiler. The Itanium 2 processors ran kernels optimized by John Dupuis of Hewlett-Packard. The memory access patterns are as follows: “in-cache” - operands are reused repeatedly; “sequential” - operands are assigned from consecutive locations in pre-allocated memory regions; “strided” - operands are assigned from non-adjacent locations from pre-allocated memory regions as determined by a stride constant significantly larger than a cache line; “mapped” - operands are assigned from pre-allocated memory regions via a randomly assigned indirect mapping.

The memory bandwidth of the commodity systems we have investigated is a function of both the processor and the chipset. For x86-type processors, the highest memory bandwidths currently available are on Pentium 4 systems with 800 MHz front side buses with interleaved DDR chipsets, such as Intel’s i875. Table 3 lists memory bandwidths as measured by the “copy” section of McCalpin’s Stream benchmark.

Two optimizations are very effective in improving the performance of the MILC inverter. First, re-ordering the data structures so that they are “field major” rather than “site major” has the effect of packing SU3 matrices together and achieving more efficient memory bus utilization. This is particularly important on Pentium 4 class processors. These processors have 64 byte cache lines. Single precision SU3 matrices are 72 bytes in length, so two cache lines must be loaded for a single matrix operand. If the extra 56 bytes loaded are not used before they are flushed out of the L2 cache, effective memory bandwidth is greatly reduced. Field major ordering, which essentially places the next matrix to be used immediately after a given matrix in memory, substantially increases performance.

The second optimization is the use of SSE instructions to vectorize the SU3 matrix algebra. Inline versions of the fundamental MILC kernels used in the Dirac inverter have been implemented sse . Because code which utilizes GCC assembler macros is very difficult to debug and maintain, we instead have written these routines using a conventional assembler (NASM nasm ). The assembly code is translated to inline GCC assembler macros by simple parser implemented in Perl.

The effects of these optimizations are shown in Figure 2. We note that in main memory, the field major optimization greatly enhances the effect of the SSE optimization.

In Figure 3 are shown the measured performances of the fully optimized code run on a variety of processors. The top performer is the 533 MHz front side bus version of the Pentium 4. At the time of these measurements, Xeon systems were limited to 400 MHz front side buses. However, 533 MHz FSB Xeons and 800 MHz FSB Pentium 4 systems are now available. Note that Pentium 4 systems have limited I/O buses, with only 32 bit, 33 MHz PCI slots, whereas Xeon systems typically have one or more PCI-X buses (133 MHz, 64 bits).

## V SMP Performance

Commodity systems based on x86 processors are available in either single or dual processor form. Because high performance network adapters, such as Myrinet, are expensive, SMP systems may be more cost effective than single processor systems if codes scale well and the communications bandwidth supplied by a single network interface is sufficient. Unfortunately, as shown in Section IV, lattice codes are often memory bandwidth bound, and dual processor x86 systems generally have the same available bandwidth as single processor systems. For the MILC code, this effect is shown in Figure 4. The scaling of the MILC code on a dual Xeon processor system was found to be about 65% with two independent processes, and about 55% with two cooperative (MPI) processes.

Pentium 4 class processors have the best single node performance for MILC codes. For these processors, only SMP motherboards based on the Xeon processor have fast, wide PCI buses, which are critical for I/O performance. These systems can be operated with only a single processor. However, since the incremental cost of the second processor is a small fraction of the total cost of the systems including the high performance network fabric, we elected to purchase the second processor for each node of our clusters. We note that a number of other applications run by users on our clusters are not as sensitive to memory bandwidth and so benefit from the availability of the second processor.

## VI Cluster Performance

In Figures 5 and 6 we show the performance of the MILC code on the 2.4 GHz Xeon Fermilab cluster as the number of nodes assigned to a parallel run varies from 1 to 128. On these runs, the sublattice assigned to each process was kept constant. The number of dimensions with communications varies with the number of processes. Communications along 1, 2, 3, and 4 directions occurs, respectively, for runs with 2, 4, 8, and 16 or greater nodes.

The performance of the PCI bus on x86 systems varies with chipset. Unfortunately, the best Xeon chipset for memory bandwidth, the Intel i860, has poor I/O performance. In Table 4, we list the burst transfer rates between a Myrinet interface and main memory, as measured by the GM driver. Since the Myrinet “wire rate” is 250 MBytes/sec, I/O performance will be constrained on motherboards with PCI burst transfer rates below this rate. Indeed, the large message size transfer rate asymptotes of i860 and E7500 chipset systems are, respectively, approximately 190 and 225 MBytes/sec, as measured by the Netpipe benchmark netpipe .

To investigate the effect on lattice codes of PCI chipsets, we measured the performance of the Dirac inverter on 32-node runs. Since frequent barrier synchronizations occur during the global sums used in the inverter, by substituting just one slower system for one of the computers, we can closely estimate the performance of a cluster consisting of slower systems. For this study we used the 2.0 GHz Xeon cluster based on the E7500 chipset. We measured performance with three slower system substitutions: a 1.7 GHz i860-based system, a 1.7 GHz E7500-based system, and a 2.26 GHz i850E-based system. The latter computer has only a narrow, slow PCI bus. The results are shown in Figure 7. The i860-based system does not substantially lower performance, whereas the slow, narrow PCI bus of the i850E system has a large effect.

In addition to bandwidth studies, we also investigated the effects of latency on the MILC code. We started with a modified version of the GM driver for Myrinet from D. K. Panda et al panda2 . His group’s modifications add quality of service functionality to GM, restricting the bandwidth available to a given connection by delaying each Myrinet packet. We added the ability to delay the initial packet, rather than subsequent packets, allowing control over the latency of communications. The results of a series of MILC runs are shown in Figure 8, where we have plotted the performance of the entire Dirac inverter (“CONGRAD”) and the underlying Dirac operator (“D-slash”) against the zero-length message latency. D-slash includes only the large message communications required to access a neighboring node’s surface lattice sites. CONGRAD includes D-slash as well as the global sums required to test convergence. D-slash is insensitive to latency; the full CONGRAD, because of the global sums, is much more sensitive.

## VII Future Work at Fermilab

Lattice gauge codes require networks with good bandwidth and low latencies. Commercial solutions, such as Myrinet, SCI and Quadrix, are readily available but at substantial cost. Indeed, the fraction of the total cost of the three Intel-based lattice gauge clusters at Fermilab represented by the Myrinet fabrics has been about 45%. Gigabit ethernet offers a significantly lower priced alternative. Unfortunately, compared to the other networks, switched gigabit ethernet suffers from lower bandwidth, higher latencies (because switches typically are of the “store and forward” type), limited switch sizes, and immature software (low overhead and reduced software latencies require the use of non-TCP/IP protocol stacks such as VIA). Gigabit ethernet meshes, such as studied by Fodor et al fodor , are very inexpensive in comparison to switched networks. However, clusters built with these meshes have rigid configurations, with lattice layouts fixed by wiring. They are also not tolerant of node failures; switched networks handle failures with ease. Finally, although communications with neighbor nodes have low latency with protocols like VIA, non-nearest-neighbor communications require interaction with the operating system of the intervening nodes, giving poor latency performance.

Recently several FPGA (field programmable gate array) manufactures have introduced so-called “platform” gate arrays, which include substantial I/O capabilities and in some cases embedded PowerPC processors. Fermilab has investigated these FPGAs for applications in data acquisition, building prototype PCI interface cards with as many as eight bidirectional high speed serial links. We will study the application of these cards to lattice gauge computing. The long term goal is to assess the feasibility of using these interfaces to build 4-D or higher order meshes. Like gigabit ethernet meshes, FGPA-based meshes would provide high, scalable bandwidth (each link will have peak bandwidth of 2 GByte/second in each direction). Unlike gigabit ethernet meshes, non-nearest-neighbor communications would not incur substantial latency penalties, as routing would be performed within the fabric itself.

During FY2004, Fermilab will purchase an additional cluster. Depending upon final funding and the choice of CPU architecture and network fabric, this cluster will be as large as 256 nodes. The choices available for the processor including Xeons with 533 MHz or higher front side buses, the new Pentium 4 processor (“Prescott”, which includes additional SSE instructions designed for complex arithmetic) with 800 MHz or higher front side bus, the next generation Itanium processor (“Madison”), AMD Opteron, and the PowerPC 970. As discussed in this paper, dual Xeon motherboards suffer from inadequate total memory bandwidth, but have good I/O buses. Pentium 4 motherboards have superior memory bandwidth, but to date have only been available with slow, narrow PCI buses. There is hope that new Pentium 4 motherboards with high performance I/O buses (“PCI Express”) will be available in time. Itanium processors offer superior performance, but only when considerable manpower is expended in hand tuning the code. Both Opteron and PPC970 offer the promise of superior, scalable memory bandwidth and excellent floating point capabilities. We will choose the processor with the best system performance to price ratio. For network fabrics, we will choose between Myrinet, gigabit ethernet over Myrinet, gigabit ethernet mesh, and Infiniband.

###### Acknowledgements.

Work supported by Universities Research Association Inc. under Contract No. DE-AC02-76CH03000 with the United States Department of Energy.


text:
# Login times in e-mail servers are scale-free

## Abstract

The login time distribution in e-mail servers are studied. The distribution seems to be decaying in a scale-free manner with a decay exponent lying between 1.2 and 2.2, varying from server to server. A simple mean-field theory shows that this distribution has a decay exponent 2.

Internet has drastically changed the mode of communication in the last few decades. E-mail has made communication efficient, fast and easy across the globe and helps business, scientific and social interactions. The structure of this communication network is dynamic, where each user is a node and the packet of information trasmitted (an e-mail) constitutes a link between a pair of nodes (users, as sender and recipient). The structure of such a network has been found to be scale-free ebel  and is somewhat similar to network structure of the internet alb1 . E-mail has been looked upon as a social network well , whose degree distribution and mode of information flow has also been studied info . Given the network topology, how the individual nodes respond temporally in terms of their activity has not yet been studied. A very important and immediate example in this respect is a study of lifetimes of different logins in an e-mail server. The fundamental question stands as: How much time does an user spend during a particular login in his/her e-mail account? It can range from very small time-scales (bounded by the time required to log out immediately after logging in) to practically infinity. But is practice, it ranges from the order of seconds to many days, depending on the purpose and necessity of the login.

Analyzing data from stock-market bbm1 , Internet (number of links to a page alb1 , number of pages within a site ah1 ), as well as natural phenomena like earthquake gut , rainfall rain  etc. yielded interesting observations, and our study on the distribution of login times is not an exception. Login time data from different countries are distributed in time according to an universal power law: there are many logins of short durations, whereas, very few extend for very large times. This universal distribution can be explained by a simple mean-field theory.

We collected login time data from e-mail servers of all major research centres in India and a few places in different parts of the world. The computers register login times to the extent of detailed particulars as user, date, login entry time, login exit time and net login time. The registered login times are in units of minutes, as it is the standard time resolution set to measure login times in the machines. Hence, all data collected were in terms of integral minutes. The login time distribution (see Fig 1) turned out to be scale-free, i.e,

$$P(t)t^\alpha $$
(1)
with the decay exponent $`\alpha `$ varying between 1.2 to 2.2 from server to server. Table 1 shows the decay exponents for different e-mail servers.

The data suggest that there are more logins of shorter durations than of longer durations, which is quite natural. Speculations could only predict the decaying nature of the distribution, but not why it is scale-free. The data shows that there is no typical time-scale of a login.

The origin of different login times are due to human behavior to check e-mails frequently between work, while those of larger time scales are due to logins which take place once or twice in a day to check e-mails coming in the previous few hours. The even larger time scales arise due to file-transfer processes, downloads and even due to reasons like ‘carelessness’!

We explain the scale-free nature of the login times by a simple mean field theory. We define $`p`$ as the probability of any user to log out at any time. We, however, consider no interaction between users for simplicity, and effects of heterogeneity and frustration are only realized by assigning a distribution of $`p`$. Hence, we treat this as a single-particle problem, and the cumulative effect in the many-particle situation is realized by assuming a distribution of $`p`$.

As per our assumption, the probability of an user to have not logged out at time $`t1`$ is $`(1p)^{t1}`$ and logs out at time $`t`$ is $`\rho (p,t)=p(1p)^{t1}`$. For a multi-user system, where $`p`$ follows a distribution $`f(p)`$, the probability of logging out at time $`t`$ is

$$P(t)=_0^1\rho (p,t)f(p)𝑑p.$$
(2)
We assume that the distribution of $`p`$ among the users is fairly uniform, i.e, $`f(p)=1`$, which gives,

$$P(t)=\frac{1}{t(t+1)}$$
(3)
which in asymptotic limit gives $`P(t)t^2`$. Fig 1 shows a comparison of the real data and Eqn.(3).

We have studied the empirical distribution of login times in e-mail servers. This observed behavior is unique for e-mail servers and differ from general purpose computers and computational servers, which show a peaked distribution with an exponential (or stretched exponential) decay. The distribution of login times in e-mail servers, however, also has a peak which is not observed due to the fact that the minimum time resolution there is greater than the typical time-scale at which this peak occurs. The origin of the observed distribution is but a result of human behavior. Human behavior giving rise to frustration and fluctuations have been approximated by a simple mean-field theory. Our mean-field result is compatible with the real data and can reasonably explain the origin of scale-free nature of login time distribution. Apart from the e-mail systems studied, this interesting behavior of a large number of components may also lead to similar observations in systems with similar dynamical topology as in other modes of social communication like telephone network and verbal conversation.

###### Acknowledgements.

The author acknowledges B. K. Chakrabarti, P. Bhattacharyya for discussions and K. Das for technical assistance with the data. The author also acknowledges datasrc  H. Scherrer and A. K. Chattopadhyay for data from MPIPKS Dresden, P. A. Sreeram for data from SNBNCBS, S. Chattopadhyay, A. DeSarkar, B. Muthukumar and I. Sarkar and System Administrators of SINP e-mail servers for data from SINP, P. Ray and G. Subramoniam for data from IMSc, Chennai, V. Shenoy and C. Karthik for data from IISc-MRC & Srinivas V. C. for data from IISc-Physics, Bangalore.


text:
# A search for periodic modulations of the solar neutrino flux in Super-Kamiokande-I

## I Introduction

It is widely accepted that neutrinos have masses based on observation of their oscillation ref:oscillation  which suggests the presence of physics beyond Standard Model. Neutrino spin orientation would be rotated by the magnetic field in the Sun if neutrinos are Dirac particles with non-vanishing magnetic moments ref:sfp . Such neutrino spin-flavor precession will result in a reduction of the observed solar neutrino fluxes because the right-handed neutrinos are sterile and will be unobserved in a detector. On the other hand, if neutrinos are Majorana particles with a flavor-changing transition magnetic moment, the solar magnetic field could flip $`\nu _e`$ into $`\overline{\nu }_\mu `$ or $`\overline{\nu }_\tau `$ and lower the observed solar neutrino fluxes due to the reduced elastic scattering cross section. However, any such possibility must not allow for a flux of solar $`\overline{\nu }_e`$ above current experimental limits ref:sk-nuebar . Based on the above context, any time variation of the solar magnetic fields could introduce periodic modulations in the observed solar neutrino fluxes. There could be semiannual (seasonal) variations of the observed solar neutrino flux because of the changing magnetic field caused by the 7.25 degree inclination of solar axis with respect to the ecliptic plane. A short-time variation might be expected due to the 27-day rotation of the Sun.

Independent of magnetic effects, the 8B neutrino flux is sensitive to the core temperature of the Sun, approximately proportional to $`T_c^{25}`$ ref:jnb-temp . If the solar core temperature changes in time, this might induce a time variation in the observed solar neutrino flux.

There has been several attempts to search for possible periodic variations in measured solar neutrino fluxes. The Homestake experiment appeared to exhibit 11-year periodic modulations on the data which anticorrelated with solar activity ref:homestake . However, Kamiokande and other experiments have not provided any evidence for a time variation of the neutrino flux outside of statistical fluctuations ref:kamiokande . Observation of a periodic modulation in the solar neutrino flux would provide a significant addition to our understanding of solar dynamics and the magnetic properties of neutrinos, and possibly require modification of the Standard Solar Model.

Super-Kamiokande (SK) is a 50 kton water Cherenkov detector located in Kamioka, Japan. Solar neutrino data were collected at SK from May 31st, 1996 to July 15th, 2001 yielding a total detector live time of 1,496 days. This data taking period is known as SK-I. A detailed description of SK can be found elsewhere ref:solar ref:sknim2003 . The solar neutrino signal is extracted from the data using the cos$`\theta _{sun}`$ distribution, the angular deviation between the Sun and the reconstructed direction of events with total energies ranging between 5 and 20 MeV ref:sk1258flux . From the strong forward peak due to elastic scattering of 8B solar neutrinos on electrons, $`22,400\pm 200`$(stat.) solar neutrino interactions were observed in 22.5 ktons of fiducial volume. The relatively high yield of real-time events in SK, 15 events per day, allows a search for short-time periodic modulations in the observed neutrino fluxes.

## II SK 10-day long sampled solar neutrino data

The solar neutrino data, acquired over 1,871 elapsed days from the beginning of data-taking, are divided into roughly 10-day long samples as listed in Table 1. The time period of each 10-day sample is chosen from consecutive 10-day periods. Note that the start time of a 10-day sample is determined by the beginning of data-taking on the first day of the sample, and the stop time by the finish of data-taking on the last day. There are on and off periods of data-taking in the 10-day interval and thus the timing of each sample is calculated as a mean of the start and end times and corrected by SK livetime. Hence the mean time is not necessarily an exact division of the time interval for each 10-day sample. Figure 1 shows the measured solar neutrino fluxes of the 10-day samples. All given uncertainties are statistical and estimated by asymmetric Gaussian approximation of the unbinned maximum likelihood fit to the cos$`\theta _{sun}`$ distributions. The measured solar neutrino fluxes are corrected for the 1/$`R^2`$ (squared average distance in units of A.U.) variation caused by the eccentricity of the Earth’s orbit around the Sun.

## III Search for Periodicity

The Lomb periodogram method ref:nr13.8 , a spectral analysis for unevenly sampled data, is applied to search for possible periodicities in the measured 10-day long fluxes. The method finds periodicities based on maximum deviation of data relative to a constant in time. The normalized Lomb power is given by

$$P_N(w)\frac{1}{2\sigma ^2}\left\{\frac{[\underset{j}{}(\varphi _j\overline{\varphi })\mathrm{cos}\omega (t_j\tau )]^2}{_j\mathrm{cos}^2\omega (t_j\tau )}+\frac{[\underset{j}{}(\varphi _j\overline{\varphi })\mathrm{sin}\omega (t_j\tau )]^2}{_j\mathrm{sin}^2\omega (t_j\tau )}\right\}$$
(1)
where $`\varphi _j`$ is the measured flux in the $`j`$-th bin, $`t_j`$ is the mean time with livetime correction in the $`j`$-th bin, $`\omega `$ is the frequency being tested, $`\overline{\varphi }_{i=1}^N\varphi _i/N`$, $`\sigma ^2_{i=1}^N(\varphi _i\overline{\varphi })^2/(N1)`$, and offset time $`\tau `$ is defined by $`\mathrm{tan}(2\omega \tau )=_j\mathrm{sin}2\omega t_j/_j\mathrm{cos}2\omega t_j`$. Under the null hypothesis, it is expected that the data are Gaussian random values and the Lomb normalized power ($`P_N(w)`$) is distributed exponentially with unit mean. Here we define the Confidence Level (C.L.) of the given Lomb Power $`P_N`$ as $`(1e^{P_N})^M\times 100`$%, the probability of no frequencies being scanned giving values larger than $`P_N`$ due to random fluctuation. The number of independent frequencies scanned, M, is given by approximately twice of the number of data points, i.e. M = 2 $`\times `$ 184 = 368. The frequency ranges from 0.00020 day-1 to 0.09870 day-1 and 18,363 frequencies are scanned in this range. The resulting periodogram is shown in Figure 2 together with its confidence level diagram. The maximum power appears at frequency f=0.0726 day-1 (or time period T=13.76 days) with Lomb power 7.51 corresponding to 81.70% C.L.

Figure 3 shows a normal distribution of the powers in the Lomb periodogram for SK solar neutrino data. Here the powers are calculated for 368 independent frequencies from 0 day-1 to 0.1 day-1.

As a consistency check for the confidence level, 10,000 MC experiments are generated based on the observed timing information and the measured solar neutrino fluxes of the 10-day long data samples. The measured solar neutrino flux values are simulated according to a random Gaussian fluctuation. For making these null modulation samples the average measured flux ($`2.33\times 10^6`$ cm-2 s-1) is taken as a Gaussian mean, and the standard deviation of the the measured flux ($`0.32\times 10^6`$ cm-2 s-1) is taken as a Gaussian error ref:horne\_ApJ1986 . The Lomb method is applied to each MC experiment to obtain a periodogram. The Lomb power distribution has an exponential distribution as expected for the null modulation. The maximum power of each periodogram is selected. Figure 4 shows a distribution of maximum powers for the MC experiment sets. Out of 10,000 simulated experiments, 19.58% have maximum powers larger than 7.51. This demonstrates that the confidence level for the T=13.76 day period of SK data is consistent with that of no modulation.

In order to extend the investigation to shorter time modulation more precisely, we repeat the search for periodicities in a time interval of 5 days (see Figure 5 and Table 2) using the Lomb method. Figure 6 shows a periodogram for the 5-day long samples. For the 5-day long sample, the frequency ranges from 0.00020 day-1 to 0.19187 day-1 and 35,763 frequencies are scanned in this range. The number of independent frequencies scanned is M=716. The maximum Lomb power of 7.35 (63.09% C.L.), consistent with no modulation, is found at frequency f=0.1197 day-1 (T=8.35 days). In this case, the Lomb power for the period of 13.76 days (f=0.0726 day-1) is 1.35 and is no longer the maximum. This provides additional confirmation that the 13.76 day period in 10-day long sample is a statistical artifact.

Our results are in clear disagreement with the previous analyses of SK data by Milsztajn ref:milsztajn2003  and Caldwell and Sturrock ref:sturrock\_ApJ2003 , in which they claim to find a significant periodic variation of our data with 13.76 days’ periodicity. The main difference between their analyses and ours is that they have used uniform 10-day binned data in which the time of each data point is exactly centered between the start and stop times of each 10 day period ref:milsztajn2003  or integrating expected flux from the start time to the end time of each 10-day bin without taking into account breaks of runs ref:sturrock\_ApJ2003 . On the other hand, we use the correct time for the data contained in each 10 day bin based on the livetime weighted mean time as described above. We investigate 5-day binned data as well to find no significant periodic variation.

## IV Sensitivity of Finding a Periodicity vs. Modulation Amplitude

We have studied the sensitivity of the SK-I solar neutrino data to find if a true periodicity indeed exists. The sensitivity will depend on the amplitude of modulation because of experimental uncertainties in the SK-I data. One thousand MC experiment sets are generated to simulate the SK 10-day long solar neutrino fluxes for a given period and modulation amplitude. While introducing modulations of periods from 5 days to 500 days and varying the magnitude of modulation from 1% to 90% of the averaged total flux ($`2.33\times 10^6`$ cm-2 s-1), we repeatedly generate 1,000 MC experiments at each amplitude. An individual MC experiment set is generated in the following way. The modulated solar neutrino flux is generated for entire SK-1 data-taking period, and the mean flux is calculated for the individual 10-day long time periods. The MC neutrino flux of every 10-day long sample is then generated by a random Gaussian fluctuation from the mean flux value. The Gaussian error is taken as the error in the measured neutrino flux of each 10-day sample as listed in Table 1. The Lomb periodogram method is applied to every MC experiment. A reconstructed period is determined by the primary modulation having the maximum Lomb power larger than 9.82 (98% C.L.). We calculate the probability of correctly finding the generated period at each modulation amplitude. The observed period is assigned to a correct one if its corresponding frequency is consistent with 1/(input period)$`\pm `$ 0.1/368/2 for the 10-day long sample, to consider the error on the reconstructed period. The probability of reconstructing the input periodicity becomes significantly reduced as the modulation amplitude becomes less than the magnitude of the measured neutrino flux error. The MC study is repeated for 5-day long samples. In the 5-day long sample case, the selection criterion of correct frequency is 1/(input period)$`\pm `$ 0.2/716/2. Figure 7 shows 95% sensitive contours of modulation amplitude vs. period for 10-day and 5-day long samples. The sensitivity is defined as a probability of finding the correct periodicity as a primary modulation.

According to Figure 7, the sensitivity of finding the correct period varies rapidly near the sampling time periods of 5 $``$ 20 days. In general, shorter modulation periods have worse sensitivity for a given sampling frequency, while the 5-day long sampled data have better sensitivity than those of the 10-day data, especially for shorter modulation periods. Based on the MC study, the Lomb method is expected to find the true periodicity of the SK-I data, if any, with an efficiency of 95% at a modulation amplitude of 10% of the averaged flux and a period of longer than 40 days (20 days) in case of the 10-day (5-day) long sample.

## V Summary

We have presented the measured solar neutrino fluxes of 10-day long samples and of 5-day long samples using all 1,496 days of SK-I data. No significant periodicity was found in the SK-I solar neutrino data when a search was made to look for periodic modulations of the observed fluxes using the Lomb method. Based on a MC study, we have obtained the probability of finding a true periodicity in the SK-I data as a function of the modulation magnitude. The Lomb method should have found a periodic modulation in the SK-I solar neutrino data of 10-day (5-day) long samples if the modulation period were longer than 40 days (20 days) and its magnitude was larger than 10% of the average measured neutrino fluxes. Based on the observation of no significant periodicity, SK-I data exclude modulations greater than 10% of the 8B neutrino flux arising as a result of more than 0.4% changes in the solar core temperature, allowing a new measure of the solar core’s stability.

###### Acknowledgements.

The authors acknowledge the cooperation of the Kamioka Mining and Smelting Company. The Super-Kamiokande has been built and operated from funding by the Japanese Ministry of Education, Culture, Sports, Science and Technology, the U.S. Department of Energy, and the U.S. National Science Foundation. This work was partially supported by the Korean Research Foundation (BK21) and the Korea Ministry of Science and Technology.


text:
# The Oklo Constraints on Alpha-Decay Half-Lives

## 1 Introduction

Beginning with Shlyakhter  in 1976, Oklo natural reactor data have been used by many authors \[2-6\] to put limits on the time variation of ”constants” such as the fine structure constant and the strong coupling constant. As noted by Damour and Dyson , such constraints generally refer to the value of the ”constants” at the time of the reactor shut down, since most prior evidence is diminished by the reaction process. Also, it cannot be ruled out that two or more factors may change concurrently, and one may cancel the effect of the other on whatever we are able to calculate for comparison with experiment.

In this work, we will use the constraints provided by the samarium-149 cross section to provide a self-consistent limit on the value of the uranium-238 decay constant at the time of reactor shutdown. This implies a limit on possible differences between the strong coupling constant at the time of reactor shutdown and today. Pierronne and Marquez  developed the complex eigenvalue method for calculating alpha-decay half-lives by matching the wavefunction of the alpha particle, and its derivative, across a boundary between a square well representing the nuclear interior and a coulomb potential representing the exterior of the nucleus outside the range of nuclear forces. Chaffin, Banks, Gothard, and Tuttle  have further developed this method into Mathematica notebooks permitting easy calculation of relevant variations.

Experimental data on alpha particle scattering \[11-14\] have long been known to validate a cluster model in which an alpha particle moves in the average field provided by the rest of the nucleus. Four-nucleon transfer reactions were observed on target nuclei that are the daughters of an alpha emitter. The comparison of reaction rates with the corresponding alpha-decay rates indicates the tendency for the four transferred particles to form a true alpha particle. The data give evidence that nodes in the alpha particle wavefunction do in fact exist. When the nuclear potential seen by the alpha particle is represented by a simple square well, the effect of changing the depth of the well is to change the alpha particle wavefunction. Usually this change is slight, but at certain critical values the number of nodes in the wavefunction can change precipitously, with a corresponding change of more than an order of magnitude in the decay constant. Hence, we find that, although the Oklo data provide Shlyakhter’s very stringent limits on the variation in, for instance, the Sm-149 absorption cross-section for neutrons, this may not translate into such stringent limits on nuclear alpha-decay half-lives.

In section 2 of this paper we will present the method of calculation of the Samarium cross-section at the time of the Oklo reactions, and indicate how this may be used to find the depth of the square well representing the nuclear potential. Then in section 3 we review the Pierronne and Marquez  method for calculating the decay constant from the depth of the well and other parameters. We then present the results of a calculation for the case of alpha-decay of U-238.

## 2 Using Oklo Data to find the Samarium-149 Cross-Section

We have used actual data of element concentrations from Oklo samples to find the implied Sm-149 absorption cross-section for neutrons. Procedures very similar to those of Damour and Dyson  and Fujii et al..  were used. As Fujii et al. explained in their Appendix B, coupled linear differential equations may be written down for the time derivatives of the Nd-143, Sm-147, U-235, and U-238 concentrations. These equations may be solved for the restitution factor $`C_{rz}`$, epi-thermal index $`r_{rz}`$ , and the neutron fluence appropriate for each sample. For each sample the ratio k = $`N_{235}^0`$/$`N^{nat}`$, the ratio of the initial U-235 abundance to the natural abundance of samarium, is also found using Fujii et al. ’s procedure.

As explained in the Damour and Dyson  and/or Fujii et al.  papers, coupled linear differential equations may also be written for the Sm-147, Sm-148, and Sm-149 concentrations. To obtain a working equation we use equation (16) of Fujii et al. , corrected for a typographical error in the subscript for the Sm-149 absorption cross-section, yielding the result for the abundance ratios at the end of the reactor activity. If the time for the end of reactor activity is $`t_1`$, and $`N_{147}`$ , $`N_{148}`$, and $`N_{149}`$ are the concentrations of Sm-147, Sm-148, and Sm-149, respectively, the appropriate equation is:

$$\begin{array}{c}\frac{N_{147}\left(t_1\right)+N_{148}\left(t_1\right)}{N_{149}\left(t_1\right)}=\hfill \\ \frac{\frac{\widehat{\sigma }_{f235}kY_{147}}{\widehat{\sigma }_a}\left[1\mathrm{exp}\left(\widehat{\sigma }_a\widehat{\varphi }t_1\right)\right]+\left(R_{147}^{nat}+R_{148}^{nat}\right)}{\frac{k\widehat{\sigma }_{f235}Y_{149}}{\widehat{\sigma }_a\widehat{\sigma }_{149}}\left[\mathrm{exp}\left(\widehat{\sigma }_{149}\widehat{\varphi }t_1\right)\mathrm{exp}\left(\widehat{\sigma }_a\widehat{\varphi }t_1\right)\right]+R_{149}^{nat}\mathrm{exp}\left(\widehat{\sigma }_{149}\widehat{\varphi }t_1\right)}\hfill \end{array}$$
(1)
Here the R’s are the relative fractional natural abundances of the isotopes of samarium indicated by the subscripts, $`Y_{147}`$ and $`Y_{149}`$ are the fission yields of Sm-147 and Sm-149, respectively. Using this equation, we used an iteration procedure called FindRoot in Mathematica to obtain the samarium-149 absorption cross-section.

Results are given in Table I, where they are compared to Fujii et al. ’s results.

The table gives the cross-sections obtained for samples inside the core of reactor zone 10, and average cross-section of $`90.29\pm 7.41`$ kilobarns. The errors are the one standard deviation statistical errors. The temperature of the Oklo ores at the time of the reactions is unknown but thought to have been around 300 to 400 degrees Celsius. This temperature would, of course have varied over the lifetime of the nuclear reactor, and is an important systematic source of error in determining, from the data, the Sm-149 cross-section at the time of the reactions.

## 3 Finding the U-238 Decay Constant

In a simple nuclear model in which a square well is used for the potential seen by the alpha particle inside the nucleus, the Sm-149 cross-section is related to the depth of the potential well. The cross-section as a function of energy is given by the Breit-Wigner shape:

$$\sigma _r=A\frac{\mathrm{\Gamma }_n\mathrm{\Gamma }_\gamma }{\left(EE_r\right)^2+\mathrm{\Gamma }_{tot}^2/4}$$
(2)
In this expression the widths $`\mathrm{\Gamma }`$ do not depend directly on the depth $`V_0`$ of the potential well, but the resonant energy $`E_r`$ is a constant minus $`V_0`$, as was shown, for example, by Weber, Hammer, and Zidell . The spread in values of the samarium cross-section found in section 2, thus can be related to a spread in values of $`V_0`$. For U-238, the half-life of $`4.47\times 10^9`$ years is reproduced by a $`V_0`$ of 110.5 Mev. From the data of section 2, the integrated Sm-149 cross section is found to be $`90.3\pm 7.4`$ kilobarns (see Table I). The samarium resonant energy is $`E_r`$ =0.0958325 eV and the corresponding one standard deviation spread is $`\pm 0.0102`$ eV. While this spread might seem small, the small change in $`V_0`$ can lead to a large change in the alpha-decay half life.

To test the variability of the decay constant, Mathematica notebooks were written. For the square well potential on the inside of the nucleus, the Mathematica notebook gave answers essentially equivalent to those of Chaffin, Gothard and Tuttle , which used a harmonic oscillator potential for the interior region, where the nuclear potential is felt by the alpha particle. In the course of this work, it was discovered that, as the nuclear potential well depth is changed, and the nuclear radius changes slightly, it is possible to have a sudden change in the number of nodes of the real part of the alpha particle’s wavefunction. This was modeled for both the harmonic oscillator and square well potentials, with nearly the same results for either notebook.

The change in the number of nodes causes the probability of tunneling to change by about a factor of 13.5, as shown by the discontinuity in the graph of Figure 2. The graph shows a large, precipitous rise in the value of the decay constant $`\lambda `$, at certain critical values of $`V_0`$, similar to the postage-stamp function.

## 4 Conclusions

The calculations of the decay constant versus well depth $`V_0`$ show that even a small changes in the strong coupling constant could result in a U-238 decay constant change by a factor of 13.5, more than one order of magnitude. The effect is tied to the change in the number of nodes in the alpha-particle wavefunction when the well depth passes critical values.
References

 Shlyakhter, A. I., Nature 264, 340 (1976).

 Damour, T. and F. J. Dyson, Nuclear Physics B480, 37-54 (1996).

 Fujii, Y., A. Iwamoto, T. Fukahori, T. Ohnuki, M. Nakagawa, H. Hidaka, Y. Oura, and P. Möller, Nuclear Physics B 573, 377-401 (2000).

 Irvine, J. M., Contemporary Physics 24(5): 427-437 (1983).

 Irvine, J. M., Philosophical Transaction of the Royal Society A310, 239-243 (1983).

 Irvine, J.M. and R. Humphreys, Progress in Particle and Nuclear Physics 17, 59-86 (1986).

 Marquez, L. Phys. Rev. C25(2), 1065-1067. (1982).

 Pierronne, M., and L. Marquez, Zeitschrift für Physik A286, 19-25 (1978).

 Chaffin, E. F., and D. S. Banks, nucl-th/0206020.

 Chaffin, E. F., N. W. Gothard, and J. P. Tuttle, nucl-th/0105070.

 Anyas-Weiss, N., J.C. Cornell, P.S. Fisher, P. N. Hudson, A. Menchaca-Rocha, D.J. Millener, A.D. Panagiotou, D.K. Scott, and D. Strottman, Physics Reports C12(3), 201-272. (1974).

 Buck, B. Nuclear-structure information from several-nucleon transfer reactions, in Nuclear spectroscopy and nuclear reactions with heavy ions, International School of Physics Enrico Fermi course 62. (North Holland, Amsterdam, 1976).

 Davies, W. G., R.M. DeVries, G.C. Ball, J.S. Forster, W. McLatchie, D. Shapira, D., J. Toke, and R.E. Warner, Nuclear Physics A 269, 477-492 (1976).

 Wildermuth, K., and W. McClure, Cluster representations of nuclei  (Springer tracts in modern physics No. 41). (New York: Springer-Verlag, New York, 1966).

 Weber, T.A., C.L. Hammer, and V.S. Zidell, American Journal of Physics 50(9): 839-845 (1982).

 Green, A. E. S. and K. Lee, Phys. Rev. 99, 772-777 (1955).

Figure Captions

Figure 1. Sudden change in the number of nodes. The harmonic oscillator wave function for well depths of 58 MeV (a) and 54 MeV (b). The x-axis is the radial coordinate of the alpha particle, $`T=\rho /(2\eta )`$, where $`\rho `$ and $`\eta `$ are defined in Green and Lee . Figure 1a shows the harmonic oscillator wave function for a well depth of 58 MeV. Figure 1b shows what happens when the well depth is changed to 54 MeV, without changing the alpha particle energy. If one counts the number of nodes in Figure 1a, there are nine, not counting the ones at zero and infinity. For Figure 1b, there are eight nodes, a reduction by one.

Figure 2: The Decay Constant $`\lambda `$ of U-238 versus Well Depth.