text:
# 1 Introduction

## 1 Introduction

Deep inelastic neutral current scattering $`epeX`$ at very high squared momentum transfer $`Q^2`$ allows one to study the structure of $`eq`$ interactions at short distances and to search for new phenomena beyond the Standard Model (SM). The concept of four-fermion contact interactions (CI) provides a convenient method to investigate the interference of any new particle field associated to large scales with the $`\gamma `$ and $`Z`$ fields of the Standard Model. This paper considers conventional contact interactions, such as general models of compositeness and the exchange of heavy leptoquarks and supersymmetric quarks, as well as low scale quantum gravity effects, which may be mediated via gravitons coupling to Standard Model particles and propagating into large extra spatial dimensions.

The present analysis is a continuation of previous studies  based on $`e^+p`$ data . New H1 cross section data on $`e^{}p`$ scattering  and on $`e^+p`$ scattering data  lead to a substantially higher sensitivity to new physics phenomena. Similar studies of $`eq`$ contact interactions have been performed by other experiments at Hera , Lep  and Tevatron , providing results comparable to those presented here.

## 2 Data and analysis method

The data have been collected with the H1 detector at Hera and correspond altogether to an integrated luminosity of $`_{int}=117.2\mathrm{pb}^1`$. They consist of three data sets of different lepton charge and centre of mass energy $`\sqrt{s}`$:

| reaction | $`_{int}[\mathrm{pb}^1]`$ | $`\sqrt{s}[\mathrm{GeV}]`$ | run period | ref. |
| --- | --- | --- | --- | --- |
|  | | | | |
| $`e^+pe^+X`$ | $`35.6\pm 0.53`$ | 301 | $`19941997`$ |  |
| $`e^{}pe^{}X`$ | $`16.4\pm 0.3`$ | 319 | $`19981999`$ |  |
| $`e^+pe^+X`$ | $`65.2\pm 0.95`$ | 319 | $`19992000`$ |  |

The cross section measurements $`\mathrm{d}\sigma /\mathrm{d}Q^2`$ extend over a large range of squared momentum transfers $`200\mathrm{GeV}^2<Q^2<30000\mathrm{GeV}^2`$ for inelasticity $`y<0.9`$. Details can be found in the quoted references. In this paper a contact interaction analysis is presented for each of the two recent data sets and the combined data including the previous measurements.

The phenomenological models under study and the analysis method are described in more detail in ref. . The analysis investigates the measured cross sections $`\mathrm{d}\sigma /\mathrm{d}Q^2`$ and performs quantitative tests of the SM or CI models, applying a minimisation of the $`\chi ^2`$ function

$$\chi ^2=\underset{i}{}\left(\frac{\widehat{\sigma }_i^{\mathrm{exp}}\widehat{\sigma }_i^{\mathrm{th}}(1_k\mathrm{\Delta }_{ik}(\epsilon _k))}{\mathrm{\Delta }\widehat{\sigma }_i}\right)^2+\underset{k}{}\epsilon _k^2.$$
(1)
Here $`\widehat{\sigma }_i^{\mathrm{exp}}`$ and $`\widehat{\sigma }_i^{\mathrm{th}}`$ are the experimental and theoretical cross sections for the measurement point $`i`$ and $`\mathrm{\Delta }\widehat{\sigma }_i`$ is the corresponding error including statistical and uncorrelated systematic errors added in quadrature. The functions $`\mathrm{\Delta }_{ik}(\epsilon _k)`$ describe correlated systematic errors for point $`i`$ associated to a source $`k`$. They depend on the fit parameters $`\epsilon _k`$, which are effectively pulls caused by systematics. In general the influence of the correlated systematic errors is small. The following dominant sources of correlations are included: an overall normalisation uncertainty of $`1.51.8\%`$ depending on the run period, experimental uncertainties on the scattered lepton energy and angle and an uncertainty on the strong coupling $`\alpha _s(M_Z)=0.118\pm 0.003`$ entering the SM prediction. In a combined $`\chi ^2`$ analysis using eq. (1) all data sets are treated as independent samples with individual normalisation and measurement uncertainties. The influence on the fit results of correlations between data sets is negligible.

The cross sections $`\mathrm{d}\sigma (e^{}pe^{}X)/\mathrm{d}Q^2`$  and $`\mathrm{d}\sigma (e^+pe^+X)/\mathrm{d}Q^2`$  from the new data are shown in figure 1.111 Corresponding figures for the $`e^+p`$ cross sections at $`\sqrt{s}=301\mathrm{GeV}`$  are shown in .  They are very well described by the Standard Model expectations over the full $`Q^2`$ range, over which they vary by six orders of magnitude. Cross section calculations in the Standard Model are performed in the DIS scheme in next-to-leading order (NLO) QCD using the corresponding CTEQ5D parton parametrisation . These parton distributions were obtained including the small amount of $`2.7\mathrm{pb}^1`$ of 1994 H1 data with $`Q^2<5000\mathrm{GeV}^2`$ only and are thus essentially uncorrelated with the present data.

Fits to the SM prediction yield $`\chi ^2=10.4`$ for $`17`$ degrees of freedom (dof) with a fitted normalisation of $`0.998`$ for the $`e^{}p`$ data and $`\chi ^2/\mathrm{dof}=13.9/17`$ with a normalisation of $`1.006`$ for the $`e^+p`$ data, using statistical errors from the experiment and including all systematics as in eq. (1). Applying other parton density functions like CTEQ6 , MRST  or GRV  leads to slight changes of the normalisation within errors, but yields equally good agreement with the data.

The cross section measurements do not show significant deviations from the Standard Model and will be used to set limits on couplings from processes mediated through contact interactions. Since the concept of contact interactions is an effective theory, theoretical expectations cannot be formulated consistently in NLO. A sensible approach is to reweight at each $`Q^2`$ value the SM NLO cross section by $`\widehat{\sigma }_i^{LO}(\mathrm{SM}+\mathrm{CI})/\widehat{\sigma }_i^{LO}(\mathrm{SM})`$, the ratio of leading order (LO) cross sections with and without inclusion of contact interactions. The $`\chi ^2`$ function of eq. (1) can be applied to evaluate the sensitivity of the data to a certain CI scenario and to determine its parameters. The corresponding limits at 95% confidence level (CL) are derived by using a frequentist approach as described in the Appendix. Systematic uncertainties are included in this procedure. This determination of limits is different from the method used in our previous publication .

## 3 Contact interaction formalism

The most general chiral invariant Lagrangian for neutral current vector-like four-fermion contact interactions can be written in the form

$`_V`$ $`=`$ $`{\displaystyle \underset{q}{}}{\displaystyle \underset{a,b=L,R}{}}\eta _{ab}^q(\overline{e}_a\gamma _\mu e_a)(\overline{q}_b\gamma ^\mu q_b).`$ (2)
For lepton $`e`$ and each up-type and down-type quark flavour $`q`$ with the corresponding currents $`e_a`$ and $`q_b`$ there are four coupling coefficients $`\eta _{ab}^q=ϵ_{ab}(g/\mathrm{\Lambda }_{ab}^q)^2`$, where $`a`$ and $`b`$ indicate the $`L`$ (left-handed) and $`R`$ (right-handed) fermion helicities, $`g`$ is the overall coupling strength, $`\mathrm{\Lambda }_{ab}^q`$ is a scale parameter and $`ϵ_{ab}`$ determines the interference sign with respect to the Standard Model currents. Any particular model, such as compositeness or the exchange of leptoquarks or supersymmetric quarks, can be constructed by an appropriate choice of the couplings $`\eta _{ab}^q`$. The phenomenological models of interest and their analytical treatment are discussed in more detail in .

## 4 Compositeness scales

In general models allowing for fermion compositeness or substructure it is convenient to choose a coupling strength of $`g^2=4\pi `$ and to assume a universal scale $`\mathrm{\Lambda }`$ for all quarks. The contact interaction coefficients then read

$$\eta _{ab}^q=ϵ_{ab}\frac{4\pi }{\mathrm{\Lambda }^2}.$$

Various scenarios of chiral structures, e.g. pure L(eft) and R(ight) or V(ector) and A(xial vector) couplings, are defined by setting particular chiral contributions to values of $`ϵ_{ab}=\pm 1`$ and setting all other combinations to zero.

The sensitivity to CI models is tested by determining the quantities $`ϵ/\mathrm{\Lambda }^2`$ in a $`\chi ^2`$ fit using the experimental statistical and systematic errors and leaving the sign of interference free. The results from a combined fit to all data sets are shown in figure 2. The data tend to prefer negative interference, but are consistent with the Standard Model within two standard deviations in all scenarios

Lower limits on compositeness scale parameters $`\mathrm{\Lambda }^\pm `$, associated to positive or negative interference, are derived from a frequentist approach and are summarised in table 1. Despite the lower integrated luminosity, the $`e^{}p`$ data exhibit for some models comparable ($`LL`$ and $`RR`$ couplings) or even higher ($`AA^+`$ coupling) sensitivity than the $`e^+p`$ data. The two lepton charges complement each other and a combined analysis of all $`e^\pm p`$ data sets yields substantial improvements in most scenarios compared with the previous publication . The results are presented in figure 3. The lower limits on compositeness scale parameters vary between $`1.6\mathrm{TeV}`$ and $`5.5\mathrm{TeV}`$ depending on the chiral structure. Choosing different parton distributions  changes the quoted limits typically by a few per cent, at most by $`15\%`$. The most restrictive bounds are observed for the $`VV`$ model, where all chiral components contribute with the same sign. As an illustration, examples of fits to the $`e^{}p`$ and $`e^+p`$ cross sections are shown in figure 4.

## 5 Leptoquarks

Leptoquarks couple to lepton–quark pairs and appear in extensions of the Standard Model which try to connect the lepton and quark sectors. They are colour triplet scalar or vector bosons, carrying lepton ($`L`$) and baryon ($`B`$) number and a fermion number $`F=L+3B`$. Since $`F=2`$ for $`e^{}q`$ and $`F=0`$ for $`e^+q`$ states, one expects different sensitivities to particular leptoquark types from electron and positron scattering. For high enough mass scales the leptoquark mass $`M_{LQ}`$ and its coupling $`\lambda `$ are related to the contact interaction coefficients via

$$\eta _{ab}^q=ϵ_{ab}^q\frac{\lambda ^2}{M_{LQ}^2}.$$

