format: text

text:
# Neutrino Physics: an Update

## 1 Introduction

Several years ago, we authored a paper in this journal entitled “Neutrino Physics,” hereafter called I, in order to encourage inclusion of material involving neutrinos into the introductory curriculum. We noted at the time that neutrinos, with new experiments about to provide first data, might continue to be a popular press item. This prediction has proven true:

* Ray Davis, the founder of the field of experimental solar neutrino physics, shared the 2002 Nobel Prize in physics with Mastoshi Koshiba, who led the Kamioka solar neutrino experiment, and Riccardo Giacconi .
* Results from the Sudbury Neutrino Observatory (SNO) resolved the solar neutrino puzzle, showing that approximately two-thirds of these neutrinos oscillate into other flavors before reaching earth .
* The KamLAND experiment, in which antineutrinos from Japanese power reactors were detected, confirmed the SNO results and further narrow the allowed range of neutrino mass differences .
* The Wilkinson Microwave Anisotropy Probe (WMAP) measured subtle temperature differences within the oldest light in the universe, from the epoch when atoms first formed 380,000 years after the Big Bang . When combined with the results of large scale structure studies , a new bound on the sums of the neutrino masses is obtained.

In addition there has been published (and disputed) evidence for the existence of neutrinoless double beta decay which, if confirmed, would show that one of the standard model’s most important symmetries, the conservation of lepton number, is violated . Thus we decided to bring our early paper up to date by explaining here the importance and implications of the new results.

In Section 2 we present a much abbreviated summary of the material presented in I. In Section 3 we discuss the SNO results and how they resolved the puzzling discrepancies Ray Davis first uncovered 30 years ago. In Section 4 we discuss KamLAND, the first terrestrial experimental to achieve the sensitivity to the small neutrino mass differences relevant to solar neutrino experiments. In Section 5 we describe marvelous new cosmological probes of large scale structure and of the time when atoms were first formed, and why the new data may soon challenge recent double beta decay claims. We summarize where we stand in neutrino physics – including the new discoveries that may be soon be within reach – in Section 6.

## 2 The Two-Minute Review

In I we summarized the basics of neutrino physics, including neutrino history, properties, and implications for contemporary physics. For the purposes of the present work, we note that in the so-called standard model of particle/nuclear physics, which is consistent with nearly all present experimental information, there exist three massless neutrino types—$`\nu _e,\nu _\mu ,\nu _\tau `$—which are produced with purely left-handed helicity in weak interaction processes. If the neutrino were shown to have a nonvanishing mass, it would be the first clear failure of this 30-year-old standard model and the first proof of the existence of the “particle dark matter” that appears necessary to explain the structure and expansion of our universe. Nonzero neutrino masses were suggested by the results of experiments measuring the flux of solar neutrinos, which are produced as a byproduct of the thermonuclear reactions occurring in the high-temperature core of our sun . Additional strong evidence comes from the study of atmospheric neutrinos, which are produced when high energy cosmic rays collide with the upper atmosphere, producing pions and other particles that then decay into neutrinos . The largest of the various atmospheric neutrino experiments is SuperKamiokande, a detector in a mine in the Japanese alps that contains 50,000 tons of ultra-pure water. SuperKamiokande’s precise data show that the flux of muon-type atmospheric neutrinos arriving from the opposite side of the earth, which have travelled a long distance to reach the detector, is depleted.

Both the solar and atmospheric neutrino results can be explained quantitatively if neutrinos are massive and if the mass eigenstates are not coincident with the weak interaction eigenstates $`\nu _e,\nu _\mu ,\nu _\tau `$, i.e., if the neutrinos produced in weak interactions are combinations of the various mass eigenstates. This is exactly what is known to occur in the analogous case of the quarks. The relation between the neutrino mass and weak-interaction eigenstates is described by a unitary matrix: for three neutrinos, the mass and weak-interaction eigenstates can be viewed as two distinct 3D coordinate systems. The unitary matrix specifies the three rotations describing the orientation of one coordinate system relative to the other.

The reduced flux of atmospheric muon neutrinos (see Figure 1) is quantitatively explained by $`\nu _\mu \nu _\tau `$ oscillations governed by maximal mixing, $`\theta _{23}\pi /4`$: the relationship between the mass eigenstate $`\nu _2`$ and $`\nu _3`$ and the flavor eigenstates $`\nu _\mu `$ and $`\nu _\tau `$ involves a rotation of 45 degrees. The magnitude of the mass difference $`\delta m_{23}=m_3^2m_2^2`$ between the two equal components making up the flavor eigenstates is $`2\times 10^3`$ eV2. When I was written, the favored solution to the solar neutrino problem was also neutrino mixing, but was described by a much smaller mixing angle $`\theta _{12}`$ specifying the relationship between mass eigenstates $`\nu _1`$ and $`\nu _2`$ and the flavor eigenstates $`\nu _e`$ and $`\nu _\mu `$. The effects of this small mixing are magnified by matter effects in the sun. This is the so-called MSW effect and is described in I. But the results of the Sudbury Neutrino Observatory provided a bit of a surprise.

## 3 The SNO Experiment

Because both charged and neutral currents contribute to the reaction important to SuperKamiokande,

$$\nu _x+e^{}\nu _x+e^{},$$
(1)
the experimentalists cannot easily distinguish $`\nu _e`$s from the $`\nu _\mu `$s and $`\nu _\tau `$s: the detector records both fluxes, though with a reduced sensitivity (0.15) for the heavy-flavor types. The reaction produces energetic recoil electrons which generate Cerenkov radiation that is recorded in an array of phototubes surrounding the detector. As the cross section is sharply forward peaked, the correlation with the position of the sun can be used to “cut” background contributions associated with cosmic rays and radioactivity in the rock walls surrounding the detector. Because the threshold for electron detection is $`6`$ MeV, only the high energy portion of the 8B solar neutrino flux is measured. These are the same neutrinos that dominate the radiochemical measurements of Ray Davis: Superkamiokande confirmed that this flux was substantially below that predicted by the standard solar model (SSM)

$$\varphi _{SSM}(\nu _x)=5.44\times 10^6\mathrm{cm}^2\mathrm{sec}^1$$
(2)
When the Davis and SuperKamiokande results were combined with those from the gallium experiments SAGE and GALLEX, the resulting three constraints on the three principal solar neutrino sources (8B, 7Be, and the low-energy pp fluxes) produced a surprising result. No combination of these fluxes could reproduce the combined data well. Though circumstantial, this evidence indicated that the solution to the solar neutrino problem would not be found in the SSM, but instead would require new particle physics.

The favored explanation became neutrino oscillations which, as we have just summarized, occur for massive neutrinos when the weak and mass eigenstates do not coincide. The SuperKamiokande discrepancy, a solar neutrino rate less than half that expected from the SSM, would require that approximately two-thirds of the high energy electron neutrinos generated in the solar interior oscillate into other species – $`\nu _\mu ,\nu _\tau `$ – before reaching earth. Low-energy $`\nu _\mu `$s and $`\nu _\tau `$s are invisible to the Davis and SAGE/GALLEX detectors and scatter electrons in the SuperKamiokande detector with a reduced cross section, as we noted previously.

The key idea behind the Sudbury Neutrino Observatory (SNO) was construction of a detector that would have multiple detection channels, recording the $`\nu _e`$s by one reaction and the total flux of all neutrinos ($`\nu _e`$s + $`\nu _\mu `$s + $`\nu _\tau `$s) by another. This was accomplished by replacing the ordinary water in a water Cerenkov detector with heavy water – D2O instead of H2O. The charged-current (CC) channel that records the $`\nu _e`$s is analogous to the reaction used in the Davis detector

$$\nu _e+dp+p+e^{}.$$
(3)
As the electron produced in this reaction carries off most of the neutrino energy, its detection in the SNO detector (by the Cerenkov light it generates) allows the experimentalists to determine the spectrum of solar $`\nu _e`$s, not just the flux. A second reaction, the neutral current (NC) breakup of deuterium, gives the total flux, independent of flavor (the $`\nu _e`$, $`\nu _\mu `$, and $`\nu _\tau `$ cross sections are identical),

$$\nu _x+dn+p+\nu _x.$$
(4)
The only signal for this reaction in a water Cerenkov detector is the neutron, which can be observed as it captures via the $`(n,\gamma )`$ reaction. SNO is currently operating with salt added to the water, as Cl in the salt is an excellent $`(n,\gamma )`$ target, producing about 8 MeV in $`\gamma `$s.

While the strategy may sound straightforward, producing such a detector was an enormous undertaking. The needed heavy water – worth about $300M – was available through the Canadian government because of its CANDU reactor program. The single-neutron detection required for the neutral current reaction is possible only if backgrounds are extremely low. For this reason the detector had to be placed very deep underground, beneath approximately two kilometers of rock, so that cosmic-ray muon backgrounds would be reduced to less than 1% of that found in the SuperKamiokande detector. The experimentalists found the needed site in an active nickel mine, the Sudbury mine in Ontario, Canada, where they worked with the miners to carve out a 10-story-high cavity on the mine’s 6800-ft level. Trace quantities of radioactivity were another background concern: if a thimblefull of dust were introduced into the massive cavity during construction, the resulting neutrons from U and Th could cause the experiment to fail. Thus, despite the mining activities that continued around them, the experimentalists constructed their detector to the strictest cleanroom standards. The detector also provided a third detection channel, neutrino elastic scattering (ES) off electrons (Eq. (1)), which we have noted is sensitive to $`\nu _e`$s and, with reduced sensitivity, $`\nu _\mu `$s and $`\nu _\tau `$s.

The ES reaction, of course, provides SNO a direct cross check against SuperKamiokande. SNO’s threshold for measuring these electrons is about 5 MeV. Assuming no oscillations, SNO’s detection rate is equivalent to a $`\nu _e`$ flux of

$$\varphi _{SNO}^{ES}=2.39\pm 0.34(\mathrm{stat})\pm 0.15(\mathrm{syst})\times 10^6\mathrm{cm}^2\mathrm{sec}^1,$$
(5)
a result in excellent accord with that from SuperKamiokande,

$$\varphi _{SK}^{ES}=2.32\pm 0.03(\mathrm{stat})\pm 0.06(\mathrm{syst})\times 10^6\mathrm{cm}^2\mathrm{sec}^1.$$
(6)
The greater accuracy of the SuperKamiokande result reflects the larger mass (50 kilotons) and longer running time of the Japanese experiment. (SNO contains, in addition to the one kiloton of heavy water in its central acrylic vessel, an additional seven kilotons of ordinary water that surrounds the central vessel, helping to shield it.)

The crucial new information provided by SNO comes from the two reactions on deuterium. The CC current channel is only sensitive to $`\nu _e`$s. Under the assumption of an undistorted 8B neutrino flux, SNO experimentalists deduced

$$\varphi _{SNO}^{CC}(\nu _e)=1.75\pm 0.07(\mathrm{stat})\pm 0.12(\mathrm{sys})\pm 0.05(\mathrm{theory})\times 10^6\mathrm{cm}^2\mathrm{sec}^1.$$
(7)
The CC flux is less than that deduced from the ES rate, indicating that $`\nu _\mu `$s and $`\nu _\tau `$s must be contributing to the later. From the difference between the SuperKamiokande ES and the SNO CC results

$$\delta \varphi =0.57\pm 0.17\times 10^6\mathrm{cm}^2\mathrm{sec}^1$$
(8)
and recalling that the $`\nu _\mu /\nu _\tau `$ ES cross section is only 0.15 that for the $`\nu _e`$, one deduces the heavy-flavor contribution to the solar neutrino flux

$$\varphi (\nu _\mu /\nu _\tau )=3.69\pm 1.13\times 10^6\mathrm{cm}^2\mathrm{sec}^1.$$
(9)
That is, approximately two-thirds of the solar neutrino flux is in these flavors.

While the first SNO analysis was done in the manner described above, a second publication gave the long awaited NC results. This allowed a direct and very accurate determination of the flavor content of solar neutrinos, without the need for combining results from two experiments. The published NC results were obtained without the addition of salt to the detector: the neutron was identified by the 6.25 MeV $`\gamma `$ ray it produces by capturing on deuterium. The resulting total flux, independent of flavor, is

$$\varphi _{SNO}^{NC}(\nu _x)=5.09\pm 0.44(\mathrm{stat})\pm 0.45(\mathrm{syst})\times 10^6\mathrm{cm}^2\mathrm{sec}^1.$$
(10)
Combining with the CC signal yields

$`\varphi _{SNO}(\nu _e)=1.76\pm 0.05(\mathrm{stat})\pm 0.09(\mathrm{syst})\times 10^6\mathrm{cm}^2\mathrm{sec}^1`$
$`\varphi _{SNO}(\nu _\mu /\nu _\tau )=3.41\pm 0.45(\mathrm{stat})\pm 0.46(\mathrm{syst})\times 10^6\mathrm{cm}^2\mathrm{sec}^1`$ (11)
The presence of heavy-flavor solar neutrinos and thus neutrino oscillations is confirmed at the 5.3$`\sigma `$ level! Furthermore the total flux is in excellent agreement with the predictions of the SSM—Eq. 2, an important vindication of stellar evolution theory.

The SNO analysis is summarized in Figure 2, which shows the three bands corresponding to the CC, NC, and ES measurements coinciding in a single region. These results can now be combined with other solar neutrino measurements to determine the parameters – the mixing angle and mass-squared difference – governing the oscillations. At the time I was written, there were several contending solutions, though the data favored one characterized by a small mixing angle (thus called the SMA solution). Figure 3 shows that the SNO result has determined an oscillation solution that, at 99% confidence level, is unique – and as in the atmospheric neutrino case, it has a large mixing angle, $`\theta _{12}30`$ degrees. This LMA oscillation is clearly distinct from that seen with atmospheric neutrinos, with $`\delta m_{12}^2=m_2^2m_1^2`$ centered on a region $`8\times 10^5`$ eV2.

The discovery that the atmospheric and solar neutrino problems are both due to neutrino oscillations has provided the first evidence for physics beyond the standard model. That neutrinos provided this evidence is perhaps not unexpected: if the standard model is viewed as an effective theory, one largely valid in our low-energy world but missing physics relevant to very high energies, beyond the reach of current accelerators, then a neutrino mass term is the lowest-order correction that can be added to that theory. But a surprise is the large mixing angles characterizing the neutrino oscillations – which contradicts the simple prejudice that neutrino mixing angles might be similar to the small angles familiar from quark mixing. Perhaps this simply reinforces something that should have been apparent at the outset: with their small masses and distinctive mixings, neutrinos likely have an underlying mechanism for mass generation that differs from that of the other standard model fermions.

## 4 The KamLAND Experiment

One remarkable aspect of the solar and atmospheric neutrino discoveries is that the derived oscillation parameters are within the reach of terrestrial experiments. This did not have to be the case – solar neutrinos are sensitive to neutrino mass differences as small as 10-12 eV2, for which terrestrial experiments would be unthinkable.

The first terrestrial experiment to probe solar neutrino oscillation parameters, KamLAND, very recently reported first results. The acronym KamLAND stands for Kamioka Liquid scintillator Anti-Neutrino Detector. The inner detector consists of one kiloton of liquid scintillator contained in a spherical balloon, 13m in diameter. The balloon is suspended in the old Kamioka cavity (where SuperKamiokande’s predecessor was housed) by Kevlar ropes, with the region between the balloon and an 18m-diameter stainless steel containment vessel filled with additional scintillator (serving to shield the target from external radiation). Several Japanese power reactors are about 180 km from the Kamioka site, and the electron antineutrinos emitted by nuclear reactions in the cores of these reactors can be detected in KamLAND via the inverse beta decay reaction

$$\overline{\nu }_e+pe^++n,$$
(12)
where the $`e^+`$ is seen in coincidence with the delayed 2.2 MeV $`\gamma `$ ray produced by the capture of the accompanying neutron on a proton. This coincidence allows the experimentalists to distinguish $`\overline{\nu }_e`$ reactions from background.

¿From the reactor operations records, which the power companies have made available, KamLAND experimentalists can calculate the resulting flux at Kamiokande to a precision of $`2`$%, in the absence of oscillations. Thus, if a significant fraction of the reactor $`\overline{\nu }_e`$s oscillate into $`\overline{\nu }_\mu `$s or $`\overline{\nu }_\tau `$s before reaching the detector, a low rate of $`e^+`$/capture-$`\gamma `$-ray coincidences will be evident: this is an example of the “disappearance” oscillation technique we described in I. For the 162 ton/yr exposure so far reported by the KamLAND collaboration, the number of events expected in the absence of oscillations is $`86.8\pm 5.6`$. But the number measured is 54 – just 61% of the no-oscillation expectation. From the two-neutrino-flavors oscillation survival probability

$$P(\overline{\nu }_e\overline{\nu }_e)1\mathrm{sin}^22\theta _{12}\mathrm{sin}^2\frac{\delta m_{12}^2L}{4E_\nu }$$
(13)
one obtains the oscillation parameters of Figure 4. KamLAND confirms the LMA solution and significantly narrows SNO’s allowed region (the red area in Figure 4). KamLAND has excellent sensitivity to $`\delta m_{12}^2`$ but less sensitivity to $`\mathrm{sin}^22\theta _{12}`$ (due to uncertainties in the shape of the reactor $`\overline{\nu }_e`$ spectrum). The result is the separation of the SNO LMA allowed region into two parts, with the best-fit $`\delta m_{12}^27\times 10^5`$ eV2, but with a larger mass difference $`1.5\times 10^4`$ eV2 also fitting well. KamLAND is an excellent example of complementary terrestrial and astrophysical measurements: solar neutrino experiments provide our best constraints on $`\theta _{12}`$, but KamLAND places the tightest bounds on $`\delta m_{12}^2`$.

## 5 WMAP, Double Beta Decay, and Neutrino Mass

Despite the wonderful recent discoveries in neutrino physics, there remain quite a number of open questions. Several have to do with the pattern of neutrino masses. All of the experiments described above probe mass differences, not absolute neutrino masses. Furthermore the atmospheric neutrino experiments only constrain the magnitude of $`\delta m_{23}^2`$, and not its sign. As a result, there exist several mass patterns fully consistent with all known data. One choice would assign $`m_3`$ to be the heaviest neutrino, split by the atmospheric mass difference $`\delta m_{23}^22\times 10^3`$ eV2 from a lighter, nearly degenerate pair of neutrinos responsible for solar neutrino oscillations. (This pair is split by solar neutrino mass difference $`\delta m_{12}^210^5`$ eV2.) However, as the sign of $`\delta m_{23}^2`$ is not known, it is also possible that $`m_3`$ is the lightest neutrino, with the nearly degenerate $`m_1`$ and $`m_2`$ heavier. Finally, the best direct laboratory constraint on absolute neutrino masses comes from studies of tritium beta decay, as described in I. Studies of the tritium spectrum near its endpoint energy places a bound of 2.2 eV  on the $`\overline{\nu }_e`$ mass (or more properly, on the principal mass eigenstate contributing to the $`\overline{\nu }_e`$). Consequently, one can add an overall scale of up to 2.2 eV to the mass splittings described above. That is, no terrestrial measurement rules out three nearly degenerate neutrinos, each with a mass $``$ 2.2 eV, but split by requisite $`\delta m_{atmos}^2`$ and $`\delta m_{solar}^2`$.

As discussed in I, the absolute neutrino mass is crucial in cosmology, as a sea of neutrinos produced in the Big Bang pervades all of space. If these neutrinos carry a significant mass, they would constitute an important component of particle dark matter, invisibly affecting the structure and expansion of our universe. Light neutrinos, such as those we have been discussing, decouple from the rest of the matter as relativistic particles, about one second after the Big Bang. Their number density and temperature can be calculated at the time of decoupling and today. Their contribution to the universe’s mass-energy budget is now dominated by their masses,

$$\rho _\nu =0.022\rho _{crit}\underset{i}{}\frac{m_\nu (i)}{\mathrm{eV}},$$
(14)
where $`\rho _{crit}`$ is the critical density that will just close the universe. A variety of cosmological probes suggest that our universe is very close to the critical density—$`\rho /\rho _{crit}=1.0\pm 0.04`$.

¿From the discussion above, we know that at least one neutrino must have a mass of at least $`\sqrt{\delta m_{atmos}^2}`$. Similarly, the tritium beta decay limit places an upper bound on the sum of the masses of 6.6 eV (corresponding to three nearly generate neutrinos of mass 2.2 eV). It follows that the neutrino contribution to dark matter is bounded above and below