Within the model of  the production and decay modes are fixed and the relative contributions $`ϵ_{ab}^q`$ have been calculated . The notation, quantum numbers and couplings of the various leptoquarks are given in table 2.

The analyses of the cross section measurements do not show an indication of a leptoquark signal. The results of fits for each type of leptoquark are interpreted in terms of limits on the ratio $`M_{LQ}/\lambda `$ and are summarised in table 2. In some cases, e.g. $`S_1^L`$ and $`V_1^L`$, the $`e^{}p`$ data give more restrictive bounds than $`e^+p`$ scattering despite the lower integrated luminosity. This sensitivity is illustrated in figure 4, which shows possible contributions of the leptoquarks $`S_1^L`$ and $`V_1^L`$ to the $`e^{}p`$ and $`e^+p`$ cross sections. The two leptoquarks differ by their spin and couple with the same chiral structure but different strength and sign to $`u`$ and $`d`$ quarks. This emphasises the complementary role and importance of both electron and positron beams. The combined analyses of all $`e^\pm p`$ data further constrain the search for leptoquarks, reaching exclusion values of up to $`M_{LQ}/\lambda =1.4\mathrm{TeV}`$. Note that upper bounds on the coupling strength $`\lambda `$ can only be set for leptoquark masses exceeding the accessible centre of mass energy of Hera.

## 6 Squarks in $`𝑹_𝒑`$ violating supersymmetry

In the most general formulation of supersymmetry there exist operators which couple a lepton-quark pair to a squark, the scalar superpartner of a quark. Such couplings violate $`R`$-parity (via lepton number violation), defined as $`R_p=(1)^{3B+L+2S}`$ with S being the spin. Thus $`R_p=+1`$ for SM particles and $`R_p=1`$ for superpartners. This interaction allows single squarks to be produced or exchanged in deep inelastic scattering  via

$`e^+d_R`$ $``$ $`\stackrel{~}{u}_L,\stackrel{~}{c}_L,\stackrel{~}{t}_L\mathrm{coupling}\lambda _{1j1}^{},`$ (3)
$`e^+\overline{u}_L`$ $``$ $`\overline{\stackrel{~}{d}}_R,\overline{\stackrel{~}{s}}_R,\overline{\stackrel{~}{b}}_R\mathrm{coupling}\lambda _{11k}^{}.`$ (4)
The subscripts $`ijk`$ of the coupling $`\lambda _{ijk}^{}`$ describe the generation indices of the left-handed leptons, the left-handed quarks and the right-handed down-type quarks of the superfields, respectively.

The $`e^+q`$ coupling of the left-handed up-type squark of reaction (3) is the same as that of the scalar leptoquark $`\stackrel{~}{S}_{1/2}^L`$, and the coupling of the right-handed down-type squark of reaction (4) is the same as that of the scalar leptoquark $`S_0^L`$. Therefore the formalism and results of the leptoquark analysis can be directly applied. Limits on the ratio $`M_{\stackrel{~}{q}}/\lambda ^{}`$ for $`R_p`$ violating squarks, assuming branching ratios $`_{\stackrel{~}{q}eq}=1`$, are given in table 3. Note that the squark generations cannot be distinguished in this analysis.

## 7 Large extra dimensions

It has been suggested that the gravitational scale $`M_S`$ in $`4+n`$ dimensional string theory may be as low as the electroweak scale of order TeV . The relation to the Planck scale $`M_P10^{19}\mathrm{GeV}`$ and the size $`R`$ of the $`n`$ compactified extra dimensions is given by $`M_P^2R^nM_S^{2+n}`$. In some models with large extra dimensions the SM particles reside on a four-dimensional brane, while the spin 2 graviton propagates into the extra spatial dimensions and appears in the four-dimensional world as a tower of massive Kaluza-Klein states with a level spacing $`\mathrm{\Delta }m=1/R`$. The gravitons couple to the SM particles via the energy-momentum tensor with a tiny strength given by the inverse Planck scale. However, the summation over the enormous number of Kaluza-Klein states up to the ultraviolet cut-off scale, taken as $`M_S`$, leads to an effective contact-type interaction  with a coupling

$$\eta _G=\frac{\lambda }{M_S^4}.$$

The coefficient $`\lambda `$ depends on details of the theory and is expected to be of order unity. However, by convention, one also allows for a negative coupling and thus sets $`\lambda =\pm 1`$. The scale dependence of gravitational effects is quite different from that of conventional contact interactions discussed in the previous sections.

The cross section formulae for virtual graviton exchange in DIS have been calculated in  using the conventions and formalism of . The interference of the graviton with the photon and $`Z`$ fields has opposite sign for electron and positron scattering, as illustrated in figure 4. Lower limits on the scale parameter $`M_S`$, derived from fits to the $`\mathrm{d}\sigma /\mathrm{d}Q^2`$ distributions including gravitational contributions, are summarised in table 4. There is similar sensitivity to the effective gravitational scale for positive and negative interference, resulting in lower limits of $`M_S>0.82\mathrm{TeV}`$ for $`\lambda =+1`$ and $`M_S>0.78\mathrm{TeV}`$ for $`\lambda =1`$.

In other scenarios the gauge bosons $`\gamma `$ and $`Z`$ are also assumed to propagate into extra dimensions. This leads to analogous Kaluza-Klein states which couple to matter with electroweak strength and interfere with the ordinary gauge boson fields. In a specific model of $`4+1=5`$ dimensions  with compactification radius $`R=1/M_C`$, the sums over the exchange of Kaluza-Klein gauge bosons essentially modify the photon propagator $`1/Q^21/Q^2+\pi ^2/(3M_C^2)`$ and the $`Z`$ propagator $`1/(Q^2+M_Z^2)1/(Q^2+M_Z^2)+\pi ^2/(3M_C^2)`$. A fit to the data yields $`M_C>1.0\mathrm{TeV}`$ as a lower limit at 95% CL for the compactification scale.

## 8 Form factors

A fermion substructure can also be formulated by assigning a finite size to the electroweak charge distributions. It is convenient to introduce electron and quark form factors $`f(Q^2)`$ which reduce the Standard Model cross section at high momentum transfer :

$`f(Q^2)`$ $`=`$ $`1{\displaystyle \frac{r^2}{6}}Q^2,`$ (5)
$`{\displaystyle \frac{\mathrm{d}\sigma }{\mathrm{d}Q^2}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{d}\sigma ^{\mathrm{SM}}}{\mathrm{d}Q^2}}f_e^2(Q^2)f_q^2(Q^2).`$ (6)
Fits to the data yield upper limits on the particle size $`R=\sqrt{r^2}`$, taken as the root of the mean squared radius of the electroweak charge distribution. Assuming a point-like electron, i.e. setting $`f_e1`$, the radius of the light $`u`$ and $`d`$ quarks can be constrained to $`R_q<1.010^{18}\mathrm{m}`$ at 95% CL. If both electrons and quarks are assumed to have common form factors one obtains limits on fermion sizes of $`R_{e,q}<0.710^{18}\mathrm{m}`$. The results of the analysis are given in table 5.

## 9 Summary

Neutral current deep inelastic $`e^{}p`$ and $`e^+p`$ scattering cross section measurements are analysed to search for new phenomena mediated through $`(\overline{e}e)(\overline{q}q)`$ contact interactions. The data are well described by the Standard Model expectations. The use of electrons and positrons provides complementary information and a combined analysis based on an integrated luminosity of $`117.2\mathrm{pb}^1`$ yields improved limits on scales of new physics. The present analysis supersedes previous results .

Lower limits at 95% CL on $`eq`$ compositeness scale parameters $`\mathrm{\Lambda }^\pm `$ are derived within a model independent analysis. They range between $`1.6\mathrm{TeV}`$ and $`5.5\mathrm{TeV}`$ depending on the chiral structure. A study of virtual leptoquark exchange yields lower limits on the ratio $`M_{LQ}/\lambda `$ between $`0.3\mathrm{TeV}`$ and $`1.4\mathrm{TeV}`$. Squarks in $`R`$-parity violating supersymmetry with masses satisfying $`M_{\stackrel{~}{u}}/\lambda _{1j1}^{}<0.43\mathrm{TeV}`$ and $`M_{\stackrel{~}{d}}/\lambda _{11k}^{}<0.71\mathrm{TeV}`$ can be excluded. Possible effects of low scale quantum gravity with gravitons propagating into extra spatial dimensions are searched for. Lower limits on the effective Planck scale $`M_S`$ of $`0.780.82\mathrm{TeV}`$ are found. Allowing for Kaluza-Klein states of the SM gauge bosons results in a lower limit on the extra dimension compactification scale $`M_C>1.0\mathrm{TeV}`$. A form factor approach yields an upper limit on the size of light $`u`$ and $`d`$ quarks of $`R_q<1.010^{18}\mathrm{m}`$ assuming point-like electrons.

### Acknowledgements

We are grateful to the Hera machine group whose outstanding efforts have made this experiment possible. We thank the engineers and technicians for their work in constructing and now maintaining the H1 detector, our funding agencies for financial support, the Desy technical staff for continual assistance, and the Desy directorate for support and the hospitality extended to the non–Desy members of the collaboration.

## Appendix Setting limits

For the present analysis several methods have been studied to calculate limits and confidence levels and the final results quoted are based on the frequentist approach (see for example ). The $`\chi ^2`$ function of eq. (1) is used as the quality measure of agreement between data and contact interaction models. It allows an easy implementation of systematic uncertainties. The other methods investigated differ in the definitions of the statistical error entering the total error $`\mathrm{\Delta }\widehat{\sigma }_i`$ in the $`\chi ^2`$ expression.

Taking the statistical errors of the experiment, one can test the compatibility of the data with a certain model hypothesis and determine the corresponding parameter, for example $`1/\mathrm{\Lambda }^2`$. Each data point contributes with a fixed weight $`\mathrm{\Delta }\widehat{\sigma }_i^2`$ to the $`\chi ^2`$ function, independent of the model parameter. However, downward fluctuations with respect to the model predictions get a larger weight (smaller statistical error) than upward fluctuations (larger statistical error). This property may enforce asymmetric situations and errors, as is observed in the present data, see figure 2. If the fitted parameter is found to be not significantly different from zero, the result can be converted into limits at a given confidence level (CL). A convenient way is to take the values corresponding to a change of $`\chi ^2(1/\mathrm{\Lambda }^{\pm 2})\chi _{SM}^2=\mathrm{\Delta }\chi _{CL}^2`$, for instance $`\mathrm{\Delta }\chi _{95\%}^2=3.84`$ for a 95% CL limit. This method was used in the previous publication  and is illustrated for several compositeness models by the $`\chi ^2`$ distributions in figure 5. This simple and robust definition is certainly meaningful for parabolic curves, but the probabilistic interpretation becomes ambiguous if there are secondary minima or other structures.