$$0.0012\text{ }<\rho _\nu /\rho _{crit}\text{ }<0.15.$$
(15)
This broad range implies that the amount of mass in neutrinos could easily exceed all the familiar baryon matter – stars, dust, gas clouds, and us – visible or invisible: big-bang nucleosynthesis and precision measurements of the cosmic microwave background (CMB) both indicate that $`\rho _{baryons}/\rho _{crit}0.042`$.

However, in the past few years a series of extraordinarily precise measurements have been made in cosmology. One of these is the recent WMAP full-sky map of the CMB and its subtle (few millionths of a degree) temperature anisotropies. This is the oldest light in the universe, the photons that decoupled from matter at the time atoms formed, about 380,000 years after the Big Bang. The CMB temperature anisotropies tell us about the structure of the universe – its clumpiness – at this very early epoch. Measurements of distant SNIa supernovae – a sort of “standard candle” by which astronomers can measure cosmological distances – have constrained the expansion rate and mass/energy budget of the universe. Large-scale surveys, such as the 2dF Galaxy Redshift Survey, have mapped the distribution of visible matter in the universe today and in recent times. The result from combining these and other cosmological probes is a rather sharp constraint on the amount of hot dark matter – in particular, the mass density in neutrinos – that can be allowed, given that our universe has evolved to its present state. One finds

$$\rho _\nu /\rho _{crit}\text{ }<0.022,$$
(16)
that is, an upper bound on the sum of neutrino masses of about 1.0 eV. (Some analyses claim even tighter upper bounds.) This bound is significantly tighter than that of Eq. (14), derived from laboratory data only.

Cosmology now demands that neutrino dark matter can make up no more than 2-3% of the universe’s mass-energy budget and, in particular, is less important than other forms of matter we know about (e.g., nucleons). (We will not go into the disconcerting fact that at least 93% of the universe’s mass-energy budget appears to be dark energy and cold dark matter that we have not yet adequately characterized!) It also tells us that neutrino mass at the $``$ 1 eV level now effects cosmological analyses: such analyses would constrain other cosmological parameters more tightly if neutrino masses were measured, rather than being free parameters that one must dial in cosmological models until unacceptable deviations are found. (The one eV bound was derived in this way.) This underscores how important it is to significantly improve laboratory mass limits.

One possibility is a new-generation tritium experiment: there is a serious effort underway to improve the current bound on the $`\overline{\nu }_e`$ mass to about 0.3 eV, which would then place an upper bound on the sum of the masses of about 0.9 eV . Another possibility – less definitive, perhaps, but with even greater reach – is offered by new-generation neutrinoless double beta decay experiments.

The phenomenon of neutrinoless double beta decay, described in I, tests not only mass, but also whether a standard model symmetry called lepton number conservation is violated. In neutrinoless double beta decay a nucleus spontaneously decays by changing its charge by two units while emitting two electrons,

$$(A,Z)(A,Z+1)+e^{}+\overline{\nu }_e(A,Z+2)+e^{}+e^{},$$
(17)
where the intermediate nuclear state $`(A,Z+1)`$ is virtual. The emitted electrons carry off the entire nuclear energy release, allowing this process to be distinguished from the standard-model-allowed process of two-neutrino double beta decay (where the energy is shared between two electrons and two $`\overline{\nu }_e`$s in the final state). The neutrinoless process clearly violates lepton number, as two leptons (the electrons) are spontaneously produced. (By contrast, as the $`e^{}`$ carries $`l=+1`$ and the $`\overline{\nu }_e`$ $`l=1`$, two-neutrino double beta decay conserves lepton number.)

What conditions will lead to neutrinoless double beta decay? The necessary lepton number violation is present if the neutrino is a Majorana particle—i.e., is identical to its antiparticle. Most theoretical models include Majorana neutrinos: this is part of the mechanism that allows us to understand why neutrinos have masses much smaller than those of the other standard-model fermions, such as electrons and quarks. But the existence of a Majorana neutrino alone is not sufficient because of the exact handedness of neutrinos, which we discussed in I. In the neutrinoless double beta decay reaction written above, the $`\overline{\nu }_e`$ appearing in the intermediate nuclear state was produced in the nucleus by the neutron $`\beta `$ decay reaction $`np+e^{}+\overline{\nu }_e`$. To complete the decay, the antineutrino must be reabsorbed on a second neutron, $`\nu _e+np+e^{}`$. At first glance, if $`\overline{\nu }_e=\nu _e`$, i.e., if the neutrino is its own antiparticle, this reabsorption looks possible. However, this conclusion overlooks the neutrino handedness. In the first step, the $`\overline{\nu }_e`$ produced is righthanded, while the second reaction only proceeds if the $`\nu _e`$ is lefthanded. Thus it would appear that the neutrinoless process is forbidden, even if $`\overline{\nu }_e=\nu _e`$.

This argument, however, overlooks the effects of neutrino mass: a small neutrino mass breaks the exact neutrino handedness, allowing neutrinoless $`\beta \beta `$ decay to proceed, though the amplitude is suppressed by the factor $`m_\nu /E_\nu `$, where $`E_\nu 30`$ MeV is the typical energy of the exchanged neutrino. It follows that neutrinoless double beta decay measures the neutrino mass – at least the Majorana portion of that mass. (Making this statement more precise, unfortunately, takes us beyond the limits of this paper.)

In the simplest case – a single Majorana mass eigenstate dominating the $`\beta \beta `$ decay – the neutrinoless amplitude is proportional to $`U_{ei}^2m_i`$, where $`U_{ei}^2`$ is the mixing probability of the $`i`$th mass eigenstate in the $`\nu _e`$ and $`m_i`$ is the mass. Currently the best neutrinoless $`\beta \beta `$ decay limits are those obtained by the Heidelberg-Moscow and IGEX enriched (86%) 76Ge experiments, which are both probing lifetimes beyond $`10^{25}`$ years – corresponding to roughly one decay per kg-year! These experiments employ Ge crystals – the Ge is both source and detector – containing about 10 kg of active material. Next-generation experiments, using a variety of double beta decay sources (76Ge, 136Xe, 100Mo), have been proposed at the one ton scale. These have as their goals sensitivities to neutrino masses of 10-50 milli-eV, corresponding to lifetimes well in excess of $`10^{27}`$ y. One important motivation for these heroic proposals is the number $`\sqrt{\delta m_{atmos}^2}`$ 55 milli-eV: in several scenarios accommodating the solar and atmospheric oscillation results, this scale plays a role in determining the level at which neutrinoless $`\beta \beta `$ decay might be observed.

A few members of the Heidelberg-Moscow collaboration have claimed that their present results are not a limit, but rather a detection of neutrinoless $`\beta \beta `$ decay, with a best value for the Majorana neutrino mass of $``$ 0.4 eV . This claim has also been strongly criticized by a group that argues that the claimed peak, shown in Figure 5, is not statistically significant. Regardless, this claim will clearly be tested soon, in other $`\beta \beta `$ decay experiments and in future cosmological tests, which promise to soon be probing masses $``$ 0.3 eV.

## 6 Conclusion

In the three years since the publication of I, several very significant neutrino discoveries have been made:

* SNO has shown that approximately two-thirds of the 8B neutrinos that arrive on earth have oscillated into $`\nu _\mu `$s or $`\nu _\tau `$s, thus demonstrating that new neutrino physics is responsible for the solar neutrino puzzle first uncovered by Ray Davis, Jr. Together with the atmospheric neutrino discoveries of SuperKamiokande, this discovery of an effect requiring massive neutrinos and neutrino mixing is the first evidence for physics beyond the standard model. The SNO results for the total solar neutrino flux, independent of flavor, are in excellent agreement with the predictions of the standard model – despite the challenge of calculating a flux that varies as $`T_c^{22}`$(!), where $`T_c`$ is the solar core temperature. The SNO results, when added to other solar neutrino data, isolate a single oscillation scenario, the LMA solution.
* The first terrestrial experiment to probe solar neutrino oscillation parameters, KamLAND, has confirmed the SNO results and further narrowed the LMA range of allowed $`\delta m_{12}^2`$.
* Both the absolute scale of neutrino masses and the detailed pattern of the masses remain unknown, as the results in hand measure only mass differences (and leave the sign of $`\delta m_{23}^2`$ undetermined). The most stringent current bound on the absolute scale of neutrino mass comes from recent precision tests of cosmology (notably WMAP and the 2dF Galaxy Redshift Survey). This limit, a bound of about 1 eV for the sum of neutrino masses, is likely to improve as new surveys are done. In addition, much improved tritium $`\beta `$ decay and neutrinoless $`\beta \beta `$ decay experiments are being planned. There is one controversial claim of an observation of neutrinoless $`\beta \beta `$ decay that must be checked soon.

We stress, as we did in I, that this field is producing many new results that promise to impact physics broadly. The most common mechanisms for explaining neutrino mass suggest that current experiments are connected with phenomena far outside the standard model, residing near the Grand Unified energy scale of $`10^{16}`$ GeV. Thus there is hope that, by fully determining the properties of neutrinos – a few of the unresolved problems have been mentioned here – we may equip theorists to begin constructing the next standard model. Neutrino physics is also crucial to astrophysics – not just the standard solar model, but also in supernovae and in high-energy astrophysical environments – and to cosmology. We now have identified the first component of particle dark matter, though the significance of neutrino mass is still unclear due to our ignorance of the overall scale. Neutrinos could prove central to one of cosmology’s deepest questions, why our universe is matter dominated, rather than matter-antimatter symmetric. But this is a story for another paper and another time.

Acknowledgement

The work was supported in part by the U.S. Department of Energy and by the National Science Foundation under award PHY-98-01875.


text:
# Parallel netCDF: A Scientific High-Performance I/O Interface

## 1 Introduction

Scientists have recognized the importance of portable and efficient mechanisms for storing large datasets created and used by their applications. The Network Common Data Form (netCDF)  is one such mechanism used by a number of applications.

NetCDF intends to provide a common data access method for atmospheric science applications to deal with a variety of data types that encompass single-point observations, time series, regularly spaced grids, and satellite or radar images . Today several organizations have adopted netCDF as a data access standard .

The netCDF design consists of both a portable file format and an easy-to-use application programming interface (API) for storing and retrieving netCDF files across multiple platforms. More and more scientific applications choose netCDF as their output file format. While these applications become computational and data intensive, they tend to be parallelized on high-performance computers. Hence, it is highly desirable to have an efficient parallel programming interface to the netCDF files. Unfortunately, the original design of netCDF interface is proving inadequate for parallel applications because of its lack of a parallel access mechanism. In particular, there is no support for concurrently writing to a netCDF file. Hence, parallel applications operating on netCDF files must serialize access. Traditionally, parallel applications write to netCDF files through one of the allocated process which easily becomes a performance bottleneck. The serial I/O access is both slow and cumbersome to the application programmer.

To facilitate parallel I/O operations, we have defined a parallel API for concurrently accessing netCDF files. With minimum changes to the names and argument lists, this interface maintains the look and feel of the serial netCDF interface while the implementation underneath incorporates well-known parallel I/O techniques such as collective I/O to allow high-performance data access. We implement this work on top of MPI-IO, which is specified by MPI-2 standard  and is freely available on most platforms. Since MPI has become de facto parallel mechanism for communication and I/O on most parallel environments, this approach is portable across different platforms.

Hierarchical Data Format version 5 (HDF5)  also provides a portable file format and programming interfaces for storing multidimensional arrays together with ancillary data in a single file. It supports parallel I/O and its implementation is also built on top of MPI-IO. However, the HDF5 API is too flexible and cumbrous to become an easy-to-use standard since it adds more programming features and completely re-designs the API from its previous version. Our parallel netCDF interface, on the other hand, is more concise, closer to the original API, and goes much closer to MPI-IO interface, which introduces less overhead while providing more optimization opportunities for performance enhancement. The goal of this work is to make the parallel netCDF interface a data access standard for parallel scientific applications.

We run a couple of benchmarks using parallel netCDF and parallel HDF5, exploring both artificially made access patterns from our own benchmark and the ones from a real application called FLASH . In our experiments, parallel netCDF brings significant I/O improvement and shows better performance than parallel HDF5 in the FLASH I/O benchmark .

The rest of this paper is organized as follows. Section 2 reviews some related work. Section 3 presents the design background of netCDF and points out its potential usage in parallel scientific applications. The design and implementation of our parallel netCDF is described in Section 4. Experimental performance results are given in Section 5. Section 6 concludes the paper.

## 2 Related Work

Considerable research has been done on data access for scientific applications. The work has focused on data I/O performance and data management convenience. Two projects, MPI-IO and HDF, are most closely related to our research.

MPI-IO is a parallel I/O interface specified in the MPI-2 standard. It is implemented and used on a wide range of platforms. The most popular implementation, ROMIO  is implemented portably on top of an abstract I/O device layer  that enables portability to new underlying I/O systems. One of the most important features in ROMIO is collective I/O operations, which adopt a two-phase I/O strategy  and improve the parallel I/O performance by significantly reducing the number of I/O requests that would otherwise result in many small, noncontiguous I/O requests. However, MPI-IO reads and writes data in a raw format without providing any functionality to effectively manage the associated metadata. Nor does it guarantee data portability, thereby making it inconvenient for scientists to organize, transfer, and share their application data.

HDF is a file format and software, developed at NCSA, for storing, retrieving, analyzing, visualizing, and converting scientific data. The most popular versions of HDF are HDF4  and HDF5 . The design goal of HDF4 is mainly to deal with sequential data access and its APIs are consistent with its earlier versions. On the other hand, HDF5 is a major revision in which its APIs are completely re-designed. Both versions store multidimensional arrays together with ancillary data in portable, self-describing file formats. The support for parallel data access in HDF5 is built on top of MPI-IO, which ensures its portability since MPI-IO has become a de facto standard for parallel I/O. However, the fact that HDF5 file format is not compatible with HDF4 can be inconvenient for existing HDF4 programmers to migrate their applications to HDF5. Furthermore, HDF5 adds several new features, such as hierarchical file structure, to describe more metadata, but it also increases the difficulties for the implementation of parallel data access underneath. And the overhead involved may make HDF5 perform much worse than its underlying MPI-IO. By using a number of scientific applications, this problem is addressed in .

## 3 NetCDF Background

NetCDF is an abstraction that supports a view of data as a collection of self-describing, portable, array-oriented objects that can be accessed through a simple interface. It defines a file format as well as a set of programming interfaces for storing and retrieving data in the form of arrays in netCDF files. We first describe the netCDF file format and its serial API and then consider various approaches to access netCDF files in parallel computing environments.

### 3.1 File Format

NetCDF stores data in an array-oriented dataset, which contains dimensions, variables, and attributes. Physically, the dataset file is divided into two parts: file header and array data. The header contains all information (or metadata) about dimensions, attributes, and variables except for the variable data itself, while the data part contains arrays of variable values (or raw data).

The netCDF file header first defines a number of dimensions, each with a name and a length. These dimensions are used to define the shapes of variables in the dataset. One dimension can be unlimited and is used as the most significant dimension (record dimension) for growing-size variables.

Following the dimensions, a list of named attributes are used to describe the properties of the dataset (e.g., data range, purpose, associated applications ). These are called global attributes and are separate from attributes associated with individual variables.

The basic units of named data in a netCDF dataset are variables, which are multidimensional arrays. The header part describes each variable by its name, shape, named attributes, data type, array size, and data offset, while the data part stores the array values for one variable after another, in their defined order.

To support variable-size arrays (e.g., data growing with time stamps), netCDF introduces record variables and uses a special technique to store such data. All record variables share the same unlimited dimension as their most significant dimension and are expected to grow together along that dimension. The rest, less significant dimensions all together define the shape for one record of the variable. For fixed-size array, each array is stored in a contiguous file space starting from a given offset. For variable-size arrays, netCDF first defines a record of an array as a subarray comprising all fixed dimensions and the records of all such arrays are stored interleaved in the arrays’ defined order. Figure 1 illustrates the storage layouts for fixed-size and variable-size arrays in a netCDF file.

In order to achieve network transparency (machine-independence), both the header and data parts of the file are represented in an well-defined format similar to XDR (eXternal Data Representation) but extended to support efficient storage of arrays of non-byte data.

### 3.2 Serial NetCDF API

The original netCDF API was designed for serial codes to perform netCDF operations through a single process. In the serial netCDF library, a typical sequence of operations to write a new netCDF dataset is to create the dataset; define the dimensions, variables and attributes; write variable data; and close the dataset. Reading an existing netCDF dataset involves first opening the dataset; inquiring about dimensions, variables, and attributes; reading variable data; and closing the dataset.

These netCDF operations can be divided into the following five categories. Refer to  for details of each function in the netCDF library.

* Dataset Functions: create/open/close a dataset, set the dataset to define/data mode, and synchronize dataset changes to disk
* Define Mode Functions: define dataset dimensions and variables
* Attribute Functions: manage adding, changing, and reading attributes of datasets
* Inquiry Functions: return dataset metadata: dim(id, name, len), var(name, ndims, shape, id)
* Data Access Functions: provide the ability to read/write variable data in one of the five access methods: single value, whole array, subarray, subsampled array (strided subarray) and mapped strided subarray

The I/O implementation of the serial netCDF API is built on the native I/O system calls and has its own buffering mechanism in user space. Its design and optimization techniques are suitable for serial access but are not efficient or even not possible for parallel access, nor do they allow further performance gains provided by modern parallel I/O techniques.

### 3.3 Using NetCDF in Parallel Environments

Today most scientific applications are programmed to run in parallel environments due to the increasing requirements on data amount and computational resources. It is highly desirable to develop a set of parallel APIs for accessing netCDF files that employs appropriate parallel I/O techniques. In the meantime, programming convenience is also important, since scientific users may desire to spend minimum effort on dealing with I/O operations. Before presenting our design on parallel netCDF, we would like to discuss current approaches for using netCDF in parallel programs in a message-passing environment.

The first and most straightforward approach is described in the scenario of Figure 2(a) in which one process is in charge of collecting/distributing data and performing I/O to a single netCDF file using the serial netCDF API. The I/O requests from other processes are carried out by shipping all the data through this single process. The drawback of this approach is that collecting all I/O data on a single process can easily cause an I/O performance bottleneck and may overwhelm its memory capacity.

To avoid unnecessary data shipping, an alternative approach is to have all processes perform their I/O independently using the serial netCDF API, as shown in Figure 2(b). In this case, all netCDF operations can proceed concurrently, but over multiple files, one for each process. However, using multiple files to store a single netCDF dataset results the complexity and difficulty of data management. This approach also destructs the purpose of netCDF design on easy data integration and management.

A third approach introduces a new set of APIs with parallel access semantics and optimized parallel I/O implementation such that all processes perform I/O operations cooperatively or collectively through the parallel netCDF library to access a single netCDF file. This approach, as shown in Figure 2(c), both frees the users from dealing with details of parallel I/O and provides more opportunities for employing various parallel I/O optimizations in order to obtain higher performance. We discuss the details of this parallel netCDF design and implementation in the next section.

## 4 Parallel NetCDF

To facilitate convenient and high-performance parallel access to netCDF files, we define a new parallel interface and provide a prototype implementation. Since a large number of existing users are running their applications over netCDF, our parallel netCDF design retains the original netCDF file format (version 3) and introduces minimum changes from the original interface. We distinguish the parallel API from the original serial API by prefixing the C function calls with “ncmpi\_” and the Fortran function calls with “nfmpi\_”.

### 4.1 Interface Design

Our parallel netCDF API is built on top of MPI-IO. The parallel netCDF built on MPI-IO can benefit from several well-known optimizations already used in existing MPI-IO implementations, such as data sieving and two-phase I/O strategies  in ROMIO. Figure 3 describes the overall architecture for our design.

In parallel netCDF, a file is opened, operated, and closed by the participating processes in a communication group. In order for these processes to operate on the same file space, especially the structural information contained in the file header, a number of changes have been made to the original serial netCDF API.

For the function calls that create/open a netCDF file, an MPI communicator is added in the argument list to define the participating I/O processes within the file’s open and close scope. MPI\_Info object is also added to pass user access hints to the MPI-IO for further optimizations. By describing the collection of processes with a communicator, we provide the underlying implementation with information that can be used to ensure file consistency. The MPI\_Info hint provides users the ability to deliver the high level access information to netCDF and MPI-IO libraries, such as file access patterns and file system specifics to direct optimization.

We keep the same syntax and semantics for the parallel netCDF define mode functions, attribute functions, and inquiry functions as the original ones. These functions are also made collective to guarantee consistency of dataset structure among the participating processes in the same MPI communication group. For instance, the define mode functions is required to be called by the processes with the same values.