An alternative possibility is to assume the validity of a certain CI model and to calculate the probability to observe the measured value for a given model parameter. Here the statistical error entering the $`\chi ^2`$ function is taken from the model prediction. This approach is completely equivalent to using the log-likelihood function, and it has been verified that both methods lead to the same results. Each data point gets a varying weight depending on the model parameter, but independent of fluctuations in the data. For the present data this leads to more symmetric situations and an unbiased evaluation. In general, the resulting $`\chi ^2`$ curves show a stronger sensitivity to details of the CI model, e.g. interference patterns, as can be seen in figure 5. The widths of the distributions are often, but not always, wider than in the previous case. The problem to extract limits from these distributions persists.

The final limits presented here have been determined by applying a frequentist approach. Starting from a specific model with a scale parameter $`\mathrm{\Lambda }_{true}`$ the cross section $`\mathrm{d}\sigma /\mathrm{d}Q^2`$ is calculated and then smeared according to the statistical error given by the predicted number of events. Distortions due to all uncorrelated and correlated systematic uncertainties (except for the parton distributions) are included, assuming Gaussian behaviour of the errors. The Monte Carlo experiment is then analysed in the same way as the data, i.e. with statistical errors from the prediction and including all systematics, resulting in a fitted value $`\mathrm{\Lambda }_{fit}`$. This procedure is repeated numerous times and the fit results are recorded in a probability distribution. The scale parameter is then varied and the 95% confidence lower limit $`\mathrm{\Lambda }^+`$ ($`\mathrm{\Lambda }^{}`$) is defined as that value $`\mathrm{\Lambda }_{true}`$ where 95% of the Monte Carlo experiments produce values of $`\mathrm{\Lambda }_{fit}`$ which are larger (smaller) than the parameter $`\mathrm{\Lambda }`$ actually obtained from the data. Examples of CL distributions are shown in figure 5. For the present analysis the frequentist approach provides in general slightly weaker and more symmetric limits than the statistical method used previously .


format: text

text:
# Rock-scissors-paper game on regular small-world networks

## 1 Introduction

The rock-scissors-paper like cyclic dominance among three states (modes, species, strategies, opinions) are widely studied in different spatial systems. For example, the Rayleigh-Bernard convection in fluid layers rotating around a vertical axis exhibits the Küpper-Lortz instability  resulting in a cyclic change of the three possible directions of parallel convection rolls . Such a situation can appear in biological (ecological) systems too . Very recently, Kerr et al.  have justified experimentally that the cyclic dominance between the toxic, sensitive and resistant microbes maintains their coexistence as predicted previously by several theoretical works . The emergence of cyclic invasions has also been observed for some three-strategy evolutionary games where this phenomenon promotes the cooperation among selfish individuals .

The above three-state systems exhibit some general features. Namely, spiral formation (or rotating three-edge vortices and antivortices) can occur on the two-dimensional backgrounds  as it is observed for the Belousov-Zhabotinskii reaction as well as for the numerical investigation of excitable activator-inhibitor media . This self-organizing structure can provide a stability against some external invaders and thereby it plays crucial role in ecological systems . Furthermore the cyclic dominance yields a paradoxical response to the external support of one of the species  and global oscillation can occur when varying either the model  or the structural parameters  if long range interactions are allowed.

In this work the spatial rock-scissors-paper game will be extended for such regular networks where each site has four neighbours. We discuss how the quenched and annealed randomness of the network affects the spatiotemporal distribution of species. Such a comparison has already been performed for a rumor propagation model . Now, the randomness is introduced in small-world manner  leaving the degree of sites unchanged. This restriction does not influence the small-world feature of the network and it simplifies the application of some theoretical approximation at least for the case of annealed randomness. However, the main conclusion remains valid for both types of randomness. The Monte Carlo (MC) simulations indicate that these systems undergo a transition (Hopf bifurcation) from a stationary state (with fixed average concentrations) to a global oscillation. The amplitude of oscillation increases with the measure of randomness and the increasing oscillation terminates at one of the absorbing (homogeneous) states above a second threshold value for annealed randomness. We show that neither the mean-field nor the pair approximations can explain these transitions. The failure of these descriptions has inspired us to use the more sophisticated dynamical cluster techniques where one determines all the configuration probabilities on a $`k`$-site cluster. At the levels of $`k=1`$ and 2 this technique is equivalent to the mentioned mean-field and pair approximations. The essence of this simple method is described in Refs. .

## 2 The model

We consider a very simple model where the sites of a regular graph are occupied by one of the three species ($`s_i=1,2,3`$) that dominate cyclically each other. The evolutionary process is governed by the sequence of elementary invasions along the randomly chosen edges of the graph. Namely, first we choose a site and one of its linked (neighbouring) sites at random. If these two sites are occupied by different species then the predator occupies the prey’s site, i.e., the $`(1,2)`$ and $`(2,1)`$ pairs transform into $`(1,1)`$, $`(2,3)`$ and $`(3,2)`$ into $`(2,2)`$, and $`(1,3)`$ and $`(3,1)`$ into $`(3,3)`$. Starting from a random initial state the above process is repeated until the system reaches the stationary state or limit cycle we study. This system was already analysed systematically by several authors if the sites form a $`d`$-dimensional lattice .

Two types of random structures will be contrasted with each other. In the first case the structure is quenched after the generation of edges. Figure 1 illustrates an example whose creation is similar to those suggested by Watts and Strogatz . Notice, that the present algorithm (explained in the caption of Fig. 1) conserves the degree of sites, i.e., the number of neighbours remains fixed ($`z=4`$) for each site. In the second case the “neighbourhood” is not fixed in time, that is a randomly chosen new site can be replaced for any standard neighbours during the elementary invasions explained above. These types of random networks can characterize the interaction among individuals in social systems . In both cases the random links connect two sites being arbitrary distance from each other in the original structure.

The measure of quenched randomness is characterized by the $`Q`$ portion of random links substituted for the nearest neighbour bonds. If $`Q=0`$ this structure reproduces the square lattice and for $`Q=1`$ it is equivalent to a random regular graph  where the typical local structure is analogous to a Cayley tree for large number of sites $`N`$. Previous investigations indicated that some phenomena on the random regular graphs can be well described analytically if the background is assumed to be a Bethe lattice in the limit $`N\mathrm{}`$ . In other words, the structure with $`Q=1`$ can be considered as a possible realization of the Bethe lattice in the simulations for large $`N`$.

Unfortunately, the present random regular structures with $`0<Q<1`$ are not yet investigated rigorously, although many other classes of networks are well studied . We think that the constraint of regularity leaves the relevant features unchanged and the present structure remains similar to those introduced by Watts and Strogatz  on a square lattice. When increasing $`Q`$ the present structure exhibits a continuous transition from the square lattice to the random regular graph.

Besides the above quenched randomness we will investigate the consequences of the annealed (temporal) randomness in the structure. In this case the standard links are defined by the bonds between the nearest neighbours sites forming a square lattice. Occasionally the standard link is replaced by a random one with a probability $`P`$ characterizing the strength of annealed randomness. Evidently, for $`P=0`$ the structure is equivalent to the square lattice. On the other hand, in the limit $`P1`$ this system satisfies the conditions of mean-field approaches.

## 3 Simulations

The MC simulations are performed on a network consisting of $`N=L\times L`$ sites where the linear size of the square lattice ($`L`$) is varied from 400 to 3200. The regular small-world networks are constructed from a square lattice repeating the rewiring steps $`2N\mathrm{ln}(1Q)`$ times as explained in the caption of Fig. 1. (The logarithmic correction takes into account that the substitution of a random link for another random one does not increase the quenched randomness.)

The above model has two parameters characteristic to the quenched ($`Q`$) and annealed ($`P`$) randomness of the background. The present analysis is restricted to those cases when one of them is chosen to be zero.

For small sizes ($`L<10`$) this system evolves into one of the three absorbing states where all the sites are occupied by the same species. For sufficiently large system sizes, however, the three species can coexist and after some transient time the state can be well described by the current concentration of the three species \[$`c_1(t)+c_2(t)+c_3(t)=1`$\]. In order to observe these states we have to choose such a large $`L`$ where the amplitude of fluctuations becomes significantly less than the minimum value of concentrations. The above choices satisfy this condition.

On the square lattice ($`Q=P=0`$) this system develops into a stationary state where all three species are present with the same average concentration, i.e., $`c_1=c_2=c_3=1/3`$ corresponding to the central point in the ternary phase diagram as plotted in Fig. 2. In this case the three species alternate cyclically at each site and the short range interactions are not able to synchronize these local oscillations.

The corresponding self-organizing spatiotemporal patterns are well investigated previously by several authors . In these patterns the rotating vortices (spirals) and antivortices are not recognizable because of the interfacial roughening. The absence of smooth interfaces (surface tension) is caused by the fact that the elementary invasions between two neighbouring sites are not affected by their neighbourhood .

Global oscillation (synchronization) occurs when the frequency of random (long range) links exceeds a threshold value dependent on whether quenched or annealed randomness is considered. To illustrate the time-dependence of the species distributions during a period of this global oscillation six consecutive snapshots are plotted in Fig. 3.

The amplitude of oscillation increases with the measure of randomness in both cases. If the annealed randomness exceeds a second threshold value then the trajectories approach the edges of the triangle and sooner or later the evolution is terminated at one of the homogeneous (absorbing) states (where the system stays for ever). It is emphasized, however, that the homogeneous states are not stable against the invasion of their predators. For example, the offsprings of a single species 3 will occupy the whole territory of species 1 in the absence of species 2.

In the ternary phase diagram the shape of the limit cycles reflects the cyclic symmetry between the species and its extension is described by the area $`A`$ compared to its maximum value. The average value of this quantity can be easily determined by numerical integration after a suitable relaxation time for either the MC simulations or the dynamical cluster techniques. Evidently, $`A`$ vanishes if the system tends toward the central fixed point, and it goes to 1 when the trajectories approach the edges of triangle (see Fig. 2).