The major effort of this work is the parallelization of the data access functions. We provide two sets of data access APIs: a high-level API that mimics the serial netCDF data access functions and serves an easy path for original netCDF users to migrate to the parallel interface, and a flexible API that provides a more MPI-like style of access. Specifically, the flexible API uses more MPI functionality in order to provide better handling of internal data representations and to more fully expose the capabilities of MPI-IO to the application programmer. The major difference between the two is the use of MPI derived data types. We believe using MPI derived datatypes can better illustrate the access patterns than the subarray mapping methods used in original API.

The most important change from the original netCDF interface with respect to data access functions is the split of data mode into two distinct modes: collective and non-collective data modes in which collective function names end with “\_all”. Similar to MPI-IO, the collective functions are synchronous across the processes in the communicator associated to the opened netCDF file, while the non-collective functions are not. Using collective operations can provide the underlying parallel netCDF implementation an opportunity to further optimize access to the netCDF file. These optimizations are performed without further intervention by the application programmer and have been proven to provide dramatic performance improvement in multidimensional dataset access . Figure 4 gives an example code of using our parallel netCDF API to write and read a dataset using collective I/O.

### 4.2 Parallel Implementation

The parallel API implementation is discussed in two parts: header I/O and parallel data I/O. We first describe out implementation strategies for dataset functions, define mode functions, attribute functions, and inquiry functions that access the netCDF file header.

#### 4.2.1 Access to File Header

Internally, the header is read/written only by a single process, although a copy is cached in local memory on each process. The define mode functions, attribute functions, and inquiry functions all work on the local copy of the file header. Since they are all in-memory operations not involved in any file I/O, they bear few changes from the serial netCDF API. They are made collective, but this feature does not necessarily imply inter-process synchronization. In some cases, however, when the header definition is changed synchronization is needed to verify that the values passed in by all processes match. In all possible cases we allow inter-process communications.

The dataset functions, unlike the other functions cited, need complete reimplementation because they are in charge of collectively opening/creating datasets, performing header I/O and file synchronization for all processes, and managing inter-process communication. We build these functions over MPI-IO so that they have better portability and provide more optimization opportunities. The basic idea is to let the ROOT process fetch the file header, broadcast it to all processes when opening a file, and write the file header at the end of definition if any modification occurs in the header part. Since all define mode and attribute functions are collective and require all processes in the communicator to provide the same arguments when adding/removing/changing definitions, the local copies of the file header shall be the same across all processes once the file is collectively opened and until it is closed.

#### 4.2.2 Parallel I/O for Array Data

Since the majority of time spent accessing a netCDF file is in data access, the data I/O must be efficient. By implementing the data access functions above MPI-IO, we enable a number of advantages and optimizations.

For each of the five data access methods in the flexible data access functions, we represent the data access pattern as an MPI file view (a set of data visible and accessible from an open file ), which is constructed from the variable metadata (shape, size, offset, etc.) in the netCDF file header and start\[\], count\[\], stride\[\], imap\[\], mpi\_datatype arguments provided by users. For parallel access, particularly for collective access, each process has a different file view; and all processes in combination can make a single MPI-IO request to transfer large contiguous data as a whole, thereby preserving useful semantic information that would otherwise be lost if the transfer were expressed as per process noncontiguous requests.

The high-level data access functions are implemented in terms of the flexible data access functions, so that existing users migrating from serial netCDF can also benefit from the MPI-IO optimizations. However, the flexible data access functions are closer to MPI-IO and hence incur less overhead. They accept a user-specified MPI derived datatype and pass it directly to MPI-IO for optimal handling of in-memory data access patterns.

In some cases (for instance, in record variable access) the data is stored interleaved by record and the contiguity information is lost, so the existing MPI-IO collective I/O optimization may not help. In that case, we need more optimization information from users, such as the number, order, and record indices of the record variables they will access consecutively. With such information we can collect multiple I/O requests over a number of record variables and optimize the file I/O over a large pool of data transfers, thereby producing more contiguous and larger transfers. This kind of information is passed in as an MPI\_Info hint when a user opens or creates a netCDF dataset. We implement our user hints in parallel netCDF for all such specific optimization points, while a number of standard hints are passed down for MPI-IO to take control of optimal parallel I/O behaviors. Thus experienced users have the opportunity to tune their applications for further performance gains.

### 4.3 Advantages and Disadvantages

There are a number of advantages within the design and implementation of our parallel netCDF, as compared to other related work, like HDF5.

First of all, the parallel netCDF design and implementation is optimized for the netCDF file format so that the data I/O performance is as good as the MPI-IO. The NetCDF file chooses linear data layout, in which the data arrays are either stored in contiguous space and in a predefined order or interleaved in a regular pattern. This regular and highly predictable data layout enables the parallel netCDF data I/O implementation to simply pass the data buffer, metadata (fileview, mpi\_datatype, etc.), and other optimization information to MPI-IO, and all parallel I/O operations are carried out in the same manner as when MPI-IO alone is used. Thus, there is very little overhead, and the parallel netCDF performance should be nearly the same as MPI-IO if only raw data I/O performance is compared. On the other hand, parallel HDF5 uses tree-like file structure that are similar to the UNIX file system and the data is dispersedly laid out using super block, header blocks, data blocks, extended header blocks and extended data blocks. This irregular layout pattern may make it difficult to pass user access pattern directly to MPI-IO especially for the case of variable-size arrays. Instead, parallel HDF5 uses dataspace and hyperslabs to define the data organization, map and transfer data between memory space and the file space and does buffer packing/unpacking in a recursive way, while these can otherwise be directly handled by MPI-IO in a more efficient and optimized way.

Secondly, the parallel netCDF implementation manages to keep the overhead involved in header I/O as low as possible. In the netCDF file, there is only one header which contains all necessary information for direct access of each data array and each array is associated with a predefined, numerical ID that can be efficiently inquired when it is needed to access the array. So, by maintaining a local copy of the header on each process, our implementation saves a lot of inter-process synchronization as well as avoids repeated access of the file header each time the header information is needed to access a single array. All header information can be accessed directly in local memory and inter-process synchronization is needed only during the definition of the dataset. And once the definition of the dataset is created, each array can be identified by its permanent ID and accessed at any time by any process, without any collective open/close operation. However, in HDF5, the header metadata is dispersed in separate header blocks for each object and, in order to operate on an object, it has to iterate through the entire namespace to get the header information of that object and then open, access and close it. This kind of access method may be inefficient for parallel access, since the parallel HDF5 designs the open/close of each object as collective operations, which force all participating processes to communicate when accessing one single object, not to mention the cost of file access to locate and fetch the header information of that object.

Lastly, the programming interface of the parallel netCDF is concise and designed for easy usage, and the file format is fully compatible with serial netCDF. Porting existing serial netCDF application to parallel netCDF should be straightforward because the parallel API contains nearly all functions of the serial API with parallel semantics but with minimum change of function names and argument lists.

However, there are also limitations in parallel netCDF. Unlike HDF5, netCDF does not support hierarchical group based organization of data objects and since it lays out the data in a linear order, adding fixed-size array or extending the file header may be very costly once the file is created and has existing data stored, though moving the existing data to extended area is performed in parallel. Also, parallel netCDF does not provide functionality to combine two or more files in memory through software mounting, as HDF5 does. Nor does netCDF support data compression within its file format. Fortunately, these features can all be achieved by external software, sacrificing some manageability of the files.

## 5 Performance Evaluation

To evaluate the performance and scalability of our parallel netCDF with that of serial netCDF, we ran some experiments and compared the results. We also compared the performance of parallel netCDF with that of parallel HDF5, using the FLASH I/O benchmark.

The experiments were run on an IBM SP-2 machine. This system is a teraflop-scale clustered SMP with 144 compute nodes. Each compute node has 4 GB of memory shared among its eight 375 MHz Power3 processors. All the compute nodes are interconnected by switches and also connected via switches to the multiple I/O nodes running the GPFS parallel file system. There are 12 I/O nodes, each with dual 222 MHz processes. The aggregate disk space is 5 TB and the peak I/O bandwidth is 1.5 GB/s.

### 5.1 Scalability Analysis

We wrote a test code (in C language) to evaluate the performance of the current implementation of parallel netCDF. This test code was originally developed in Fortran by Woo-sun Yang and Chris Ding at Lawrence Berkeley National Laboratory (LBL). Basically it reads/writes a three-dimensional array field tt(Z,Y,X) from/into a single netCDF file, where Z=level is the most significant dimension and X=longitude is the least significant dimension. The test code partitions the three dimensional array along Z, Y, X, ZY, ZX, YX, and ZYX axes, respectively, as illustrated in Figure 5. All data I/O operations in these tests used collective I/O. For comparison purpose, we prepared the same test using the original serial netCDF API and ran it in serial mode, in which a single processor reads/writes the whole array.

Figure 6 shows the performance results for reading and writing 64 MB and 1 GB netCDF datasets. Generally, the parallel netCDF performance scales with the number of processes. Because of collective I/O optimization, the performance difference made by various access patterns is small, although partitioning in the Z dimension generally performs better than in the X dimension because of the different access contiguity. The overhead involved is inter-process communication, which is negligible comparing to the disk I/O when using large file size. The I/O bandwidth does not scale in direct proportion because the number of I/O nodes (and disks) is fixed so that the dominating disk access time at I/O nodes is almost fixed. As expected, the parallel netCDF outperforms the original serial netCDF as the number of processes increases. The difference between the serial netCDF performance and the parallel netCDF performance with single processor is because of their different I/O implementations and different I/O caching/buffering strategies. In the serial netCDF case, if, as in Figure 2(a), multi-processors were used and the ROOT processor needed to collect partitioned data and then perform the serial netCDF I/O, the performance would be much worse and decrease with the number of processors because of the additional communication cost and division of a large I/O request into a series of small requests.

### 5.2 FLASH I/O Performance

The FLASH I/O benchmark simulates the I/O pattern of an important scientific application called FLASH . It recreates the primary data structures in the FLASH code and produces a checkpoint file, a plotfile with centered data, and a plotfile with corner data, using parallel HDF5. Basically, these three output files contains a series of multidimensional arrays, and the access pattern is simple (Block, \*, …), which is similar to the Z partition in Figure 5. In each of the files, the benchmark writes the related arrays in a fixed order from contiguous user buffers, respectively. The I/O routines in the benchmark are identical to the routines used by FLASH, so any performance improvements made to the benchmark program will be shared by FLASH. In our experiments, in order to focus on the data I/O performance, we modified this benchmark, removed the part of code writing attributes, ported it to parallel netCDF, and observed the effect of our new parallel I/O approach.

Figure 7 shows the performance results of the FLASH I/O benchmark using parallel netCDF and parallel HDF5. We tested both small data size and large data size. The parameters used in these two experiments are: (a) nxb = nyb = nzb = 8, nguard = 4, number of blocks = 80, and nvar = 24; (b) nxb = nyb = nzb = 16, nguard = 8, number of blocks = 80, and nvar = 24. Although both I/O libraries are built above MPI-IO, the parallel netCDF has much less overhead and outperforms parallel HDF5 by almost doubling the overall I/O rate. The extra overhead involved in the current release of HDF5 (version 5-1.4.3) includes inter-process synchronizations and file header access performed internally in parallel open/close of every dataset (analogous to a netCDF variable) and recursive handling of the hyperslab used for parallel access, which makes the packing of the hyperslabs into contiguous buffers take a relatively long time.

## 6 Conclusion and Future Work

In this work, we extend the serial netCDF interface to facilitate parallel access, and we provide an implementation for a subset of this new parallel netCDF interface. By building on top of MPI-IO, we gain a number of interface advantages and performance optimizations users can benefit from by using this parallel netCDF package, as shown by our test results. So far, a number of users from LBL, ORNL, and University of Chicago are using our parallel netCDF library.

Future work involves developing a production-quality parallel netCDF API (for C, C++, Fortran, and other programming languages) and making it freely available to the high-performance computing community. Moreover, we need to develop a mechanism for matching the file organization to access patterns, and we need to develop cross-file optimizations for addressing common data access patterns.

### Acknowledgements

This work is sponsored by Scientific Data Management Center of DOE SciDAC ISICs program and jointly conducted at Northwestern University and Argonne National Laboratory. This research was also supported in part by NSF cooperative agreement ACI-9619020 through computing resources provided by the National Partnership for Advanced Computational Infrastructure at the San Diego Supercomputer Center.

We thank Woo-Sun Yang from LBL for providing us the test code for performance evaluation and Nagiza F. Samatova and David Bauer at ORNL for using our library and for giving us feedback and valuable suggestions.


text:
# Orbital Motion in the Radio Galaxy 3C 66B: Evidence for a Supermassive Black Hole Binary

> Supermassive black hole binaries may exist in the centers of active galactic nuclei like quasars and radio galaxies and mergers between galaxies may result in the formation of supermassive binaries during the course of galactic evolution. Using the very-long-baseline interferometer, we imaged the radio galaxy 3C 66B at radio frequencies and found that the unresolved radio core of 3C 66B shows well-defined elliptical motions with a period of 1.05 $`\pm `$ 0.03 years, which provides a direct detection of a supermassive black hole binary.

The presence of a supermassive black hole binary (or a supermassive binary; SMB) in active galactic nuclei (AGNs) has been suggested mainly by periodic optical and radio outbursts (1,2), wiggled patterns of compact radio jets (i.e., the indication of precession motions of the radio jets) (3,4), and X-shaped morphology of radio lobes (5). AGNs with those characteristics tend to be associated with strong radio jets and thus their host galaxies are often giant elliptical galaxies (6). If they were made from mergers between or among nucleated galaxies (7,8), they may contain two or more black holes in their central regions because progenitor galaxies may contain a supermassive black hole in their nuclei (9,10). Captured black holes will be settled under dynamical friction in the core of the merged stellar system and then will form a SMB (11). Therefore, it is important to establish the presence of a SMB in an AGN and then investigate the true role of a SMB in AGN activity.

The most direct evidence for a SMB can be obtained by detection of the Kepler orbital motion of some emission component close to black holes. One technical problem is that an expected separation between the two supermassive black holes (e.g., $`10^{17}`$ cm) is too small to be resolved spatially when we use optical telescope facilities because that separation corresponds to $``$100 micro-arcsecond ($`\mu `$as) if its host galaxy is located at a distance of 100 Mparsec (Mpc). However, the use of a technique of phase-referencing very-long-baseline interferometry (VLBI) at radio frequencies allows us to achieve such position measurements with the accuracy of tens of $`\mu `$as (12). We looked for the Kepler motion of a radio-emission component in a pair of radio sources, 3C 66B and 3C 66A; 3C 66B is a radio galaxy at redshift $`z=0.0215`$ (13), whereas 3C 66A is a more distant BL Lac object at $`z=0.44`$ (14). This pairing allows us to use this source as the stationary position reference to 3C 66B. They are one of best pairs for the phase-referencing VLBI, because they are bright radio sources and their spatial separation is only 6 arc min.

We observed this pair at 2.3 and 8.4 GHz with the very-long-baseline array (VLBA) of the National Radio Astronomy Observatory (NRAO) over six epochs between 13 March 2001 and 14 June 2002. Because the separation between 3C 66A and 3C 66B is smaller than the beam widths of the VLBA antennas at 2.3 and 8.4 GHz (26 and 6 arc min, respectively), it is possible to observe both sources in each antenna beam simultaneously at these frequencies.

The radio core position of 3C 66B at each epoch was measured with respect to that of 3C 66A. We used the local brightest peak in the map, rather than the Gaussian peak, to determine the radio core position. The Gaussian peak may have a larger uncertainty as a result of the effects of asymmetrical structure of the sources.

Time variations of the radio core position in 3C 66B at 2.3 and 8.4 GHz (Fig. 1 and Fig. 2, respectively) can be fitted by an elliptical motion (Table 1). The reduced chi-square value obtained when calculated with the observation error is $``$ 0.6 at 2.3 GHz and $``$ 0.05 at 8.4 GHz (12). Because these values are less than unity, especially at 8.4 GHz, the position error might have been overestimated. Assuming the reduced chi-square value of unity, we find that the revised position errors, $`\mathrm{\Delta }\alpha ^{}`$ and $`\mathrm{\Delta }\delta ^{}`$, are 71 and 60 $`\mu `$as at 2.3 GHz and 9 and 7 $`\mu `$as at 8.4 GHz, respectively (12). Although the orbital periods at 2.3 and 8.4 GHz are nearly the same, the major axis at 2.3 GHz is about five times as long as that at 8.4 GHz and the position angle of the major axis of the two frequencies differs by $`24^{}`$.

The averaged period of the fitted core motion is estimated to be $`1.05\pm 0.03`$ years. Although this period might be interpreted as due to some geodetic effects related to Earth’s orbital motion around the sun, such possibilities can be rejected for the following reasons. The annual parallax of 3C 66B cannot explain the observed motion, because the distance of 3C 66B (85 Mpc for the Hubble constant $`H_0=75`$ km s-1 Mpc-1) leads to the amplitude of the annual parallax of only $`0.01`$ $`\mu `$as. Furthermore, the large difference in the length of the major axis between 2.3 and 8.4 GHz can not be understood from our knowledge of the annual parallax. Only a possible variability of the total electron content (TEC) in the ionosphere may explain this frequency dependence, because the excess path due to the ionosphere is proportional to TEC $`\nu ^2`$, where $`\nu `$ is the observing frequency (15). The typical position error due to the ionosphere for the separation of 6 arc min is estimated to be 0.005 $`\left(\frac{\nu }{\mathrm{GHz}}\right)^2`$ $`\mu `$as, which is negligible (16). Another possibility is that the radio core in 3C 66B is subject to microlensing effects by a star in the Milky Way galaxy near 3C 66B on the celestial plane (17) and is moved by the annual parallax of the star. We attempted to fit the radio core motion by this microlensing model and obtained the reduced chi-square value of 1.6, corresponding to a chi-square probability of $`10`$ %. This result indicates that the microlensing model does not explain the sinusoidal core motion as well as it explains the elliptical motion of the radio core.

Because the radio core is located at the root of the jet where the optical depth is unity (18, 19), the radio core at 2.3 GHz is expected to be located at a greater distance from the central engine than that at 8.4 GHz. Thus, a precession motion of the jet provides a natural explanation for the observation that the major axis of the radio core orbit is longer at 2.3 GHz than at 8.4 GHz. Although the orbital major axis should be always perpendicular to the mean jet axis, the observed inclination angle between the mean jet axis and the orbital major axis is estimated to be 63$`\pm 9^{}`$ at 2.3 GHz and 39$`\pm 6^{}`$ at 8.4 GHz (calculated on the basis of a $`50^{}`$ position angle of the jet, which is estimated from the global jet structure at 2.3 GHz). This fact indicates that the observed core motion is not dominated by only the simple precession motion. The most plausible explanation for this result is that the elliptical motion of the radio core is attributable to the combination of the orbital and precession motions of a SMB in 3C 66B, whereas the core motion at 2.3 GHz is dominated by the precession motion of the jet, which at 8.4 GHz is dominated by the orbital motion of the SMB (Fig. 3).

According to Kepler’s third law, the observed period can be expressed as,

$$P=2\pi G^{1/2}r^{3/2}(M+m)^{1/2},$$
(1)
where where $`G`$ is the gravitation constant, $`r`$ is the separation of two orbiting black hole, and $`M`$ and $`m`$ ($`Mm`$) are the masses of the black holes, respectively. Eq. 1 allows us to determine the mass density $`\rho `$ as a function of only $`P`$,

$$\rho \frac{M+m}{(4/3)\pi r^3}=\frac{3\pi }{GP^2}.$$
(2)
The observed result leads to $`\rho 2\times 10^{15}`$ $`M_{}`$pc-3 (where $`M_{}`$ is the mass of the Sun), which is larger than those of massive dark objects in NGC 4258 and Sgr A (20, 21). Because the radio core motion at 8.4 GHz is believed to reflect the orbital motion of the SMB more than does that at 2.3 GHz, we regard the length of its major axis as the upper limit of the orbital radius of the SMB. Assuming that the less massive black hole emanates the jet, we obtain the upper limit of the separation and the mass of the SMB of $`r_{\mathrm{max}}5.4(1+q)\times 10^{16}`$ cm and $`M_{\mathrm{max}}4.4(1+q)^2\times 10^{10}M_{}`$, respectively, where $`q\frac{m}{M}`$ is the mass ratio between the objects. We note that these parameters are similar to those of a radio-loud AGN, OJ 287, which is believed to be an archetypical AGN with a SMB (22). This conclusion is consistent with 3C 66B’s being a giant elliptical galaxy which suggests a galaxy merger.

The SMB in 3C 66B will merge into one as a result of the gravitational radiation loss. The lifetime in this evolution phase can be estimated by the relation in Eq. 3 (23).