Systematic MC simulations are carried out to determine the average value of $`A`$ on the regular small-world structure for different $`Q`$ values. The results in Fig. 4 refer to a transition from a fixed point to the limit cycle. For weak randomness \[$`Q<Q_1=0.067(1)`$\] the system always tends to the central fixed point. Conversely, oscillating behaviour occurs and the area $`A`$ (as well as the amplitude) of the limit cycle increases monotonously with $`Q`$ and tends to the value $`A(Q=1)=0.980(1)`$. The limit $`Q1`$ is investigated on a random regular graph created by using a different algorithm . In the vicinity of the transition point $`A`$ vanishes linearly with $`QQ_1`$ in agreement with the prediction of Hopf bifurcation describing the emergence of a limit cycle in a mean-field type system if the model parameters are varied. Further rigorous investigations are required to quantify the variations in the spatiotemporal patterns when approaching the transition point.

As it is already mentioned the system size should be large enough to avoid the finite size effects. The sharpe transition to the global oscillation is wiped out by fluctuations on smaller systems as discussed previously by Kuperman and Abramson  who considered a three-state (cyclic) epidemiological model on small world networks suggested by Watts and Strogatz . The finite size effects are more dangerous when the value of area $`A`$ approaches 1. In this case the evolution may be easily traped by one of the absorbing states. This difficulty is avoided by choosing sufficiently large system size. As a result, the plotted $`A(Q)`$ function may be considered as the infinite size limit.

The above analysis was repeated for the annealed randomness when varying $`P`$ for $`Q=0`$. The results of MC simulations are summarized in Fig. 5. The oscillating behaviour can be observed for $`P_1<P<P_2`$ where $`P_1=0.020(1)`$ and $`P_2=0.170(1)`$. If $`P>P_2`$ then the increasing spiral trajectory terminates in one of the absorbing states as demonstrated in Fig. 2. In the vicinity of the first threshold value $`A(PP_1)`$ in agreement with the expectation. On the contrary, $`A`$ approaches 1 very smoothly when $`P`$ goes to $`P_2`$. Surprisingly, our MC data can be well approximated by a power law behaviour \[$`1A(P_2P)^\gamma `$\] as indicated in the inset of Fig. 5. The numerical fit yields $`\gamma =3.3(3)`$.

It is worth mentioning that the emergence of global oscillation has already been observed in the above mentioned epidemiological model when varying the quenched randomness without the consraints of regularity . The similar behavior refers to the robustness of this type of transition.

## 4 Extended mean-field analysis

The predictions of the traditional mean-field analysis are well discussed in the textbook by Hofbauer and Sigmund . According to this approach four stationary states exist, namely, the above mentioned three absorbing states (e.g. $`c_1=1`$ and $`c_2=c_3=0`$) and the well-mixed symmetric state where the three species are present with the same concentration (1/3). Besides these stationary solutions the mean-field analysis shows the existence of set of oscillating states whose closed trajectories draw concentric orbits around the centrum in the ternary phase diagram. Along these orbits the product $`c_1c_2c_3`$ is conserved. In agreement with the expectation the MC simulations reproduce this behaviour for $`P=1`$ as shown in Fig. 2.

The application of the pair approximation is inspired by its success when investigating an evolutionary game with three (cyclically dominated) strategies on a random regular graph . This model exhibits both transitions mentioned above when varying the parameter(s) of payoff matrix. It is underlined that here the adoption of the neighbouring strategies depends on the neighbourhood.

In the pair approximations one determines the $`p_2(s_1,s_2)`$ probability of finding $`(s_1,s_2)`$ configuration on two nearest neighbour (or linked) sites. In this case equations of motion are derived for these quantities taking into account the contribution of all the elementary invasion processes (details are given in the textbook by Marro and Dickman ). The numerical integration of the corresponding equations predicts growing spirals approaching the boundaries (and resulting in the maximum value $`A=1`$ in the limit $`t\mathrm{}`$) as indicated in Fig. 2. This prediction is contrary to the results of MC simulations obtained on either the square lattice ($`A=0`$) or the random regular graphs ($`A=0.98`$). We have to emphasize that although the equations of motion involve explicitly the number of neighbours the simple pair approximation can not distinguish the structure of the Bethe and square lattices. In the light of previous experiences  it is expected that the pair approximation can well describe the behaviour observed by MC simulations on the random regular graph because in this structure the average loop size increases with with $`\mathrm{ln}N`$ . Thus, for large $`N`$, the local structure is tree-like and the short range correlations between two sites can be well approximated by a product of the $`p_2(s,s^{})`$ configuration probabilities involved along the single path connecting the two sites. The comparison of the above values of $`A`$ does not confirm this expectation. More precisely, the pair approximation can not account for the effect preventing the unlimited growth of the area of limit cycle. Preliminary results suggest that the limit value of $`A`$ depends on the degree of the random regular graph. In the near future we wish to study this effect by a suitable extension of the dynamical cluster techniques. Henceforth, however, our efforts will be concentrated on the annealed randomness because its investigation can be easily builted into the dynamical cluster technique. At the level of pair approximation the corresponding results predict that the ”spiral pitch” decreases when $`P`$ is increased and vanishes in the mean-field limit ($`P=1`$) as expected.

On the square lattice, as mentioned above, the pair approximation is not capable to describe the appearance of self-organizing patterns maintained by the interfacial motions due to cyclic invasions. This striking shortage can be eliminated by choosing larger and larger clusters on which we determine all the possible configuration probabilities . For example, at the level of four-site approximation we determine the configuration probabilities $`p_4(s_1,s_2,s_3,s_4)`$ on a $`2\times 2`$ cluster assumed to be translation invariant on the square lattice. A recent summary of larger-cluster approximations can be found in Ref. .

This approach takes explicitly into account some topological features of the square lattice. A nearest-neighbor invasion (e.g. $`s_4s_5`$ as demonstrated in Fig. 6) will simultaneously affect all the four four-site configuration probabilities involved. The spatial effect is taken into account more rigorously if the probability of such a nine-site constellation is approximated as

$`p_9(s_1,\mathrm{},s_9)=`$ $`{\displaystyle \frac{p_4(s_1,s_2,s_4,s_5)p_4(s_2,s_3,s_5,s_6)}{p_2(s_2,s_5)p_2(s_4,s_5)}}`$
$`\times {\displaystyle \frac{p_4(s_4,s_5,s_7,s_8)p_4(s_5,s_6,s_8,s_9)}{p_2(s_5,s_6)p_2(s_5,s_8)}}p_1(s_5),`$
where the configuration probabilities satisfy the following compatibility conditions:

$`p_1(s_1)=`$ $`{\displaystyle \underset{s_2}{}}p_2(s_1,s_2)={\displaystyle \underset{s_2}{}}p_2(s_2,s_1),`$
$`p_2(s_1,s_2)`$ $`=`$ $`{\displaystyle \underset{s_3,s_4}{}}p_4(s_1,s_2,s_3,s_4)={\displaystyle \underset{s_3,s_4}{}}p_4(s_3,s_4,s_1,s_2)`$
$`=`$ $`{\displaystyle \underset{s_3,s_4}{}}p_4(s_1,s_3,s_2,s_4)={\displaystyle \underset{s_3,s_4}{}}p_4(s_3,s_1,s_4,s_2).`$
The time derivative of $`p_4(s_1,s_2,s_3,s_4)`$ can be described by a master equation that summarizes the contribution of all the possible elementary invasions. Namely,

$`{\displaystyle \frac{d}{dt}}p_4(s_1,s_2,s_3,s_4)=`$
$`P{\displaystyle \underset{s_x}{}}p_4(s_1,s_2,s_3,s_4)p_1(s_x)[\mathrm{\Gamma }(s_x,s_1)+\mathrm{\Gamma }(s_x,s_2)+\mathrm{\Gamma }(s_x,s_3)+\mathrm{\Gamma }(s_x,s_4)]`$
$`+P{\displaystyle \underset{s_x}{}}p_4(s_x,s_2,s_3,s_4)p_1(s_1)\mathrm{\Gamma }(s_1,s_x)`$
$`+P{\displaystyle \underset{s_x}{}}p_4(s_1,s_x,s_3,s_4)p_1(s_2)\mathrm{\Gamma }(s_2,s_x)`$
$`+P{\displaystyle \underset{s_x}{}}p_4(s_1,s_2,s_x,s_4)p_1(s_3)\mathrm{\Gamma }(s_3,s_x)`$
$`+P{\displaystyle \underset{s_x}{}}p_4(s_1,s_2,s_3,s_x)p_1(s_4)\mathrm{\Gamma }(s_4,s_x)`$ (3)
$`{\displaystyle \frac{1P}{4}}{\displaystyle \underset{s_5,s_6,s_7,s_8,s_9}{}}p_9(s_1,s_2,s_5,s_3,s_4,s_6,s_7,s_8,s_9)\mathrm{\Gamma }(s_2,s_4)`$
$`{\displaystyle \frac{1P}{4}}{\displaystyle \underset{s_5,s_6,s_7,s_8,s_9}{}}p_9(s_5,s_1,s_2,s_6,s_3,s_4,s_7,s_8,s_9)\mathrm{\Gamma }(s_1,s_3)`$
$`{\displaystyle \frac{1P}{4}}{\displaystyle \underset{s_5,s_6,s_7,s_8,s_9}{}}p_9(s_5,s_6,s_7,s_1,s_2,s_8,s_3,s_4,s_9)\mathrm{\Gamma }(s_6,s_2)`$
$`{\displaystyle \frac{1P}{4}}{\displaystyle \underset{s_5,s_6,s_7,s_8,s_9}{}}p_9(s_5,s_6,s_7,s_8,s_1,s_2,s_9,s_3,s_4)\mathrm{\Gamma }(s_6,s_1)`$
$`+{\displaystyle \frac{1P}{4}}{\displaystyle \underset{s_x,s_5,s_6,s_7,s_8,s_9}{}}p_9(s_1,s_2,s_5,s_3,s_x,s_6,s_7,s_8,s_9)\delta (s_2,s_4)\mathrm{\Gamma }(s_2,s_x)`$
$`+{\displaystyle \frac{1P}{4}}{\displaystyle \underset{s_x,s_5,s_6,s_7,s_8,s_9}{}}p_9(s_5,s_1,s_2,s_6,s_x,s_4,s_7,s_8,s_9)\delta (s_1,s_3)\mathrm{\Gamma }(s_1,s_x)`$
$`+{\displaystyle \frac{1P}{4}}{\displaystyle \underset{s_x,s_5,s_6,s_7,s_8,s_9}{}}p_9(s_5,s_6,s_7,s_1,s_x,s_8,s_3,s_4,s_9)\delta (s_6,s_2)\mathrm{\Gamma }(s_6,s_x)`$
$`+{\displaystyle \frac{1P}{4}}{\displaystyle \underset{s_x,s_5,s_6,s_7,s_8,s_9}{}}p_9(s_5,s_6,s_7,s_8,s_x,s_2,s_9,s_3,s_4)\delta (s_6,s_1)\mathrm{\Gamma }(s_6,s_x)`$
$`+\mathrm{}`$
where $`\delta (s_x,s_y)`$ denotes the Kronecker’s delta and the constraint of invasion is expressed as

$$\mathrm{\Gamma }(s_x,s_y)=\{\begin{array}{cc}1,\hfill & \text{if }s_y1=s_x\text{mod}\mathrm{\hspace{0.33em}3}\text{ ;}\hfill \\ 0,\hfill & \text{otherwise .}\hfill \end{array}$$
(4)
The terms proportional to $`P`$ describe the contributions coming from the invasions from arbitrary distance while the contributions from one of the four nearest-neighbor sites are proportional to $`(1P)/4`$. Equation (3) involves explicitly only those terms coming from the downward invasions. The derivation of the missing terms is straightforward.

At the level of six-site approximation the probability of a nine-site configuration (as shown in Fig. 6) is expressed by the product of configuration probabilities on $`3\times 2`$ clusters as

$$p_9(s_1,\mathrm{},s_9)=\frac{p_6(s_1,s_2,s_3,s_4,s_5,s_6)p_6(s_4,s_5,s_6,s_7,s_8,s_9)}{p_3(s_4,s_5,s_6)}$$
(5)
where the $`p_3(s_1,s_2,s_3)`$ indicates the configuration probabilities on a $`3\times 1`$ cluster. Evidently, in this case the invasion of the central site will influence some other six-site configuration probabilities that we can handle in a similar way.

In both cases the corresponding master equations are solved by numerical integration and the results are summarized in Fig. 5. In agreement with the expectation, the dynamical cluster techniques (at such a high level) reproduce qualitatively well the results obtained by MC simulations. Namely, both descriptions confirm the stability of the central stationary state, that is the area $`A`$ tends to zero, if $`P<P_1^{(4p)}=P_1^{(6p)}=0.011(1)`$. Above this transition point the present approaches predict the appearance of a limit cycle within a range $`P_1<P<P_2`$. The area $`A`$ increases linearly with $`PP_1^{(4p)}`$ in the close vicinity of the transition point. Similarly, $`A`$ approaches to 1 linearly for both the four- and six-site approximations, although these methods predict different transition points, i.e., $`P_2^{(4p)}=0.097(3)`$ and $`P_2^{(6p)}=0.109(3)`$. These sophisticated techniques have significantly improved the description and the deviations from the MC results reflect the relevant role of topological features.

## 5 Conclusions

In summary, cyclic invasions like the rock-scissors-paper in the three-state systems can maintain a rich variety of spatiotemporal patterns that depend on the quenched and annealed randomness of the background. According to the MC simulations a self-organizing pattern occurs on the square lattice. The classical mean-field and pair approximations are not capable to reproduce this behaviour. It is demonstrated, however, that the extended versions of this approach, called dynamical cluster technique at the level of four- and six-site approximations, can describe the appearance of this self-organizing pattern.

The quenched randomness is generated by a modified rewiring technique that conserves the degree at each site. For weak quenched randomness the above spatiotemporal pattern remains stable. When the measure of quenched randomness exceeds a threshold value this system undergoes a transition from the symmetric stationary state (central fixed point) to a synchronized oscillation (limit cycle). For annealed randomness this model exhibits similar behaviour with a higher sensitivity to the variation of annealed randomness and above a second threshold value the increasing oscillation terminates at one of the homogeneous (absorbing) states. These features are reproduced qualitatively well by the dynamical cluster technique considering the configuration probabilities on four- and six-site clusters. We think that further systematic analyses can clarify how the transitions are affected on those systems where the quenched and annealed randomness occur simultaneously. More significant modifications of this technique are required if we wish to study the cyclic invasions on networks with arbitrary degree distributions.

## 6 Acknowledgement

This work was supported by the Hungarian National Research Fund under Grant Nos. T-33098, F-30499 and Bolyai Grant No. BO/0067/00.

## References


text:
# John Clifford Barton

## John Clifford Barton

John Barton was one of the founders of particle astrophysics. This experimental science, born in the aftermath of the Second World War and fuelled by developments in electronics and computing, seeks answers to fundamental cosmological problems. Barton was one of the greatest experimentalists, working in deep underground locations all over world and inspiring a generation of physicists to follow him.

Born in 1923, he was a wartime student at University College London during its evacuation to Bangor. Clearly, wartime requirements influenced the curriculum and he gained a “Certificate of Proficiency in Radio-Physics”, which was to mark out his future path. He graduated in 1943, receiving the Granville Prize, annually awarded to the best Physics graduate of London University. He did his National Service at Marconi in Chelmsford, Essex, working on military electronics. Here he acquired a thorough grounding in the basics of electronic design and construction techniques. After the war he became a PhD student at Birkbeck College, London University, working with E.P. George on cosmic rays. His first experiments were performed at an altitude of 3,457 metres, high on the Jungfraujoch in Switzerland. An amusing aside in the report thanks the director of the Jungfraujoch railway for the loan of 15 tonnes of coal that had been used for absorber. This use of available materials was to become something of a trademark of Barton’s experiments.

In 1954 he began four years at the nascent University College of the West Indies, Jamaica. As well as his teaching, Barton continued his research, measuring the cosmic muon flux as a function of depth in the sea. For this he built a transistorised cosmic-ray “telescope” using Geiger tubes contained in a pressure vessel, which was lowered to depths of 3,000m. Calibrations were made underground at Norton Hill Colliery in Somerset which were Barton’s first experiments in a mine. To record the data he used his own design of tape recorder, capable of recording eight tracks of digital data across standard quarter-inch audio tape. This was the first ever use of digital recording on a scientific experiment.

In 1958, Barton returned to London and became lecturer at Northern Polytechnic, almost immediately publishing his first paper describing a cosmic-ray detector using photomultipliers, which had just become commercially available. Photomultipliers are vacuum tubes which detect faint flashes of light, and their large area of sensitivity and high gain meant that large, robust particle detectors could be built. They are still found at the heart of innumerable physics experiments and huge numbers are used in medical imaging. Barton rapidly became one of the world experts in photomultipliers and their applications. He also had an almost intuitive feel for the “non-imaging” optics needed to carry light to the photomultipliers, and used to say to his students, “Light doesn’t go down a funnel like water does”, when they came up with ideas that didn’t work.

He began a series of experiments to determine the nature of cosmic rays that could penetrate deep underground. These experiments were performed in the “Holborn Laboratory”, a series of rooms deep in Holborn Underground station. A spare platform at Holborn had been converted to offices during the war. Immediately after the war it was used as a staff hostel, and later many of the rooms were used by physicists for experiments needing a deep location. The laboratory rooms were reached through a service door on one of the Piccadilly Line platforms. They were linked by an extremely narrow corridor, only wide enough for a single person, running along the edge of what had once been the platform. It was a dry and dusty environment and there were occasional problems caused by rodents chewing cables, but it was none the less an extraordinarily convenient site to work. For many measurements Holborn was not deep enough and Barton and his colleagues also ran experiments in Tilmanstone Colliery in Kent and later in the Woodhead Tunnel, a disused railway tunnel under the Pennines.

In the early 1980s Barton started on a series of studies on meteorites. His low background laboratory was ideal for identifying trace radioisotopes produced in the meteorites in space before they hit the earth. This work led on to a search for “superheavy elements”. Theoretical analysis suggested that while nuclides heavier than Uranium were unstable, there would be an “island of stability” around element 114 which would have half-lives long enough to exist in nature. Others had already undertaken searches in a range of samples, particularly meteorites, and some had claimed positive results. Together with a group from Leeds, Barton repeated the experiments and, despite having more sensitive equipment, saw no superheavy elements. Years later, element 114 was made artificially at Darmstadt and was found to have a half-life of 30 seconds, a full 15 orders of magnitude smaller than the original predictions. For Barton this was a vindication of his belief that theoretical predictions must be tested by experiment and that theoreticians are often just plain wrong.

When the Physics Department of what was now the Polytechnic of North London closed in 1984, Barton officially retired, devoting himself to research. He held honorary posts at Birkbeck and at Queen Mary, London University. In 1993 the Holborn Laboratory was closed, following increasing safety concerns in the wake of the King’s Cross fire. Barton transferred his underground laboratory to the Eisenhower Centre, a wartime control centre near Goodge Street, and, when the lease on this expired, to the basement at Queen Mary, not really deep enough but workable. Increasing frailty did not deter him – an ingenious assembly of car jacks enabled him single-handed to move several tonnes of lead shielding, no mean feat in a cramped laboratory packed with chemicals, electronics, computers, domestic appliances such as freezers and all the latest state-of-the-art instruments that he could get his hands on.

He became a member of the team that built the Sudbury Neutrino Observatory (SNO) deep in a Canadian nickel mine. The observatory relies heavily on the use of photomultiplier tubes in a harsh environment that must be as free as possible from radioactive contamination. Barton’s contributions were pivotal and without them SNO’s evidence that neutrinos emitted by the Sun change their “flavour” on their way to the Earth would be much less convincing. In 1988, Neil Spooner, then at Imperial College, London, and Professor Peter Smith of the Rutherford Laboratory were forming a new collaboration to study dark matter and to hunt for the elusive “WIMPS” (weakly interacting massive particles). They knew that they needed a deep site to shield the detectors from cosmic rays and that the Boulby potash mine in North Yorkshire was the deepest mine in Britain. Barton was enthusiastic about this new project and went with Spooner on the site visit to help persuade the mine management to accept an underground laboratory. They were successful and the Boulby Dark Matter Collaboration came into being, operating a range of dark matter detectors in the rocksalt seams, one kilometre underground. While he never visited the site again, Barton continued to make vital contributions to the collaboration. The new surface building there has recently been named the John Barton Building.

Much of Barton’s pioneering work was made with relatively cheap equipment he built himself, using the very latest technology available to him, but always on a shoestring budget. Once, when asked why he never applied for grants from the Science Research Council, he replied that, if one applied for a grant, then one had to write reports on the grant, irrespective of the scientific results, and that then one ended up believing what one had written in the reports. Most of his career was spent at Northern Polytechnic, in a period when polytechnics rarely did fundamental research. Barton managed to, despite the environment. In 1968 he put together a pack of 50 research papers and submitted them for a DSc at London University, because he wanted to show that it was possible to do good science in a polytechnic with a supportive head of department.

John Barton was a shy, private and unassuming man but you knew within the first minute of meeting him that you were in the presence of an exceptionally talented and intelligent person. He was an enthusiastic walker and always took an annual walking holiday, most recently a strenuous traverse of Corsica. He loved the cinema and chose to live in Hampstead because of the proximity to the Everyman Cinema. Barton worked in his laboratory almost daily until his final illness. On becoming housebound, he bought the first television he had ever owned, typically finding even Bang & Olufsen’s superior specification left much to be desired.

*John Clifford Barton, physicist: London 29 September 1923; Senior Lecturer in Physics, University College of the West Indies 1954-58; Lecturer in Physics, Northern Polytechnic (later Polytechnic of North London) 1958-84, Head of Physics Department 1971-84; died London 14 October 2002.*

John McMillan
This obituary first appeared in the daily newspaper “The Independant”,
London, 29th November 2002.
j.e.mcmillan@sheffield.ac.uk
Department of Physics and Astronomy,
The University of Sheffield, Sheffield, S3 7RH, Great Britain.


text:
# Optical Response of Solid CO2 as a Tool for the Determination of the High Pressure Phase

## Abstract

We report first-principles calculations of the frequency dependent linear and second-order optical properties of the two probable extended-solid phases of CO2–V, i.e. $`I\overline{4}2d`$ and $`P2_12_12_1`$. Compared to the parent $`Cmca`$ phase the linear optical susceptibility of both phases is much smaller. We find that $`I\overline{4}2d`$ and $`P2_12_12_1`$ differ substantially in their linear optical response in the higher energy regime. The nonlinear optical responses of the two possible crystal structures differ by roughly a factor of five. Since the differences in the nonlinear optical spectra are pronounced in the low energy regime, i.e. below the band gap of diamond, measurements with the sample inside the diamond anvil cell are feasible. We therefore suggest optical experiments in comparison with our calculated data as a tool for the unambiguous identification of the high pressure phase of CO2.

A few years after the synthesis of the extended-solid phase of CO2 at 40 GPa the crystal structure of this phase is still not known. This phase was found to be similar to one of the noncentrosymmetric polymorphs of quartz iota99 ; serra99  and was designated as CO2–V. Since solid CO2 has very high tensile strength (comparable to that of boron nitride yooc99 ) this discovery could be very useful in device applications. There is also considerable interest in the role of CO2 in planetary geology hoffman01 .

Experimental work done to identify the exact crystal structure of CO2 at high pressure suggests a low symmetry tridymite phase ($`P2_12_12_1`$) yoo99 . The presence of a nonlinear optical signal also indicated the absence of inversion symmetry in the crystal structure iota99 . Most of the theoretical work has been concentrated around the lattice relaxations and total energy calculations in search of the minimum energy structure dong00I ; holm00 ; dong00II . Theory predicts the $`\beta `$-cristobalite phase ($`I\overline{4}2d`$) to be lower in energy than $`P2_12_12_1`$ by around 0.19 eV per formula unit dong00II . However, none of the SiO2-like polymorphs of CO2 gave a satisfactory fit to the experimental XRD data yoo99 ; holm00 . Other possible quartz polymorphs have been eliminated by their much higher total energies dong00I . Since all experiments were performed at elevated temperatures (1800K) and the theoretical energy difference between the two phases is small, the energy minimization is not sufficient to solve the discrepancy between theory and experiments and to predict the actual high pressure phase. In order to understand the nature of this phase transition there is the need for analyzing other physical properties and finding a satisfactory agreement between theory and experiment. Since the two probable phases have very different electronic structure it is possible to distinguish them by their optical response. So far, neither a theoretical or experimental study of the linear optical response have been performed. Although second harmonic generation (SHG) by CO2 has been seen experimentally, there exist no experimental data or calculations of the frequency dependent SHG susceptibility. The aim of the present work is to show that the linear and nonlinear optical spectra indeed allow an unambiguous identification of the CO2–V crystal structure. To this extent we have calculated the linear and second-order optical response of solid carbon dioxide, CO2–V, in the two possible phases. The linear optical spectra of these probable phases are compared to each other and also to that of the $`Cmca`$ phase, which under high pressure and temperature conditions leads to CO2–V. Since quartz is often used as a standard in second-harmonic (SH) experiments, a detailed comparison is made between the second-order optical response of the two CO2–V phases and that of SiO2 ($`I\overline{4}2d`$ phase). We propose that on basis of our results, future experimental work will help to resolve the present disagreement between theory and experiment and hence clarify the high pressure phase of CO2.

Total energy calculations are performed using the state-of-the-art full-potential linearized augmented plane wave FPLAPW method implemented in the WIEN2k code WIEN . For the exchange-correlation potential we use the generalized gradient approximation (GGA) derived by Perdew and Wang perdew92 . The detailed formalism for the determination of the linear susceptibility tensor $`\chi ^{(1)}(\omega )`$ within the FPLAPW formalism has been presented before cad . The susceptibility for the SHG $`\chi ^{(2)}(2\omega ,\omega ,\omega )`$ has been calculated using an extension to this program sl . The SHG susceptibility obtained satisfies all the theoretical sum rules presented by Scandolo and Bassani scandolo95 . All the calculations are converged in terms of basis functions as well as in the size of the $`k`$-point mesh representing the Brillouin zone. The optical properties are calculated on a mesh of 1500 $`k`$-points in the irreducible wedge of the Brillouin zone.

The phases of CO2 under investigation in the present work are the $`Cmca`$ phase and the two probable extended-solid phases for CO2–V, namely $`I\overline{4}2d`$ and $`P2_12_12_1`$. For the $`P2_12_12_1`$ phase the starting point is the experimental lattice parameters yoo99  with the atomic positions of the $`C222_1`$ phase (which is the supergroup of the space group $`P2_12_12_1`$). In a first step, the lattice parameters are kept constant at the experimental values and the internal parameters are fully relaxed. This structure is referred to as $`P2_12_12_1`$$`expt`$ in the present work. Next, the lattice and the internal parameters are simultaneously relaxed (keeping a constant volume remark\_vol ) to attain the minimum energy configuration. This structure is referred to as $`P2_12_12_1`$$`rel`$. The energy lowering on going from $`P2_12_12_1`$$`expt`$ to the $`P2_12_12_1`$$`rel`$ structure is 0.81 eV per formula unit, which is around 0.16 eV more than the previous theoretical value dong00II . The larger energy in the present calculations compared to the previous theoretical value could be due to the difference in the method of calculation, since in the previous work the pseudopotential approach was used dong00II . In case of the $`I\overline{4}2d`$ phase, due to the lack of experimental data, the starting lattice constants and internal parameters are taken to be those of SiO2 in the same phase. Then the atoms as well as the cell volume are fully relaxed to attain the minimum energy configuration. The lattice parameters for all the structures are presented in Table I. In accordance with the previous theoretical work holm00 ; dong00II  we find $`\beta `$-cristobalite ($`I\overline{4}2d`$) to be lower in energy compared to the $`P2_12_12_1`$$`rel`$ phase. The difference in energy between the $`P2_12_12_1`$$`rel`$ and the theoretically determined $`I\overline{4}2d`$ phase at 40 GPa is 0.38 eV/CO2.

The results for the linear optical spectra of CO2 in various phases are presented in Fig. 1. None of the optical spectra are scissors corrected. Both the $`P2_12_12_1`$ and the $`I\overline{4}2d`$ phase are taken at 40 GPa. The $`Cmca`$ phase shows a strong anisotropy with $`|\chi _{yy}^{(1)}(\omega )|`$ much larger in amplitude than $`|\chi _{xx}^{(1)}(\omega )|`$, and $`|\chi _{zz}^{(1)}(\omega )|`$ lying in between the two. The reason for this anisotropy is the orientation of the molecules within the crystal. The neighbouring molecules are parallel to each other in the $`ab`$ plane and form a herringbone pattern in the $`ac`$ and $`bc`$ planes \[Fig. 2(b) in Ref. yoo02, \]. This causes the $`yy`$ and $`zz`$ polarized transition matrix elements (TMEs) to be large. On the other hand the $`xx`$ polarized TMEs are close to zero, leading to small values for $`|\chi _{xx}^{(1)}(\omega )|`$. $`|\chi _{zz}^{(1)}(\omega )|`$ is smaller than $`|\chi _{yy}^{(1)}(\omega )|`$ because the lattice parameter $`b`$ is smaller than $`c`$. On going from the low pressure $`Cmca`$ phase to any other phase of CO2-V, $`|\chi _{yy}^{(1)}(\omega )|`$ and the anisotropy between $`|\chi _{yy}^{(1)}(\omega )|`$ and $`|\chi _{xx}^{(1)}(\omega )|`$ decrease dramatically. The latter fact can again be explained by the orientation of the molecules inside the unit cell \[Fig. 3E and 3F in Ref. holm00, \]. The neighbouring CO2 molecules are no more parallel or perpendicular in any direction leading to a weakening of the very strong transitions and an enhancement of the TMEs which were forbidden in the $`Cmca`$ phase. For all the probable phases of CO2–V the spectral weights are shifted to higher energies. The optical band gap of the $`P2_12_12_1`$$`expt`$ structure is smaller than that of the $`Cmca`$ phase, while the other two structures show almost equal optical band gaps, which are higher than that of the $`Cmca`$ phase. Compared to the $`Cmca`$ phase the maximum peak height is around 3 times smaller in the $`I\overline{4}2d`$ phase, while it is about 4 times smaller in the $`P2_12_12_1`$ structures. Beyond 11 eV, the response of all the probable phases of CO2–V is on average similar in magnitude, but differs in the amount of anisotropy.

For $`I\overline{4}2d`$ there are two independent tensor components. Compared to the $`zz`$ component, the $`xx`$ component is approximately upto 2 times as large in the energy range below 10 eV. Upto 10.5 eV and between 11.5 and 13 eV the effect is opposite, while the two are nearly the same for higher energies. For the $`P2_12_12_1`$ structures there are three inequivalent components and the anisotropy between the $`xx`$ and $`zz`$ components is very small compared to the $`I\overline{4}2d`$ phase. We also find the lattice relaxation effects which have already been described in detail by Dong et al. dong00II . They have some consequences for the optical spectra. On going from the $`Cmca`$ phase to $`P2_12_12_1`$$`expt`$, the optical band gap decreases, while it increases in case of $`P2_12_12_1`$$`rel`$. Another noticeable feature is that the anisotropy is much larger in the $`P2_12_12_1`$$`expt`$ structure. Since CO2–V is not quenchable into ambient pressure at room temperature the optical measurements need to be performed with the sample inside the diamond anvil cell (DAC). The band gap of diamond is 5.5eV. Therefore the energy regime below 5.5eV is most interesting, since the optical response in this energy regime is purely from the sample inside the DAC. The difference between the linear optical spectra of various probable phases of CO2 below this energy range is not as pronounced as in the higher energy regime. The results though are of theoretical interest and help in identification of various features in the nonlinear optical spectra (which is rather difficult to do from the band structure plots), but may not be very good in phase distinction under these experimental conditions.

The nonlinear optical properties are much more sensitive to small changes in the band structure than the linear optical spectra. The calculated magnitude of the SHG susceptibility by CO2–V in various probable phases are presented in Fig. 2. Due to the presence of inversion symmetry in the $`Cmca`$ phase the SHG susceptibility is zero. The second-order response of SiO2 ($`I\overline{4}2d`$), which is often used as standard in SHG experiments, is also presented in the figure. In the energy range between 3 and 5.5 eV the magnitude of the SHG susceptibility by CO2 in the $`P2_12_12_1`$$`expt`$ phase is 1.5–2 times smaller than that of quartz. On the other hand the $`I\overline{4}2d`$ phase of CO2 shows a SHG susceptibility 3–4 times larger than that of SiO2. These pronounced differences make the phases easily identifiable. The energy range of 3–5.5 eV lies below the band gap of diamond and corresponds to laser light of wavelengths between 413.3 - 247.9 nm. Previous SH experiments were performed using laser light of 1054 nm wavelength iota99 , which correspond to 1.176 eV. At this energy our calculated value of $`|\chi _1^{(2)}(2\omega ,\omega ,\omega )|`$ for CO2 in the $`I\overline{4}2d`$ phase is 2.7 times larger than that of SiO2. While, the SHG susceptibility by CO2 in the $`P2_12_12_1`$$`expt`$ phase is 2.6 times smaller than that of SiO2. The previous experimental work yoo99  reports the SH conversion efficiency of the CO2–V phase, which is not just proportional to the square of the SHG susceptibility but also depends on the input intensity and sample parameters. So no direct comparison of the existing experimental data can be made with the present results. Since the difference in the second-order response by the probable phases of CO2-V are more pronounced in the energy range between 3 and 5.5 eV than in between 1–2 eV, experiments at higher energies would be favourable. Experimental determination of the SHG susceptibility by CO2–V and SiO2, over some frequency range, and a comparison with our work could lead to conclusive results.

The two structures of the $`P2_12_12_1`$ phase are also differentiable by their SHG response. Due to the large optical band gap of $`P2_12_12_1`$$`rel`$ the SHG susceptibility is approximately 3 times smaller than that of $`P2_12_12_1`$$`expt`$. As shown in Ref. dong00II,  the XRD patterns for neither of these structures matches the experiments. The theoretically calculated XRD pattern of $`P2_12_12_1`$$`rel`$ shows a better agreement with the experimental XRD pattern at higher values of 2$`\theta `$ (9–13 degrees), while the XRD pattern for $`P2_12_12_1`$$`expt`$ is in better agreement at low values of 2$`\theta `$ (6–8 degrees). It is possible that under the experimental conditions of high temperatures the material is a polycrystalline mixture of these two structures or has a configuration in between the two, in which case the difference between the optical spectra of the two phases would represent a measure of the uncertainty in the optical response of this phase. Despite this fact, we want to point out that in general the differences in the optical response between the two phases are much more pronounced than “the lattice relaxation effects” for the $`P2_12_12_1`$ phase making them distinguishable by optical probes.

The optical coefficients in the static limit are directly measurable quantities. The calculated inequivalent components of $`ϵ_1(0)`$ and $`|\chi ^{(2)}(0)|`$ are presented in the Table II.

In the $`Cmca`$ phase, the anisotropy between $`ϵ_1^{xx}(0)`$ and $`ϵ_1^{yy}(0)`$ is 54%, and between $`ϵ_1^{xx}(0)`$ and $`ϵ_1^{zz}(0)`$ it is 41%. In the $`I\overline{4}2d`$ phase, the anisotropy between the $`xx`$ and $`zz`$ components is 6% only. For the tridymite phase at experimental lattice parameters ($`P2_12_12_1`$$`expt`$) the maximum anisotropy found between $`ϵ_1^{yy}(0)`$ and $`ϵ_1^{zz}(0)`$ (28 %), while between $`ϵ_1^{xx}(0)`$ and $`ϵ_1^{zz}(0)`$ it is only 1%. In the case of $`P2_12_12_1`$$`rel`$ the anisotropy between the $`xx`$ and $`yy`$ components is 7%, and between the $`xx`$ and $`zz`$ components is close to zero.

Amongst the phases under investigation in the present work the anisotropy in the $`xyz`$ and $`zyx`$ components of $`|\chi ^{(2)}(0)|`$ is maximum for $`P2_12_12_1`$$`expt`$ with 86 % and minimum for $`I\overline{4}2d`$ with 7%. $`P2_12_12_1`$$`rel`$ lies in between the two with an anisotropy of 67 %. For SiO2 in the $`I\overline{4}2d`$ phase, our calculated values of $`|\chi _{xyz}^{(2)}(0)|`$ and $`|\chi _{zyx}^{(2)}(0)|`$ are 0.28 and 0.20 $`\times 10^8`$ esu, respectively. It may be noted that the anisotropy as well as the magnitude of $`\chi _1^{(2)}(0)`$ for CO2 in the $`P2_12_12_1`$$`expt`$ phase are close to that of SiO2. The fact that $`|\chi ^{(2)}(0)|`$ of CO2 is around 3 times larger in the $`I\overline{4}2d`$ phase, and 5–10 times smaller in the $`P2_12_12_1`$$`rel`$ phase compared to SiO2 makes the distinction between the two phases clear and hence the identification of the high pressure phase feasible.

In summary we find that the optical properties are excellent tools for an unambiguous determination of the high pressure phase of CO2. Particularly, the second-order response of CO2–V compared to that of SiO2 can be used for phase determination even with the sample inside the DAC. Experiments of this kind and a detailed comparison with the present theoretical results would shed light on the existing discrepancy between theory and experiments regarding the space group of the CO2–V phase.

The work is supported by the Austrian Science Fund (projects P13430, P14004 and P16227) and the EXCITING network funded by the EU (Contract HPRN-CT-2002-00317). One of the authors (SS) would also like to thank Dr. Valentin Iota (Lawrence Livermore National Lab., California) for valuable comments and suggestions.


text:
# Further Evidence for Neutrino Flux Variability from Super-Kamiokande Data

## Abstract

While KamLAND apparently rules out Resonant-Spin-Flavor-Precession (RSFP) as an explanation of the solar neutrino deficit, the solar neutrino fluxes in the Cl and Ga experiments vary with solar rotation rates. Added to this evidence, summarized here, a power spectrum analysis of the Super-Kamiokande data reveals (99.9% CL) an oscillation in the band of twice the equatorial rotation frequencies of the solar interior. An $`m=2`$ magnetic structure and RSFP, perhaps as a subdominant process, would give this effect. Solar cycle data changes are seen, as expected for convection zone modulations.

Recent results from the KamLAND experiment ref:1  seem to confirm the Large-Mixing Angle (LMA) solution to the solar neutrino deficit and rule out the Resonant-Spin-Flavor-Precession (RSFP) explanation ref:2 . On the other hand, there is increasing evidence ref:3 -ref:9  that the solar neutrino flux is not constant as assumed for the LMA solution, but rather it varies with well known solar rotation periods. The solar neutrino situation may be complex, and Spin-Flavor-Precession could be subdominant to LMA, as suggested in ref:10 , or if there is at least one light sterile neutrino even RSFP could be a subdominant process. Since this recent information on solar neutrino variability is not widely known, a brief summary is presented of analyses of radiochemical neutrino data, along with new input from the Super-Kamiokande experiment ref:11 . While the 10-day averages of Super-Kamiokande solar neutrino data ref:12  show no obvious time dependence, a power-spectrum analysis ref:13a  displays a strong peak at the frequency $`26.57\pm 0.05`$ y-1 (period 13.75 d), where the width of peaks are computed for a probability drop to 10%. The probability of finding this peak (or a stronger peak) by chance within a band specified by twice the near-equatorial rotation frequencies of the solar interior (25.36-27.66 y-1) is found to be 0.001 (99.9% CL). This frequency is also seen in the radiochemical data.

Whether or not it is the correct interpretation, the RSFP framework provides a simple way to understand the data, and hence it will be used here. The RSFP mechanism, which requires a neutrino transition magnetic moment, provides at least as good a global fit ref:13 ; ref:14  to solar data (which depends mainly on the Super-Kamiokande spectrum) and a better fit to average rates of the individual experiments than does LMA. This is because the neutrino survival probability, while having a resonance pit at a density that suppresses the 0.86 MeV 7Be line (as does the Small-Mixing-Angle (SMA) solution), tends at high energies toward 1/2 and hence fits the Super-Kamiokande spectrum, whereas the survival probability goes to unity in the SMA case. Thus, as a subdominant process, it could also improve fits to the data.

A brief review follows of the published evidence for solar neutrino flux variability, put into a coherent neutrino scheme. By analyzing $`10^3`$ simulated data sequences, it was found ref:3  that the variance of the Homestake solar neutrino data ref:15  is larger than expected at the 99.9% CL. A power spectrum analysis ref:3  of the data showed a peak at $`12.88\pm 0.02`$ y-1 (28.4 d), compatible with the rotation rate of the solar radiative zone. Four sidebands gave evidence at the 99.8% CL for a latitudinal effect associated with the tilt of the sun’s rotation axis, and the latitude dependence was also seen directly in the data at the 98% CL ref:4 . The GALLEX data ref:16  showed a peak at $`13.59\pm 0.04`$ y-1, compatible with the equatorial rotation rate of the deep convection zone. This peak is also in the Homestake data, and a combined analysis of both data sets shows that the 13.59 y-1 peak is larger than in either data set alone. Comparison ref:7  of the power spectrum for the GALLEX data with a probability distribution function for the synodic rotation frequency as a function of radius and latitude, derived from SOHO/MDI helioseismology data ref:17 , results in Fig. 1. This map shows the rotation frequency coincides with the neutrino modulation in the equatorial section of the convection zone at about 0.8 of the solar radius, $`R_{}`$. The influence of these rotation frequencies extends even to the corona, since the SXT instrument ref:18  on the Yohkoh spacecraft provides X-ray evidence for two “rigid” rotation rates, one ($`13.55\pm 0.02`$ y-1) mainly at the equator, and the other ($`12.86\pm 0.02`$ y-1) mainly at high latitudes. These values are in remarkable agreement ref:8  with the neutrino modulation frequencies and also with their equatorial location in one case and non-equatorial in the other; see Fig. 2.

That two dominant frequencies are shown by both coronal X-rays and neutrino flux is a feature of magnetic rotation well known in solar physics. For instance, an analysis ref:19  of the photospheric magnetic field during solar cycle 21 found two dominant magnetic regions: one in the northern hemisphere with synodic rotation frequency $`13.6`$ y-1, and the other in the southern hemisphere with synodic rotation frequency $`13.0`$ y-1. Similarly, an analysis ref:20  of flares during solar cycle 23 found a dominant synodic frequency of $`13.5`$ y-1 for the northern hemisphere and $`12.9`$ y-1 for the southern hemisphere. These and other studies show a strong tendency for magnetic structure to rotate either at about 12.9 y-1 or at about 13.6 y-1. Since the sun’s magnetic field is believed to originate in a dynamo process at or near the tachocline, it is possible that some magnetic flux is anchored in the radiative zone just below the tachocline, where the synodic rotation frequency is about 12.9 y-1, and some just above the tachocline in the convection zone, where the synodic rotation frequency is about 13.6 y-1. It is also possible that the 12.9 y-1 frequency results from a latitudinal wave motion in the convection zone excited by structures at or near the tachocline.

An example of latitudinal oscillatory motion of magnetic structures may be the well-known Rieger-type oscillations with frequencies of about 2.4, 4.7, and 7.1 y-1 ref:21 . These periodicities may be attributed to r-mode oscillations with spherical harmonic indices $`\mathrm{}=3`$, $`m=1`$,2,3, giving frequencies $`\nu =2m\nu _R/\mathrm{}(\mathrm{}+1)`$ which are seen in neutrino data ref:3 ; ref:5 ; ref:9 . A joint spectrum analysis ref:9  of Homestake and GALLEX-GNO data yields peaks at 12.88, 2.33, 4.62, and 6.94 y-1, indicating a sidereal rotation frequency $`\nu _R=13.88\pm 0.03`$ y-1. While $`l2`$ is required, odd$`l`$ values have nonzero poloidal velocity at the equator that could move magnetic regions in or out of the neutrino beam to earth.

The difference between the main modulations detected by the Cl and Ga experiments may be explained by the tilt of the solar axis relative to the ecliptic, along with the fact that Cl and Ga neutrinos are produced mainly at quite different radii. The Ga data comes mostly—especially as the fit requires suppressing the 7Be line—from $`pp`$ neutrinos, which originate predominantly at large solar radius ($`0.2`$ R), so that the wide beam of neutrinos detected on earth is insensitive to axis tilt. Thus the beam of neutrinos detected by the Ga experiments exhibits no seasonal variation and can be modulated by the equatorial structure indicated in Fig. 1, leading to the observed frequency at about 13.6 y-1. On the other hand, the Homestake experiment detects neutrinos produced from a smaller sphere ($`0.05`$ R), so that the axis tilt causes these neutrinos mainly to miss the equatorial structure of Fig. 1 and instead sample nonzero latitudes where the 12.9 y-1 modulation may occur. Twice a year axis tilt has no effect for these neutrinos, leading to a seasonal variation in the measured flux ref:4 .

The 13.6 y-1 frequency, located as in Fig. 1, represents a modulation of the $`pp`$ neutrinos which are at or near the steeply falling edge of the neutrino survival probability. This is close to the RSFP resonance pit suppressing the $`{}_{}{}^{7}\mathrm{Be}`$ neutrinos, which is where the largest value of $`\mathrm{\Delta }m^2/E`$ satisfies

$$\mathrm{\Delta }m^2/E=2\sqrt{2}G_FN_{\mathrm{eff}},$$
(1)
and is essentially the same as for an MSW resonance ref:23 . For MSW, $`N_{\mathrm{eff}}=N_e`$ (the electron density), and for RSFP $`N_{\mathrm{eff}}=N_eN_n`$ for Majorana neutrinos, where $`N_n`$ is the neutron density (about $`N_e/6`$ in the region of interest in the Sun). For Dirac neutrinos $`N_{\mathrm{eff}}=N_eN_n/2`$, but only Majorana neutrinos provide a fit to solar experimental rates with RSFP ref:13 ; ref:14 , a very important consequence of such a solar solution. Using the $`N_e`$ and $`N_n`$ values at the 0.8 R location of Fig. 1 in Eq. 1 results in $`\mathrm{\Delta }m^2/E10^{14}`$ eV, putting the RSFP resonance pit at the location needed to fit the solar data. $`N_{\mathrm{eff}}`$ varies exponentially with radius and would be quite different for a somewhat changed neutrino modulation frequency.

It is well known that the convection zone magnetic field changes with the solar cycle, so the neutrino modulation features described above are not permanent. It has been argued ref:24  that an RSFP effect would have to be in the radiative zone, where the field does not change with the solar cycle, under the assumption solar cycle variations are not observed. On the contrary, we show here that solar cycle changes play an important role. Variations in neutrino rates are difficult to observe, since changes in field magnitude would be undetectable if the transition remains adiabatic, and even if flux modulation results, average rates may vary only slightly. The feature that is most sensitive to field magnitude or radial variations is the intersection of the very steeply falling $`pp`$ neutrino spectrum with the falling edge of the RSFP resonance pit. As can be seen in Fig. 3, the 13.6 y-1 frequency associated with these neutrinos increased in strength from the start of data taking after the solar maximum of 1989.6 to the solar minimum of 1996.8, after which the modulation becomes weak. Also it was during that cycle that the main buildup in the strength of the 12.9 y-1 oscillation was detected by Homestake. The SXT X-ray data, with which these two frequencies had remarkable agreement ref:8 , also came from that same solar cycle.

Another feature attributable to a time variation of the intersection of the $`pp`$ spectrum with the edge of the RSFP resonance pit is the appearance of a bimodal flux distribution in the Ga data. This distribution is clearly evident ref:6  during the same solar maximum to solar minimum period when the flux is binned appropriately using individual runs, instead of averages over several data runs. This neutrino flux effect also diminishes after the solar minimum. Adding to the significance of this result is the plot of Fig. 4, which shows that when the end times of runs are reordered according to the phase of the 13.59 y-1 rotation, the flux values are low in one-half of the cycle and high in the other.

It is unfortunate that the Homestake experiment, which detected mainly the 8B neutrinos (the only neutrino component registered by Super-Kamiokande) stopped operating at about the time that Super-Kamiokande started. As a result, there is no way to predict from results of other experiments what neutrino flux variation should have been detected by Super-Kamiokande during its operation from May 1996 (near solar minimum) until July 2001 (near solar maximum).

Recently the Super-Kamiokande group released ref:12  flux measurements in 184 bins of about 10 days each. While these measurements vary by $`2`$, averaged over all bins their fractional error is $`0.14`$. Because of the regularity of the binning, the “window spectrum” (the power spectrum of the acquisition times) has a huge peak (power $`S>120`$) at a frequency of $`\nu =35.98`$ y-1 (period 10.15 d). (Note that the probability of obtaining a power of strength $`S`$ or more by chance at a specified frequency is $`e^S`$.) This regularity in binning leads to aliasing of the power spectrum of the flux measurements. In a spectrum formed by a likelihood procedure, the strongest peak in the range 0 to 100 y-1 occurs at $`\nu =26.57`$ y-1 ref:13a ; ref:25  with $`S=11.26`$. In the range 0–40 y-1, the next strongest peak is at $`\nu =9.41`$ y-1 with $`S=7.33`$. Since $`26.57+9.41=35.98`$, we infer that the weaker peak (9.41 y-1) is an alias of the stronger (26.57 y-1).

Having found this frequency in the high-statistics Super-Kamiokande data, we then examined the power spectrum from a combined analysis of the SAGE ref:26  and GALLEX-GNO data, the sampling time of the Homestake data being too long for this high a frequency. The combined experiments showed a peak at 26.54$`\pm 0.05`$ y-1 with $`S=5.75`$. When combined with the Super-Kamiokande data, we find a peak at 26.54$`\pm 0.03`$ y-1 with $`S=13.95`$, as shown in Fig. 5. There is no significant peak near 9.41 y-1, further showing this to be an alias frequency.

Since 26.5 y-1 is within the band of twice the near-equatorial rotation frequencies in the solar interior, this result indicates that some 8B neutrinos are experiencing an “$`m=2`$” structure, two circumferential regions of the magnetic field which differ in strength from the rest. This modulation, unlike the Ga $`pp`$ threshold effect, is not very deep ($``$ 10%). Such oscillations at the harmonic of the rotation frequency are not uncommon in solar data.

The analysis of Super-Kamiokande data, when combined with the results of analyses of Cl and Ga data, yields evidence for variability of the solar neutrino flux. The KamLAND experiment supports the LMA solution, but it could be combined with RSFP. A definitive result on neutrino flux time dependence can be carried out using a combined analysis from all experiments, particularly using one-day bins from all water Cerenkov detectors.

We are indebted to many friends for assistance and helpful discussions. D.O.C. is supported by a grant from the Department of Energy No. DE-FG03-91ER40618, and P.A.S. is supported by grant No. AST-0097128 from the National Science Foundation.