$$t_{\mathrm{GR}}7.2\times 10^4\left(\frac{M}{10^8M_{}}\right)^3\left(\frac{r}{10^{16}\mathrm{cm}}\right)^4\frac{1}{q(1+q)}\mathrm{yr}.$$
(3)
If we adopt $`q=0.1`$, we obtain $`r_{\mathrm{max}}=5.9\times 10^{16}`$ cm and $`M_{\mathrm{max}}=5.4\times 10^{10}M_{}`$. These values give us $`t_{\mathrm{GR}}5`$ yr for the SMB in 3C 66B. This time scale seems to be extremely short compared with a whole life time of a SMB system (13). A more reasonable lifetime would be obtained for smaller values of $`q`$ and/or $`r`$. Combining Eq. 1 and Eq. 3, we find $`t_{\mathrm{GR}}r^5`$. If the true $`r`$ is a quarter of the corresponding value for the upper limit, a lifetime $`10^3`$ times as long is obtained. The expected maximum amplitude of gravitational radiation from 3C 66B will be $`h10^{11}`$ just before the two black holes merge into one (24, 25).

> References and Notes
>
> 1. A. Sillanp$`\ddot{\mathrm{a}}\ddot{\mathrm{a}}`$, S. Haarala, M. J. Valtonen, B. Sundelius, G. G. Byrd, Astrophys. J. 325, 628 (1988).
> 2. E. Valtaoja, et al., Astrophys. J. 531, 744 (2000).
> 3. N. Roos, J. S. Kaastra, C. A. Hummel, Astrophys. J. 409, 130 (1993).
> 4. Z. Abraham, Astron. Astrophys. 355, 915 (2000).
> 5. D. Merritt, R. D. Ekers, Science 297, 1310 (2002).
> 6. C. M. Gaskell, in Jets from Stars and Galactic Nuclei, W. Kundt, Ed. (Springer-Verlag, New York, 2002), p 165
> 7. D. B. Sanders, et al., Astrophys. J. 325, 74 (1988).
> 8. Y. Taniguchi, S. Ikeuchi, Y. Shioya, Astrophys. J. 514, L9 (1999).
> 9. J. Kormendy, D. Richstone, Annu. Rev. Astron. Astrophys. 33, 581 (1995).
> 10. M. J. Rees, Annu. Rev. Astron. Astrophys. 22, 471 (1984).
> 11. M. C. Begelman, R. D. Blandford, M. J. Rees, Nature 287, 307 (1980).
> 12. Materials and Methods are available as supporting material in Science Online.
> 13. T. A. Matthews, W. W. Morgan, M. Schmidt, Astrophys. J. 140, 35 (1964).
> 14. H. R. Miller, B. Q. McGimsey, Astrophys. J. 220, 19 (1978).
> 15. A. R. Thompson, J. M. Moran, G. W. Swenson, in Interferometry and Synthesis in Radio Astronomy, (Wiley, New York, 1986) p 554
> 16. The excess path in the zenith direction due to ionosphere is $`l(\nu )=40.3\mathrm{TEC}\nu ^2`$ (15). As a result of the Doppler shift by the Earth’s spin, a frequency shift between two sources occurs, given by $`\nu ^{}=\nu (1+\omega \theta D/c)`$, where $`D`$ is the baseline length, and $`\omega `$ is the angular velocity of the Earth spin. Thus, dividing $`l(\nu )l(\nu ^{})`$ by $`c`$, we obtain the time delay due to ionospheric excess path. We assumed that TEC in the zenith direction is $`6\times 10^{17}`$ m-2.
> 17. S. Refsdal, Mon. Not. R. Astron. Soc. 128, 295 (1964).
> 18. A. K$`\ddot{\mathrm{o}}`$nigl, Astrophys. J. 243, 700 (1981).
> 19. A. P. Lobanov, Astron. Astrophys. 330, 79 (1998).
> 20. M. Miyoshi, et al., Nature 373, 127 (1995).
> 21. A. M. Ghez, M. Morris, E. E. Becklin, A. Tanner, T. Kremenek, Nature 407, 349 (2000).
> 22. H. J. Lehto, M. J. Valtonen, Astrophys. J. 460, 207 (1996)
> 23. 23. Q. Yu, Mon. Not. R. Astron. Soc. 331, 935 (2002).
> 24. K. P. Thorne, V. B. Braginsky, Astrophys. J. 204, L1 (1976).
> 25. T. Fukushige, T. Ebisuzaki, Astrophys. J. 396, L61 (1992)
> 26. Supported by a grant-in-aid for Fellows of the Japan Society for the Promotion of Science from the Japanese Ministry of Education, Culture, Sports, and Science. NRAO is a facility of the National Science Foundation, operated under cooperative agreement by Associated Universities, Inc. We acknowledge anonymous reviewers for useful suggestions.

Fig. 1. Orbital fit calculations applied to the data of the position change of the core at 2.3 GHz. Observations were carried out on 13 March, 25 June, and 9 November 2001, and 8 February, 21 February, and 14 June 2002. We show (A) spatial distribution, (B) time evolution toward right ascension direction, and (C) time evolution toward declination direction. The error bars indicate the revised position errors suggested from the chi-square value; 71 and 60 $`\mu `$as at 2.3 GHz and 9 and 7 $`\mu `$as at 8.4 GHz, in right ascension and declination, respectively (12).

Fig. 2 Same as Fig. 1, but at 8.4 GHz.

Fig. 3 A schematic view of the geometrical relation between the radio core motion and the SMB. Because the orbital velocity of the SMB is added to the intrinsic jet velocity, the jet precession can be induced with the same period as that of the SMB orbital motion (3). Therefore, if the radio core is far enough from the SMB, the orbital major axis of the radio core should be observed as perpendicular to the mean jet axis, as is the normal precession. However, if the radio core is very close to the SMB, it should be observed as almost parallel to the orbital plane of the SMB. Therefore, the position angle of the radio core orbit depends on the distance between the radio core and the SMB, if the jet direction misaligns with the angular momentum vector of the SMB orbit.


text:
# Transition to Instability in a Kicked Bose-Einstein Condensate

## Abstract

A periodically kicked ring of a Bose-Einstein condensate is considered as a nonlinear generalization of the quantum kicked rotor. For weak interactions between atoms, periodic motion (anti-resonance) becomes quasiperiodic (quantum beating) but remains stable. There exists a critical strength of interactions beyond which quasiperiodic motion becomes chaotic, resulting in an instability of the condensate manifested by exponential growth in the number of noncondensed atoms. Similar behavior is observed for dynamically localized states (essentially quasiperiodic motions), where stability remains for weak interactions but is destroyed by strong interactions.

The classical kicked rotor is a textbook paradigm for dynamical chaos reichl . The quantum kicked rotor has play an equally important role for the study of quantum chaos, for which a wide range of effects have been predicted QKR  and observed in experiments raizen . In recent years, the realization of Bose-Einstein condensation (BEC) of dilute gases BEC  has opened new opportunities for studying dynamical systems in the presence of many-body interactions. A natural question to ask is how the physics of the quantum kicked rotor is modified by interactions NLSE .

In the classical regime, chaotic motion leads to diffusive growth in the kinetic energy. In quantum dynamics, where chaos is no longer possible because of the linearity of the Schrödinger equation, the motion becomes periodic (anti-resonance), quasiperiodic (dynamical localization), or resonant (quantum resonance) hog ; anti-res . In the mean field approximation, many-body interactions in BEC are represented by adding a nonlinear term in the Schrödinger equation review . This nonlinearity makes it possible to bring chaos back into the system, leading to instability (in the sense of exponential sensitivity to initial conditions) of the condensate wave function instability . The onset of instability of the condensate can cause rapid proliferation of thermal particles castin  that can be observed in experiments ketterle . It is therefore important to understand the route to chaos with increasing interactions. This problem has recently been studied for the kicked BEC in a harmonic oscillator zoller .

In this Letter, we investigate the quantum dynamics of a BEC with repulsive interaction that is confined on a ring and kicked periodically. This system is a nonlinear generalization of the quantum kicked rotor. From the point of view of dynamical theory, the kicked rotor is more generic than the kicked harmonic oscillator, because it is a typical low dimensional system that obeys the KAM theorem liu1 . It is very interesting to understand how both quantum mechanics and mean field interaction affect the dynamics of such a generic system. We will focus most of our attention on the case of anti-resonance because it is the simplest periodic motion. Here we find, with both analytic and numerical calculations, that weak interactions will make the motion quasiperiodic in the form of quantum beating. For strong interactions, quasiperiodic motions are destroyed and we observe a transition to instability of the BEC that is also characterized by an exponential growth in the number of noncondensed atoms. Universal critical behavior for the transition is found. We have also studied nonlinear effect on dynamically localized states that may be regarded as quasiperiodic. Similar results are obtained in that localization remians for sufficiently weak interactions but become unstable beyond a critical strength of interactions.

Consider condensed atoms confined in a toroidal trap of radius $`R`$ and thickness $`r`$, where $`rR`$ so that lateral motion is negligible and the system is essentially one-dimensional ring . The dynamics of the BEC is described by the dimensionless nonlinear Gross-Pitaveskii (GP) equation,

$`i{\displaystyle \frac{}{t}}\psi `$ $`=`$ $`\left({\displaystyle \frac{1}{2}}{\displaystyle \frac{^2}{\theta ^2}}+g\left|\psi \right|^2+K\mathrm{cos}\theta \delta _t(T)\right)\psi ,`$ (1)
where $`g=8NaR/r^2`$ is the scaled strength of nonlinear interaction, $`N`$ is the number of atoms, $`a`$ is the $`s`$-wave scattering length, $`K`$ is the kick strength, $`\delta _t(T)`$ represents $`\underset{n}{}\delta \left(tnT\right)`$, $`T`$ is the kick period, and $`\theta `$ denotes the azimuthal angle. The length and the energy are measured in units $`R`$ and $`\frac{\mathrm{}^2}{mR^2}`$, respectively. The wavefunction $`\psi (\theta ,t)`$ has the normalization $`_0^{2\pi }\left|\psi \right|^2𝑑\theta =1`$ and satisfies periodic boundary condition $`\psi (\theta ,t)=\psi (\theta +2\pi ,t)`$. Experimentally, the ring-shape potential may be realized using two 2D circular “optical billiards” with the lateral dimension being confined by two plane optical billiards billiard . The $`\delta `$-kick may be realized by adding potential points along the ring with an off-resonant laser raizen . The interaction strength $`g`$ may be adjusted using a magnetic field-dependent Feshbach resonance Feshbach .

The mean energy of each particle is $`E(t)=_0^{2\pi }𝑑\theta \psi ^{}\left(\frac{1}{2}\frac{^2}{\theta ^2}+\frac{1}{2}g\left|\psi \right|^2\right)\psi `$. To determine the evolution of the energy, we numerically integrate Eq.(1) over a time span of $`100`$ kicks, using a split-operator method split , with the initial wavefunction $`\psi `$ being the ground state $`\psi (\theta ,0)=1/\sqrt{2\pi }`$. The kick period is chosen as $`T=2\pi `$ to match the condition for anti-resonance. After each kick, the energy $`E(t)`$ is calculated and plotted as shown in Fig.1. We see a remarkable difference among noninteraction (Fig.1(a)), weak interaction (Fig.1(b)), and strong interaction (Fig.1(c)) cases. For noninteraction, the energy $`E(t)`$ oscillates between two values (anti-resonance) and the period is $`2T`$; While in the weak interaction case, the amplitude of the oscillation decreases gradually to zero and then revives, similar to the phenomena of beating in classical waves. The values of the oscillation and beat frequencies are obtained by Fourier Transform and the results are presented in Fig.2. It is shown that, besides the appearance of the beat frequency, the interaction also modifies the oscillation frequency; For stronger interaction (Fig.1(c)), i.e. $`g1.96`$, we find that the energy’s evolution demonstrates an irregular pattern, clearly indicating the collapse of the quasiperiodic motion and the occurrence of instability.

For the case of weak interaction, the system is predominantly in the two energy levels. We can map the system onto a spin model with a two-mode approximation to the GP equation liu . By considering the conservation of parity we may write the wavefunction as $`\psi =\frac{1}{\sqrt{2\pi }}\left(a+\sqrt{2}b\mathrm{cos}\theta \right)`$, where the populations $`a`$ and $`b`$ at the ground and excited states satisfy the normalization condition $`\left|a\right|^2+\left|b\right|^2=1.`$ The Hamiltonian in the spin representation reads

$``$ $`=`$ $`{\displaystyle \frac{S_z}{2}}+{\displaystyle \frac{g}{2\pi }}\left(S_x^2{\displaystyle \frac{S_z}{4}}+{\displaystyle \frac{S_z^2}{8}}\right)+\sqrt{2}KS_x\delta _t(T),`$ (2)
where $`S_z`$ corresponds to the population difference $`\left|a\right|^2\left|b\right|^2`$ and arctan$`(S_y/S_x)`$ gives the relative phase $`\alpha =\mathrm{arg}(a)\mathrm{arg}(b)`$. This Hamiltonian is similar to a kicked top model kicktop , but here the evolution between two kicks is more complicated.

With the spin model, we can readily study the dynamics of the system. For the case of noninteraction, the evolution between two consecutive kicks is simply an angle $`\pi `$ rotation about the $`z`$ axis. The spin initially directing to north pole stays there for time duration $`T`$, then the first kick rotates the spin by an angle $`\sqrt{2}K`$ about the $`x`$ axis. The following free evolution rotates the spin by an angle $`\pi `$ about the $`z`$ axis. Then, the second kick will drive the spin back to north pole through another rotation of $`\sqrt{2}K`$ about the $`x`$ axis. With this the spin’s motion is two kick period recurrence and quantum anti-resonance occurs.

With interaction, the motion between two consecutive kicks is approximately described by a rotation of $`\pi +g(1+3S_z)/4`$ about the $`z`$ axis. Compared with the noninteraction case, the mean field interaction leads to an additional phase shift $`g(1+3S_z)/4`$. This phase shift results in a deviation of the spin from $`S_x=0`$ plane at time $`2T^{}`$, i.e., moment just before the second kick. As a result, the second kick cannot drive the spin back to its initial position and quantum anti-resonance is absent. However, the phase shift will be accumulated in future evolution and the spin may reach $`S_x=0`$ plane at a certain time $`mT^{}`$(beat period) when the total accumulated phase shift is $`\pi /2`$. Then the next kick will drive spin back north pole by applying an angle $`\sqrt{2}K`$ rotation about the $`x`$ axis.

The above picture is confirmed by our numerical solution of the spin Hamiltonian with fourth-order Runge-Kutta method recipe . In Fig.3 we see that, the relative phase at the moment just before the even kicks increases almost linearly and reaches $`2\pi `$ in a beat period. The slope of the increment reads, $`\gamma _{RP}=\left(\alpha \left(4T^{}\right)\alpha \left(2T^{}\right)\right)/2`$, which can be deduced analytically. With this and through a lengthy deduction, we obtain an analytic expression for the beat frequency to first order in $`g`$,

$`f_{beat}`$ $``$ $`{\displaystyle \frac{g}{4\pi }}\left(1+3\mathrm{cos}\left(\sqrt{2}K\right)\right).`$ (3)
This expression as well as the relation between the oscillation frequency and the beat frequency, $`f_{osc}=\frac{1}{2}\frac{1}{2}f_{beat}`$, is in very good agreement with the numerical results as shown in Fig.2. Therefore the beating provides a method to measure interaction strength in an experiment.

Tuning the interaction strength still larger means enhancing further the nonlinearity of the system. From our general understanding of nonlinear systems, we expect that the solution will be driven towards chaos, in the sense of exponential sensitivity to initial condition and random evolution in the temporal domain. The latter character has been clearly displayed by the irregular pattern of the energy evolution in Fig.1(c). On the other hand, the onset of instability (or chaotic motion) of the condensate is accompanied with the rapid proliferation of thermal particles. Within the formalism of Castin and Dum castin , the growth of the number of the noncondensed atom will be exponential, similar to the exponential divergence of nearby trajectories in phase space of classical system. The growth rate of the noncondensed atoms is similar to the Lyapounov exponent, turning from zero to nonzero as instability occurs.

In Castin and Dum’s formalism, the mean number of noncondensed atoms at zero temperature is described by $`\delta \widehat{N}(t)=_{k=1}^{\mathrm{}}v_k(t)|v_k(t)`$, where $`|v_k(t)`$ are governed by

$$i\frac{d}{dt}(\begin{array}{c}u_k\\ v_k\end{array})=(\begin{array}{cc}H+gQ|\psi |^2Q& gQ\psi ^2Q^{}\\ gQ^{}\psi ^2Q& HgQ^{}|\psi |^2Q^{}\end{array})(\begin{array}{c}u_k\\ v_k\end{array}),$$
(4)
where $`H=\frac{\widehat{p}^2}{2}+g|\psi |^2\mu ,`$ $`\mu `$ is the chemical potential, $`\psi `$ is the ground state of GP equation and the projection operators $`Q`$ are given by $`Q=1|\psi \psi |.`$

We numerically integrate Eq.(4) for the $`u_k`$, $`v_k`$ pairs over a time span of 100 kicks, using a split operator method, parallel to numerical integration of GP equation. The initial conditions $`|u_k(0)`$, $`|v_k(0),`$ for initial ground state wavefunction $`\psi (\theta )=1/\sqrt{2\pi },`$ are obtained by diagonalizing the linear operator in Eq.(4) castin2 . After each kick the mean number of noncondensed atoms is calculated and plotted versus time in Fig.4(a). We find that there exists a critical value for the interaction strength, i.e., $`g_c=1.96`$, above which, the mean number of noncondensed atoms increases exponentially, indicating the instability of BEC. Below the critical point, the mean number of noncondensed atoms increases polynomially. As the nonlinear parameter crosses over the critical point, the growth rate turns from zero to nonzero, following a square-root law (inset in Fig.4(a)). This scaling law may be universal for Bogoliubov excitation as confirmed by recent experiments ketterle .

The critical value of the interaction strength depends on the kick strength. For very small kick strength, the critical interaction is expected to be large, because the ground state of the ring-shape BEC with repulsive interaction is dynamically stable biao . For large kick strength, to induce chaos, the interaction strength must be large enough to compete with the external kick potential. So, in the parameter plane of $`(g,k)`$, the boundary of instability shows a ”U” type curve (Fig.4(b)).

Across the critical point, the density profiles of both condensed and noncondensed atoms change dramatically. In Fig.5, we plot the temporal evolution of the density distributions of condensed atoms as well as noncondensed atoms. In the stable regime, the condensate density oscillates regularly with time and shows clear beating pattern (Fig.5(a)), whereas the density of the noncondensed atoms grows slowly and shows main peaks around $`\theta =\pm \pi `$ and $`0`$, besides some small oscillations (Fig.5(b)). In the unstable regime, the temporal oscillation of the condensate density is irregular (Fig.5(c)), whereas the density of noncondensed atoms grows explosively with the main concentration peaks at $`\theta =\pm \pi /2`$ where the gradient density of the condensed part is maximum (Fig.5(d)). Moreover, our numerical explorations show that the $`\mathrm{cos}^2\theta `$ mode (Fig5.(b)) dominates the density distribution of the noncondensed atoms as the interaction strength is less than 1.8. Thereafter, the $`\mathrm{sin}^2\theta `$ mode grows while $`\mathrm{cos}^2\theta `$ mode decays, and finally $`\mathrm{sin}^2\theta `$ mode become dominating in the density distribution of noncondensed atoms above the transition point (Fig5.(d)). Since the density distribution can be measured in experiment, this effect can be used to identify the transition to instability.

Finally, although the above discussions have been focused on a periodic state of anti-resonance, the transition to instability due to strong interactions also follows a similar path for a dynamically localized state detail . The only difference is that we start out with a quasiperiodic rather than periodic motion in the absence of interaction. This means that it will generally be easier to induce instability but still requires a finite strength of interaction. In Fig.6, we show the nonlinear effect on a dynamically localized state at $`K=5`$ and $`T=1`$. For weak interactions ($`g=1`$) the motion is quasiperiodic with slow growth in the number of noncondensed atoms. Strong interaction ($`g=5`$) destroys the quasiperiodic motion and leads to diffusive growth of energy, accompanied with exponential growth of noncondensed atoms that clearly indicates the instability of the BEC. Notice that the rate of growth in energy is much slower than the classical diffusion rate, which means that chaos brought back by interaction in this quantum system is still much weaker than pure classical chaos.

We acknowledge the support from the NSF, and the R. A. Welch foundation, MGR also acknowledges supports from the Sid W. Richardson Foundation.