text:
# Modelling of SARS for Hong Kong

## I Susceptible-Infected-Recovered model

In the SIR model, the fixed population $`N`$ is divided into three distinct groups: Susceptible $`S`$, Infected $`I`$ and Removed $`R`$. Those at risk of the disease are susceptible, those that have it are infected and those that are either quarantined, dead, or have acquired immunity are removed. The following flow chart shows the basic procession of the SIR model Daley1 .

$$\begin{array}{ccccc}& rSI& & aI& \\ S& & I& & R\end{array}$$
(1)
Here $`r`$ and $`a`$ are the infection coefficient and removal rate, respectively. A discrete model was adapted by Gani from the original SIR model through the coarse-graining process and was applied to successfully predict outbreak of influenza epidemics in England and WalesGani1 ; Spicer1 .

$`S_{i+1}`$ $`=`$ $`S_i\left(1rI_i\right)`$
$`I_{i+1}`$ $`=`$ $`\left[rS_i+\left(1a\right)\right]I_i`$ (2)
$`R_{i+1}`$ $`=`$ $`R_i+aI_i`$
During the epidemic process, $`N=S_i+I_i+R_i`$ is fixed.

Initially, we examine the data for the first 15 days to estimate the parameters $`r`$ and $`a`$ of SARS for HK. The only data are new cases (removal $`R`$) announced everyday by HK Dept. of Health from March 12, 2003 followed by a revised version laterHK1 . To avoid inadvertently using future information we do not use the revised data at this stage. Population $`N=7.3\times 10^6`$; since $`I+RN`$ it is reasonable to let $`S_i=N`$ in right hand side of (2). $`I_1=1`$ and $`R_1=0`$ is set as the initial condition. $`I_i`$ is replaced with $`R_{i1}`$ whereas $`I_i`$ is uncertain. This assumption implies the incubation period is only one day, in spite of the fact that the true incubation period of the coronavirus is 2-7 daysRota1 . The parameters $`r`$ and $`a`$ are scanned for the best fit for the stage. For every $`(r,a)`$, a sequence of $`R_{i+1}^{}`$ is obtained by numerical simulations. An Euclidean norm of $`_{i=1}^{15}\left(R_{i+1}R_{i+1}^{}\right)^2`$, which indicates a distance between the true and simulated data, is applied to measure the fit. In Fig. 1 $`(r,a)`$ of the lowest point is the best fit parameters for this stage and this value is used for the following prediction. We get $`r=2.05\times 10^7`$ and $`a=1.444`$. The method is applied to get parameters of $`a`$ and $`r`$ in Figs. 2 and 3.

The prediction is available for the trend based on parameters $`r`$ and $`a`$ of this stage, the middle day (March 20, 2003) of the stage is applied as the first day and $`\overline{R}=\frac{1}{15}_{i=1}^{15}R_{i+1}`$ is assumed as $`\overline{I}`$. A curve of squares is plotted in Fig. 2 for the first stage prediction. In the same way the best fit and prediction is applied for the next two stages and plotted in Fig. 2. The best fit is also processed weekly for detail discussion later.

The curves for prediction clearly shows the 3 stages in Fig. 2. The first stage (March 12 - 27, 2003) exhibits dangerous exponential growth. It shows that the extremely infectious SARS coronavirus spread quickly in the public with few protections during this early stage. More seriously, it leads to a higher infection peak although in this stage the averaged number of new cases $`\overline{R}`$ is below 30. This prediction gets confirmation in the second stage (March 28 - April 11, 2003). The peak characterized by the Amoy Garden outbreak comes earlier and higher. It indicates appearance of a new transmission mode that differs from intimate contact route observed in the first stage. We name the new transmission mode $`explosive`$ growth in the contrast to transmission in the first stage which we refer to as $`burning`$ growth. Outbreaks at Prince Wales Hospital (PWH) (where SARS patients received treatment) and Amoy Gardens (a high-density housing estate in Hong Kong) represent infections of $`burning`$ and $`explosive`$ modes in the full epidemic, respectively. The revised data provided by HK Dept. of Health in Fig. 2 present a clearer indication for these two transmission modes. The outbreak of PWH sent the clear message that intimate contacts with SARS patient led to infectionTomlinson1 . Detailed investigations of the propagation at Amoy Garden suggested that faulty of sewerage pipes allowed droplets containing the coronavirus to enter neighboring units vertically in the buildingRiley1 . Furthermore, poor ventilation of lifts and rat infestation where also suggested as possible modes of contamination Normile1 ; Ng1 ; Helen1 ; Donnelly1 . However, control measures bring the spread of the disease under control in the second stage. Contrary to the increasing trend in the first stage, the prediction curve of the 2nd stage (circle) declines. And it anticipates that new cases drop below 10 before the 60th day (May 10, 2003). Also on April 12, 2003 we predicted number of whole infection cases would reach 1700. Up to April 11, 2003 there were 32 deaths and 169 recoveriesHK1 . We calculated the mortality of SARS as the ratio of death to sum of both deathes and recoveries, and it was 15.9%. Therefore we predicted approximately 270 fatalities in total. In the third stage (April 12 - 27, 2003), triangle curve in Fig. 2 refines results and gives more accurate prediction. This stage predicts that new cases per day drops to 5 before the 62th day, May 12, 2003. The travel warning for HK was cancelled by WHO because HK had kept new cases below 5 for 10 days since May 15, 2003. And, finally, we predict that the total cases reaches 1730 and nearly 287 deaths (up to April 27, 2003 there were 668 recovered and 133, the mortality increased to 16.6%). These numbers are very close to the true dataHK1 . Precisely the method drawn from the SIR model has been verified for prediction for full epidemic. However, the accuracy is only possible by the first dividing the epidemic into separate ”stages”.

This problem of determining to what extent an epidemic is under control is of greater strategic significance. Information on the efficacy of epidemic control will help determine whether to apply more control policies or not, and balance cost and benefit from them. For each individual the same question will also inform the degree to which precautions are taken: i.e. wearing a surgical mask to prevent the spread (acquisition) of SARS. Obviously a way to estimate control level is required. Actually this is a difficult problem because of significant statistical fluctuations in the data. A quite simple method to evaluate control efficiency is discussed below.

In the SIR model (2), if $`\mathrm{\Delta }I=I_{i+1}I_i0`$ a disease is regards as being controlled as new cases will decrease. This inequality leads to a control criterion for some diseases in epidemiology research. Applying the approximation of $`S_i=N`$ we get a threshold $`Nra`$ from (2). We rescale $`rNr`$ ($`r`$ is called infection rate in place of infection coefficient now) and then get the threshold that is free from population $`N`$.

$$ra$$
(3)
This indicates that the removal rate $`a`$ exceeds the infection rate $`r`$. In Fig. 3 the circles show the 3 stages evaluations with dash line of $`r=a`$. And the parameters $`r`$ and $`a`$ of SARS data for HK are estimated weekly as squares also in Fig. 3. The line of $`r=a`$ is regarded as the critical line since number of infected cases increases when $`(r,a)`$ passes through it from below. It is possible to apply the diagram of $`r`$ and $`a`$ to compare control level for different countries and areas even for different diseases. This provides organizations like WHO with a simple and standard method to supervise infection level of any disease. The limitation of the method comes from the assumptions of the SIR model. More accurate models may provide better estimate of epidemic state and future behaviours.

In summary, a discrete SIR model gives good predictions for epidemic of SARS for HK. Two distinct methods are described for the disease propagation dynamics. Particularly the $`explosive`$ mode is much more hazardous than the $`burning`$ mode. We have introduced a simple method to evaluate control levels. The method is generic and can be widely applied to various epidemiological data.

## II The small-world model

Contrary to the long established SIR model, epidemiology research using Small-World (SW) network models is young and growing area. The concept of SW was imported from the study of social network into nature sciences in 1998Watts1 . However, it provides a novel insight for networks and arouses a lot of explorations in the brain, social networks and the Internet.Laughlin1 ; Albert1 ; Ebel1 ; Liljeros1 . Some SW networks also exhibit a scale-free (power law) distribution: a node in this network has probability $`P(k)k^\gamma `$ to connect $`k`$ nodes, $`\gamma `$ is one basic character for the systemAlbert1 ; Newman1 ; Ebel1 ; Liljeros1 . Most researches of “virus” spreading with the SW model, concern the propagation of computer viruses on the Internet, however, a few studies have been published relating to epidemiologyEbel1 ; Liljeros1 ; Strogatz1 ; Kuperman1 ; Huerta1 . The dynamics of SW models enriches our realization of epidemic and possibly provides better control policiesMiramontes1 . From the first SARS patient to the last one, an epidemic chain is embedded on the scale free SW network of social contacts. It is of great importance to discover the underlying structure of the epidemic network because a successful quarantine of all possible candidates for infection will lead to a rapid termination of the epidemic. In HK e-SARS, an electronic database to capture on-line and in real time clinical and administrative details of all SARS patients, provided invaluable quarantine information by tracing contactsHK1 ; Brower1 . Unfortunately a full epidemic network of SARS for HK is still unavailable. Because data representing the underlying network structure is currently not available, we have no choice other than numerical simulations. Therefore, our analysis of the SW model is largely theoretical. The only confirmation of our model we can offer is that the data appear to be realistic and exhibit the same features as the true epidemic data.

To simulate an epidemic chain, a simple model of social contacts is proposed . The model is established on a grid network weaved by $`m`$ parallel and $`m`$ vertical lines. Every node in the network represents a person. We set $`m=2700`$ with population $`N=m^2=7.29\times 10^6`$. All nodes are initiated with a value of 0 (named $`good`$ nodes). Every node has 2, 3 or 4 nearest neighbors as short range contacts for corner, edge and center, respectively. For every node there are two long range contacts with 2 other nodes randomly selected in the whole system everyday. These linkages model the social contact between individuals (i.e. social contacts that are sufficiently intimate to bring individuals at risk of spreading the disease). One random node of the system is set to 1 (called the $`bad`$ node), through its short and long range contacts, the value of the nodes linked with it turns into 1 according to probability of $`p_1`$ and $`p_2`$, respectively. An infection happens if a node changes its value from 0 to 1. This change is irreversible. During the full simulation process, the bad node is not removed from the whole system. We make this assumption because the number of deaths is small in comparision of the population, and there is no absolute quarantine—even the SARS patient in hospital can affect the medical workers. Moreover, the treatment period for SARS is relatively length, and during this time infected indviduals are highly infectious. To reflect the true variation in control strategy and individual behaviour, the control parameters $`p_{1,2}`$ ($`0p_{1,2}1`$) vary with time.

In same model it would be interesting to compare with the complete infection epidemic. A small system with $`100\times 100`$ nodes is chosen. A fixed probability $`p_1=p_2=0.05`$ leads to eventual infection of the entire population since all nodes are linked and the infected are not removed. The epidemic curve for this process is plotted in Fig. 4 a) and is typical of many plagues. In Fig. 4 b) various sizes of clusters with infected nodes (black dots) scatter over the geographical map. It contains all infection facts in first 45 days. Comparing to the true SARS infection distribution in HK only slight similarity is observed. The short and long range linkages give good infection dynamics as we expected. The epidemic chain is easily drawn by recording infection fact from the seed to last patient. However, during simulations there is a problem: for a $`good`$ node linked by more than one $`bad`$ nodes, which $`bad`$ node infects it? A widely accepted preferential attachment of rich-get-richer is a good answerAlbert1 ; Newman1 ; Watts2 . A linear preference function is applied here. With presence of both growth and preferential attachment, it is general to ask whether the chain is a scale-invariant SW network. We plot the distribution in Fig. 5 a) and b) in log-log and log-linear coordinates, respectively. The hollow circle is for the system with $`100\times 100`$ nodes. The scaling behaviour of it looks more like a piecewise linear (i.e. bilinear) in b) rather than a power law in a). To confirm this we tried a larger system with $`1000\times 1000`$ nodes and the same fixed $`p_{1,2}`$ in the rest curves in Figs. 5. The solid square curve is distribution of the whole epidemic network. The linear fit for the solid square gives a correlation coefficient $`R`$ of $`0.956`$ in Fig. 5 a). In b) piecewise linear fits have $`R`$ of $`0.997`$ and $`0.994`$. These provide positive support for the original model. The ratio of the cases of first 30% and 1% of whole process for this bigger system are plotted as hollow triangle and solid circle curves in Figs. 5, respectively. The case of 1%, early stage of full infection epidemic, suggests a better fit for scale-free than the other cases since it has correlation coefficients $`R`$ of $`0.983`$ and ratio $`\gamma =2.89`$ for linear fit in a). For the solid circle curve piecewise linear fits give $`R`$ of $`0.995`$ and $`0.977`$ in Fig. 5 b). So, what scaling behaviour is true during full infection process? With absence of rigorous proof this problem is hard to answer correctly in simulations. We may only draw conclusions based on which simulations most closely matches the qualitative features of the observed data.

Let’s return to the system with $`2700\times 2700`$ nodes for modelling SARS for HK. The probability $`p_{1,2}`$ is fitted to the true epidemic data. However, it is fruitless to obtain an exact coincidence between the simulated results and the true data as the model evolution is highly random (moreover this would result in overfitting). The control parameters $`p_{1,2}`$ (dots and dashes curves) and a simulated epidemic curve (black dots) with column diagram of SARS for HK is plotted in Figs. 6 a) and b), respectively. For the simulated data, the total number of cases is 1830 that has a 4.3% deviation from the true data of 1755. Contrary to the above full infection with fixed parameters, $`p_{1,2}`$ are believed to drop exponentially and lead to small part infection epidemic without quarantine or removal. In any case, the high probability of infection in the early stage is indicative of the critical ability of the SARS coronavirus to attack an individual without protection. If SARS returns, the same high initial infection level is likely to occur. The only hope to avoid a repeat of the SARS crisis of 2003 is to shorten the high infection stage by quick identification, wide protection and sufficient quarantining. In other words, the best time to eliminate possible epidemic is the moment that the first patient surfaces. Any delay may lead to a worsening crisis. The long range infections in the model and the world also imply an efficient mechanism to respond rapidly to any infectious disease is required to establish global control.

Data on the geographical distribution of SARS cases in HK is much easier to collect than the full epidemic chain. Numerical simulations provide both simply. The full geographical map marked with all infected nodes (black dots) and an amplified window of a cluster is plotted in Figs. 7 a) and b), respectively. Similarity is expected and verified in cluster patterns of Figs. 4 b). The scaling behaviour of the epidemic chain is plotted in Figs. 8. The curve in log-log diagram exhibits a power-law coefficient of $`\gamma =3.55`$ and gives the linear fit correlation coefficient $`R=0.989`$. The piecewise linear fit for the log-linear case gives correlation coefficients $`R`$ of $`0.987`$ and $`0.959`$, respectively.

For a SW network often few vertexes play more important roles than othersWatts2 . SARS super-spreaders found in HK, Singapore and China are consistent with thisRiley1 ; Leo1 . Data for the early spreading of SARS in SingaporeLeo1  show definite SW structure with a small number of nodes with a large number of links. The average number of links per node also shows a scale-free structure Leo1 , but the available data is extremely limited (the linear scaling can only be estimated from three observations).

This character is also verified in our model. The first few nodes have a high chance of infecting a large number of individuals. In Fig. 8, a single node has 40 links. Clearly the index node has many long range linkages. It has been suggested that travelling in crowded public places (train, hospital, even an elevator) without suitable precautions can cause an ordinary SARS patient to infect a significant number of others. Again, this is an indication that to increase an individual’s (especially a probable SARS patient’s) personal protection is key to rapidly control an epidemic. Actually, in our model, if duration of the early stage with high probability $`p_{1,2}`$ is reduced to less than 10, the infection scale decreases sharply.

An engrossing phenomenon is the points of inflection in curves in log-linear diagrams of Figs. 8 and Figs. 5 b). All are located near linkages number of 6-7. On average a nodes has about 6 contacts (2-4 short plus 2 long range) every day, although in whole process there is no limitation for linkages. For a growing random network, a general problem always exists among its scaling behaviour, preferential attachment and dynamics, even embedding a geographical mapCohen1 ; Warren1 ; Krapivsky1 ; Dorogovtsev1 ; Moore1 . More work is required to address this issue.

In conclusion a SW epidemic network is simulated to model SARS spreading in HK. A comparison of the simulations with full infection data is presented. Our discussion of the infection probability and occurrence of super spreaders lead to the obvious conclusion: rapid response of an individual and government is a key to eliminating an epidemic with limited impact and at minimal cost.

###### Acknowledgements.

This research is supported by Hong Kong University Grants Council CERG number B-Q709.


format: text

text:
# Mutual Experimental and Theoretical Validation of Bulk Photoemission Spectra of Sr1-xCaxVO3

Present Address: \]National Institute of Advanced Industrial Science and Technology (AIST), Umezono, Tsukuba, Ibaraki 305-8568, Japan

## Abstract

We report high-resolution high-energy photoemission spectra together with parameter-free LDA+DMFT (local density approximation \+ dynamical mean-field theory) results for Sr1-xCaxVO3, a prototype 3$`d^1`$ system. In contrast to earlier investigations the bulk spectra are found to be insensitive to $`x`$. The good agreement between experiment and theory confirms the bulk-sensitivity of the high-energy photoemission spectra.

Electronic correlations in transition metal oxides (TMO) are at the center of present solid state research. Among TMO, cubic perovskites have the simplest crystal structure and, thus, may be viewed as a starting point for understanding the electronic properties of more complex systems. Typically, the 3$`d`$ states in those materials form comparatively narrow bands of width $`W23`$ eV which lead to strong Coulomb correlations between the electrons. Particularly simple are TMO with a 3$`d^1`$ configuration like Sr1-xCaxVO3 since they do not show a complicated electronic structure.

Intensive experimental investigations of spectral and transport properties of strongly correlated 3$`d^1`$ TMO started with the paper by Fujimori et al. AFSCVO  These authors observed a pronounced lower Hubbard band in the photoemission spectra (PES) measured at low photon energies ($`h\nu 120`$ eV) which cannot be explained by conventional band structure theory. A number of papers IHI95 ; Morikawa  subsequently addressed the spectral, transport and thermodynamic properties of Sr1-xCaxVO3, yielding contradictory results. While the thermodynamic and transport properties (Sommerfeld coefficient, resistivity, and paramagnetic susceptibility) do not change much with $`x`$, low photon energy (low-$`h\nu `$) PES show drastic differences for varying $`x`$. In fact, these spectroscopic data seemed to imply that Sr1-xCaxVO3 is on the verge of a Mott-Hubbard transition for $`x1`$. This is in concordance with the widespread theoretical expectation that the bandwidth $`W`$ should decrease strongly with increasing $`x`$ since the V-O-V angle is reduced from $`\theta =180^{}`$ for cubic SrVO3 Rey  to $`\theta =154.3^{}`$ (rotation angle), 171 (titling angle) MHJung  for orthorhombic CaVO3Chamberland

With the SPring-8 beamline we are now able to reexamine this puzzling discrepancy by high-resolution high-energy PES which, owing to a longer photoelectron mean free path $`\lambda `$ at $`h\nu 1000`$eV, is much more bulk-sensitive. Tanuma ; ASN ; RSI ; MaitiPRL ; MaitiEU  In this Letter, we show that these bulk-sensitive PES on fractured surfaces with $`100`$meV resolution are nearly independent of the Ca concentration $`x`$ and significantly different from high-energy PES at 0.5 eV resolution on scraped surfaces. MaitiEU  The good agreement with LDA+DMFT Metzner ; VIA97 ; AL2000 ; IAN2000 ; AIL2001 ; Savrasov ; Psik  spectra provides an important cross-validation: By means of a TMO prototype we show that the two cutting-edge methods employed here are indeed able to determine bulk spectra.

The high-resolution PES are summarized in Fig. 1ExpCond  In all spectra, the peak near $`E_F`$ and the broad peak centered at about $`1.6`$ eV correspond to the coherent (quasiparticle peak) and incoherent parts (lower Hubbard band), respectively. In contrast to low-$`h\nu `$ PES of Sr1-xCaxVO3 reported so far, including our own, the spectral variation with $`x`$ is noticeably smaller in the spectra measured at $`h\nu =`$ 900 eV as summarized in Fig. 1(b). As shown in Fig. 1(c), the coherent spectral weight decreases drastically with decreasing $`h\nu `$ for all compounds. One might consider that the relatively strong incoherent spectral weight of the $`1.6`$ eV peak at low $`h\nu `$ originates possibly from a large O $`2p`$ weight in this peak. However, both coherent and incoherent parts are drastically enhanced in a V $`2p3d`$ resonance photoemission (not shown), and are therefore assigned to the effective V $`3d`$ states around $`E_F`$O2p  Since the oxygen admixture to these predominantly $`3d`$ states is very minor, the relative V $`3d`$/O $`2p`$ cross section which also changes by only 8% from $`h\nu =`$ 275 and 900 eV Lindau  can be neglected. Also note that the effect of PES matrix elements is generally small for high-$`h\nu `$ angle-integrated spectra. Therefore, the monotonous increase of the coherent part with $`h\nu `$ is inevitably attributed to an increased bulk sensitivity.

We have obtained the bulk V $`3d`$O2p  PES, Figure 1(d), of Sr1-xCaxVO3 from the data at $`h\nu =`$ 900 and 275 eV by the following procedure HeII : (1) The mean free path $`\lambda `$ has been calculated as $``$17 and $``$7 Å at $`h\nu =`$ 900 and 275 eV. Tanuma  (2) The bulk weight $`R`$ ($`<1`$, depending on $`h\nu `$) should be determined as $`\mathrm{exp}(s/\lambda )`$ where $`s`$ is a ”surface thickness”. Therefore $`R`$s at 900 eV ($`R_{900}`$) and 275 eV ($`R_{275}`$) are related as $`R_{275}=R_{900}^{}{}_{}{}^{2.4}`$. (3) The observed V $`3d`$ PES intensity at $`h\nu =`$ 900 eV is represented as $`I_{900}(E)=I_B(E)R_{900}+I_S(E)(1R_{900})`$ while $`I_{275}(E)=I_B(E)R_{900}^{}{}_{}{}^{2.4}+I_S(E)(1R_{900}^{}{}_{}{}^{2.4})`$, where $`I_B(E)`$/$`I_S(E)`$ is the bulk/surface $`3d`$ PES and $`E`$ stands for the binding energy. (4) If $`s`$ is assumed to be 7.5 Å corresponding to about twice the V-O-V distance, IHI95  $`R_{900}`$ ($`R_{275}`$) is determined as $``$0.64 ($``$0.34). With this $`R_{900}`$ we obtain the $`3d`$ bulk PES $`I_B(E)`$ shown in Fig. 1(d), noting that $`I_B(E)`$ hardly changes ($`<15\%`$ at $`E_F`$) even when we assume $`s`$ to range from 5.4 to 11 Å. In contrast to the previous PES, the bulk $`3d`$ PES are almost equivalent among the three compounds, indicating that the V-O-V distortion does not much influence the occupied bulk $`3d`$ states in Sr1-xCaxVO3. This is consistent with the thermodynamic properties. It should be noticed that on scraped sample surfaces nearly $`x`$-independent bulk spectral functions could not be obtained. MaitiEU

Figure 2 shows the LDA density of states (DOS) for SrVO3 and CaVO3 which we obtained using the TBLMTO47 code of Andersen and coworkers LMTO  as well as a full-potential linearized augmented plane wave (FLAPW). The van-Hove-like peak of the latter can also be obtained by LMTO, Pavarini  leading to minor differences which are, however, not important for our LDA+DMFT calculations. Most importantly, the one-electron $`t_{2g}`$ bandwidth of CaVO3, defined as the energy interval where the DOS in Fig. 2 is non-zero, is found to be only 4% smaller than that of SrVO3 ($`W_{\text{CaVO}_3}`$ = 2.5 eV, $`W_{\text{SrVO}_3}`$ = 2.6 eV). The small reduction of the bandwidth well agrees with Pavarini  and will later turn out to be important.

In order to determine the genuine effect of the electron correlations on the bulk spectra, we first compare in Fig. 3(a) the measured spectrum with a V $`3d`$ partial DOS for SrVO3 obtained by using the FLAPW method. Here, the partial DOS broadened by the instrumental resolution (dashed curve) is normalized to the bulk $`3d`$ spectral function of SrVO3 by the integrated intensity from $`E_F`$ to $`2.6`$ eV. This comparison shows that the band-structure calculation does not reproduce the incoherent spectral weight at all and that the width of the observed coherent part is about 60% of the predicted value. Clearly, these deviations originate from electronic correlations which are not treated correctly in the band-structure calculation.

For a more quantitative analysis, we use LDA+DMFT as a non-perturbative approach, VIA97 ; Psik ; Pavarini  in which the local intra-orbital Coulomb repulsion $`U`$, the inter-orbital repulsion $`U^{}`$, and the exchange interaction $`J`$ are explicitly taken into account. We calculated these interaction strengths by means of the constrained LDA method OG89  for SrVO3. The resulting value of the averaged Coulomb interaction is $`\overline{U}`$ = 3.55 eV ($`\overline{U}=U^{}`$ for $`t_{2g}`$ orbitals Psik ; Zolfl ) and $`J`$ = 1.0 eV. $`U`$ is then fixed by the rotational invariance to $`U=U^{}+2J=`$ 5.55 eV. These calculated values are in agreement with previous calculations of vanadium systems Solovyev  and experiments. Makino  We use these values in our LDA+DMFT calculations for both SrVO3 and CaVO3 UJCa  and employ quantum Monte Carlo (QMC) QMC  simulations to solve the DMFT self-consistency equation. The resulting LDA+DMFT (QMC) spectra for SrVO3 and CaVO3 in Fig. 3(b) show genuine correlation effects, i.e., the formation of lower Hubbard bands at about $`2`$ eV and upper Hubbard bands at about $`2.5`$ eV, with well pronounced quasiparticle peaks (coherent part) near $`E_F`$. The results indicate that both SrVO3 and CaVO3 are strongly correlated metals, but not on the verge of a Mott-Hubbard metal-insulator transition. The 4% difference in the LDA bandwidth between the compounds is only reflected in some additional transfer of spectral weight from the coherent to incoherent parts. The many-particle spectra of the two systems are seen to be quite similar, especially in the occupied states, in agreement with Ref. Pavarini . Figure 3(c) shows that the effect of temperature on the spectrum is small for $`T`$ 700 K. The slight differences in the quasiparticle peaks between SrVO3 and CaVO3 lead to different effective masses, namely $`m^{}/m_0=2.1`$ for SrVO3 and 2.4 for CaVO3, which agree with $`m^{}/m_0=23`$ for SrVO3 and CaVO3 as obtained from thermodynamics IHI95  and de Haas-van Alphen experiments. IHI02

In Figure 3(d), we compare the experimental bulk V $`3d`$ O2p  PES to the LDA+DMFT spectra (300 K), which were multiplied with the Fermi function at the experimental temperature (20 K) and broadened with the PES resolution. The quasiparticle peaks in theory and experiment are seen to be in very good agreement. In particular, their height and width are almost identical for both SrVO3 and CaVO3. There is some difference in the positions of the incoherent part which may be partly due to uncertainties in the ab-initio calculations of $`\overline{U}`$ or in the oxygen subtraction procedure of Fig. 1. Indeed, some differences have to be expected since both the experimental determination of the bulk spectrum and the LDA+DMFT method involve approximations. Hence, our good overall agreement confirms that the new high-energy high-resolution PES indeed measures bulk spectra which can be calculated quite accurately by LDA+DMFT without free-parameters.

This work was supported by a Grant-in-Aid for COE Research from the Ministry of Education, Culture, Sports, Science and Technology (MEXT), Japan, by the Deutsche Forschungsgemeinschaft through SFB 484 and the Emmy-Noether program, the Russian Foundation for Basic Research through RFFI-01-02-17063 and RFFI-03-02-06126, the Ural Branch of the Russian Academy of Sciences for Young Scientists, Grant of the President of Russia MK-95.2003.02, and the National Science Foundation under Grant No. PHY99-07949. The photoemission at $`h\nu =`$ 275 and 900 eV was performed under the approval of the Japan Synchrotron Radiation Research Institute (1999B0076-NS-np). We thank T. Ushida, A. Shigemoto, T. Satonaka, T. Iwasaki, M. Okazaki, S. Kasai, A. Higashiya, K. Konoike, A. Yamasaki, A. Irizawa and the staff of SPring-8, especially T. Muro, T. Matsushita and T. Nakatani for supporting the experiments. We also thank A. Sandvik for making available his maximum entropy code and acknowledge valuable discussions with I. H. Inoue, A. Fujimori, R. Claessen, and M. Mayr.


text:
# Neutrino Oscillations, Fluctuations and Solar Magneto-gravity Waves

## 1 Introduction

Over the past decade a consistent picture of neutrino mixing has emerged, in which the oscillations responsible for both the solar  and atmospheric  neutrino deficits are consistent with the results of terrestrial neutrino-disappearance experiments . Experiments are now moving beyond the discovery phase and into a measurement phase during which the oscillation parameters are being determined with unprecedented precision.

Now that a coherent picture of neutrino properties seems to be in place, solar neutrinos can be used in the way the early investigators originally envisaged : as probes of the deep solar interior. Traditionally, the solar neutrino flux was believed to be exclusively sensitive to the properties of the solar medium deep inside the core since this is the region where the neutrinos are produced. The $`{}_{}{}^{8}B`$ neutrino flux predicted by solar models is very sensitive to the temperature in the solar core and the extremely weak neutrino interactions ensure that their observed energy spectrum is not affected by their passage through the solar medium (in the absence of oscillations). Indeed, we now know that the good agreement between solar-model predictions and the measured flux of neutrinos of all three species strongly constrains deviations of the core temperature from solar-model predictions. Better yet, these models are independently precisely tested by comparison with helioseismic measurements, again with good agreement between the models and observations .

### 1.1 Sensitivity to Solar Fluctuations

As was realized in the mid-nineties , the existence of resonant neutrino oscillations potentially makes the observed solar neutrino flux sensitive to other parts of the solar environment. In particular, this flux can be sensitive to fluctuations in the solar medium at the place where the neutrino resonance occurs. As a result, some of the properties of the solar medium at the neutrino resonance point may be inferred by comparing the measured neutrino energy spectrum with what is predicted by neutrino oscillations.

Of course such an inference of solar properties depends on having precise terrestrial observations of neutrino oscillation parameters, since these are required in order to cleanly predict the neutrino energy spectrum. Conversely, in the absence of terrestrial measurements, the precision of any determination of neutrino oscillation parameters using solar neutrinos can be degraded by the possible existence of solar fluctuations which affect the observed neutrino signal .

It is the purpose of the first part of this review to show that terrestrial  and solar  neutrino measurements have recently turned a corner inasmuch as they are now sufficiently accurate to allow the removal of solar uncertainties from the inference of neutrino oscillation properties. Conversely, the comparison of solar neutrino data with terrestrial observations now provides a clean window onto a new part of the deep solar interior.

### 1.2 Solar Magneto-gravity Waves

Aficionados have not been too alarmed by the necessity to assume the absence of solar fluctuations in order to infer neutrino properties, for several reasons. First, helioseismic measurements were known to constrain deviations of solar properties from Standard Solar Model predictions at better than the percent level. Second, the first studies of the implications for neutrino oscillations of radiative-zone helioseismic waves  showed that they were very unlikely to have observable effects. Third, no other known sources of fluctuations seemed to have the properties required to influence neutrino oscillations.

All three of these points have been re-examined in recent years and although it may yet turn out that solar fluctuations do not observably perturb neutrino oscillations, the possibilities are more promising than had been presumed earlier. Direct helioseismic bounds turn out to be insensitive to fluctuations whose size is as small as those to which neutrinos are sensitive  (which turn out to be those whose size is comparable to the neutrino oscillation length in matter: several hundreds of km). Furthermore, recent studies of magnetic fields deep inside the solar radiative zone  have identified potential fluctuations to which neutrinos might be sensitive after all (due to a resonance between Alfvén waves and helioseismic $`g`$-modes). It is the summary of this last observation which is the topic of the second half of this review.

We now turn to a more detailed description of these two topics.

## 2 Sensitivity to Solar Fluctuations

The standard description of MSW oscillations  amount to the use of a mean-field approximation for the solar medium. The corrections to this mean-field approximation are due to the fluctuations in the solar medium about this mean, and the leading interaction of neutrinos with these fluctuations are parameterized by the electron-density autocorrelation, $`\delta n_e(t)\delta n_e(t^{})`$, measured along the neutrino trajectory .

As fig. (1) shows, such fluctuations act to degrade the efficiency of neutrino oscillations. They can do so because successive neutrinos ‘see’ slightly different solar properties, and so in particular do not experience an equally adiabatic transition as they pass through the neutrino resonance region. The net effect is to degrade the effectiveness of the neutrino conversion because those neutrinos for which the transition is less adiabatic are more likely to survive as electron-type neutrinos. Since the criterion for the transition to be adiabatic depends on how quickly the electron distribution varies near resonance, fluctuations give observable effects for neutrinos if they occur at resonance with sufficient amplitude, and if their correlation length, $`L_0`$, is comparable to the local neutrino oscillation length, $`L_{\mathrm{osc}}100`$ km.

### 2.1 The Implications of KamLAND and SNO Salt

Ref.  has performed fits which are obtained using a global analysis of the solar data, including radiochemical experiments (Chlorine, Gallex-GNO and SAGE) as well as the latest SNO data in the form of 17 (day) + 17 (night) recoil energy bins (which include CC, ES and NC contributions, see  and the Super-Kamiokande spectra in the form of 44 bins .

The sensitivity of the solar neutrino data to fluctuations in the solar medium is summarized by figures (2) and (3). Fig. (2) is taken from ref. , and summarizes the sensitivity before the SNO salt measurements. Fig. (3) gives the same results after SNO salt. Comparing these figures shows the improvement in constraints due to the SNO salt data, and comparing the panels in each figure shows the the importance of a precise determination of the neutrino oscillation parameters for obtaining a constraint on the magnitude of fluctuations.

The importance of both the KamLAND and SNO salt measurements in these results is most easily seen from fig. (4), which compares the dependence of the fit’s $`\chi ^2`$ on the amplitude of the fluctuations for various data sets. This figure makes clear how the KamLAND results are largely responsible for localizing the best fit near zero fluctuation amplitude. This is as should be expected, since the evidence for the absence of fluctuations follows from the comparison of solar-neutrino observations with terrestrial measurements of neutrino oscillation properties.

Fig. (5) shows how the existence of solar fluctuations influences the determination of the neutrino oscillation parameters. The two panels of the figure contrast the precision of the fit with and without solar fluctuations. The left panel gives results subject to the usual prior assumption of no solar fluctuations, while the right panel leaves the amplitude of such fluctuations to be obtained from the fit. (When fluctuations are included, they are assumed to have the optimal correlation length $`L_0=100`$ km.) The lines indicate contours of fixed confidence level when the KamLAND data is not included, while the coloured regions give the same information when KamLAND is included.

The main conclusion which follows from this figure is that the precision with which the neutrino oscillations are known is now largely independent of whether a prior assumption is made about the existence of solar fluctuations. With the release of the SNO salt results the comparison of solar-neutrino with KamLAND data suffices to robustly determine the oscillation parameters independent of the assumed amplitude of solar fluctuations. The SNO salt data are crucial for reaching this conclusion, as is clear from fig. (6), which compares the right panel of fig. (5) with the same fit performed without using the SNO salt results.

We see that the SNO salt data, when combined with KamLAND results, for the first time places the determination of neutrino-oscillation parameters beyond the reach of sensitivity to prior assumptions concerning the existence of fluctuations in the solar radiative zone. Besides making more robust the determination of neutrino-oscillation parameters, this allows a much crisper determination of the kinds of solar fluctuations which can be entertained deep within the solar radiative zone. As we have seen, the resulting constraints apply to fluctuations whose spatial scales are of order 100 km, and so are complementary to those obtained from helioseismology, which are insensitive to fluctuations on such short scales.

## 3 Solar Magneto-Gravity Waves

Given the experimental sensitivity to solar fluctuations summarized above, the question remains as to whether a reasonable source of fluctuations in the solar medium might exist having sufficient amplitude and the correlation length required to observably affect solar neutrinos. Ordinary helioseismic waves are an obvious possibility, since they are known to exist deep within the sun. Physically, they might have an effect because the waves cause successive neutrinos to see different electronic density profiles and so have the effect of making the neutrino ‘jump probability’ at the neutrino resonance into a random variable.

The influence of helioseismic waves on neutrino oscillations was investigated in ref. , where it was found they are very unlikely to produce an observable effect. They do not for one of two reasons. For those waves which definitely have been observed ($`p`$-waves) the observed wave amplitude is much too small deep inside the solar radiative zone to produce observable effects. However, there are other helioseismic modes ($`g`$-modes) which have not yet been observed but which must exist within the solar radiative zone. It turns out that even if these modes are assumed to have an amplitude as large as a few percent, their wavelengths are too long to produce observable deviations from the predictions of neutrino oscillations in the absence of solar fluctuations.111A small number of potentially over-stable modes could evade both of these arguments, but only if they were to carry an inordinate amount of energy .

The remainder of this review summarizes the results of ref.  which has suggested a possible mechanism for obtaining observable fluctuations in the relevant part of the sun. In this picture fluctuations having the appropriate distance scale may arise if magnetic fields as large as of order 10 kG should exist deep in the solar core. Magnetic fields of this size would be consistent with current observational bounds . It is not yet known whether such modes could have sufficiently large amplitudes to allow observable effects on neutrino oscillations, but the main lesson to be learned from ref.  is that the detailed shape of helioseismic $`g`$-waves can be significantly changed by reasonable radiative-zone magnetic fields, and so their influence on neutrinos bears further study. The surprisingly strong sensitivity to magnetic fields turns out to be driven by the occurrence of level crossing between helioseismic $`g`$-modes and Alfvén waves.

Can 10 kG magnetic fields be present in the solar radiative zone? Very little is directly known about magnetic field strengths there – see  for early studies. A generally-applicable bound is due to Chandrasekhar, and states that the magnetic field energy must be less than the gravitational binding energy: $`B^2/8\pi \text{ }\stackrel{<}{}\text{ }GM_{}^2/R_{}^4`$, or $`B\text{ }\stackrel{<}{}\text{ }10^8`$ G. Stronger bounds are possible if one makes assumptions about the nature and origins of the solar magnetic field. For instance, if it is a relic of the primordial field of the collapsing gas cloud from which the sun formed , then it has been argued that central fields cannot exceed around 30 G . Similarly, the dynamo mechanism can only generate a global field in the radiative zone of a newly-born Sun with amplitude below 1 G . Even stronger limits, $`B\text{ }\stackrel{<}{}\text{ }10^6`$ G apply  if the solar core is rapidly rotating, as is sometimes proposed. On the other hand, it has recently been argued  that fields up to 7 MG could persist in the radiative zone for billions of years and are consistent with current observational bounds. Other authors  have recently entertained radiative-zone fields as large as 30 MG. Since the initial origin and current nature of the central magnetic field is unclear, we consider as admissible any magnetic field smaller than of order 10 MG.

### 3.1 Magneto-gravity waves

In this section we briefly review how magnetic fields influence the equations of hydrostatic equilibrium on which helioseismic analyses are based, and the approximations used to solve them.

The first relevant equations express conservation of mass:

$$\frac{d\rho }{dt}+\rho u=0,$$
(1)
where $`d/dt=/t+\stackrel{}{v}`$ is the usual convective derivative, with $`\stackrel{}{v}`$ representing the fluid velocity, and $`\rho `$ and $`p`$ are the fluid’s mass density and pressure. The variable $`u=\mathrm{div}\stackrel{}{v}`$ need not vanish if the fluid is compressible.

The second equation of interest expresses energy conservation:

$$\frac{dp}{dt}\gamma \frac{p}{\rho }\frac{d\rho }{dt}=(\gamma 1)Q,$$
(2)
where $`\gamma =c_p/c_V`$ is given by the ratio of heat capacities, and $`Q`$ is the sum of all energy density sources and losses, such as heat conductivity, viscosity and and ohmic dissipation.

The magnetic fields enter more directly into the Euler equation (conservation of momentum), which is of the form

$$\rho \frac{d\stackrel{}{v}}{dt}=p+\rho \stackrel{}{g}+\frac{1}{4\pi }\left[\mathrm{curl}\stackrel{}{B}\times \stackrel{}{B}\right],$$
(3)
with the local force of gravity given by $`\rho \stackrel{}{g}`$. The local acceleration due to gravity is related to the Newtonian potential by $`\stackrel{}{g}=\varphi `$, with $`\varphi `$ given by the Poisson equation $`^2\varphi =4\pi G\rho `$. As usual, $`G`$ here denotes Newton’s constant. The last term of equation (3) expresses the contribution of the Lorentz force to the local momentum budget.

Finally, the system is completed by Faraday’s equation,

$$\frac{\stackrel{}{B}}{t}=\mathrm{curl}[\stackrel{}{v}\times \stackrel{}{B}]+\frac{c^2}{4\pi \sigma _{cond}}^2\stackrel{}{B},$$
(4)
that governs the time evolution of the magnetic field. Here $`\sigma _{cond}`$ represents the fluid’s conductivity. In what follows we take the plasma to be an ideal conductor, meaning we take $`\sigma _{cond}`$ to be large enough to neglect the last term in this last equation.

In order to find solutions to these equations a number of assumptions are required. We list them all here for ease of reference. In deriving the differential equations to be solved, we adopt the following approximations.

1. In principle, electromagnetic processes enter into Eq. (2) through their contributions to $`Q`$, but for ideal MHD we neglect both the heat conductivity and viscosity contributions to energy losses, as well as the ohmic dissipation, $`Q=j^2/\sigma _{cond}`$, where $`\sigma _{cond}`$ is the electrical conductivity.
2. We linearize the equations about a static background configuration, i.e. a background configuration which is time independent and for which the background fluid velocity vanishes, $`\stackrel{}{v}_0=0`$.
3. We assume the thermodynamic quantity, $`\gamma =c_p/c_V`$, has vanishing derivative for a Lagrangian fluid element, $`d\gamma /dt=\gamma /t+𝐯\gamma =0`$. This would be true, in particular, for a polytropic gas.
4. We adopt the Cowling approximation, which amounts to the neglect of perturbations of the gravitational potential, (i.e.: $`\varphi ^{}=0`$).
5. We assume the fluctuations to be adiabatic, with the contributions of fluctuations to the heat source vanishing: $`Q^{}=0`$ (this is satisfied with good accuracy by high frequency oscillations).
6. We assume no background electric or displacement currents, so the background magnetic field satisfies $`\mathrm{curl}\stackrel{}{B}_0=0`$.

Finally, we are able to solve the resulting equations analytically if we make two further assumptions.

1. We assume a rectangular geometry, with background quantities varying along the $`z`$ direction (which implies the local gravitational acceleration, $`𝐠`$, is directed along the $`z`$ axis). We also take a constant, uniform background magnetic field, $`𝐁_0`$, pointing along the $`x`$ axis.
2. The background mass-density profile is assumed to be exponential, $`\rho _0=\rho _c\mathrm{exp}[z/H]`$, for constant $`\rho _c`$ and $`H`$. The conditions of hydrostatic equilibrium for the background then determine the profiles of thermodynamic quantities, and in particular imply $`\gamma `$ is a constant.

This last approximation applies to very good approximation for real mass-density profiles obtained by standard Solar models, provided we identify the $`z`$ direction with the radial direction, and focus our attention to deep within the radiative zone. The constancy of $`\gamma `$ in this region is also expected since the highly-ionized plasma satisfies an ideal gas equation of state to good approximation. The rectangular geometry provides a reasonable approximation so long as we do not examine too close to the solar centre. What is important about our choice for $`𝐁_0`$ is that it is slowly varying in the region of interest, and it is perpendicular to both $`𝐠`$ and all background gradients, $`\rho _0`$, $`p_0`$, etc.

### 3.2 The eigenvalue problem

Using the above assumptions we can express all oscillating quantities in terms of any one component, which we take to be the profile $`b_z(z)`$. This quantity turns out to satisfy a linear ordinary second-order differential equation whose solution determines all of the other variables through simple expressions . The relevant ordinary differential equation is

$$\left(1\frac{k_x^2v_A^2}{\omega ^2}\right)\frac{d^2b_z(z)}{dz^2}\frac{N^2}{g}\frac{db_z(z)}{dz}+k_{}^2\left(\frac{k_x^2v_A^2}{\omega ^2}1+\frac{N^2}{\omega ^2}\right)b_z(z)=0,$$
(5)
where $`N`$ denotes the Brunt-Väisälä frequency, as defined by

$$N^2(z)=g(z)\left(\frac{1}{\gamma p_0}\frac{dp_0}{dz}\frac{1}{\rho _0}\frac{d\rho _0}{dz}\right),$$
(6)
and $`v_A=B_0/\left(4\pi \rho _0\right)^{1/2}`$ defines the Alfvén speed. $`k_x`$ and $`k_y`$ here denote the wave-number in the transverse directions. Equation (5) contains all of the information required about the system’s normal modes. Notice that its derivation allows the quantities $`g`$, $`\gamma `$, $`\rho _0`$ and $`p_0`$ to depend generally on $`z`$, subject only to the conditions of hydrostatic equilibrium and vanishing Lagrangian derivative, $`d\gamma /dt=0`$.

The boundary conditions to be imposed are as follows. At the solar center ($`z=0`$) we use the boundary condition which would have arisen as a smoothness condition in cylindrical coordinates:

$$b_z(0)=k_xv_z(0)/\omega =0.$$
(7)
The other boundary condition is imposed at the top of the radiative zone, $`z_{}0.7R_{}`$. In this region the two solutions to the differential equation behave as $`e^{\pm k_{}z}`$ (where $`k_{}^2=k_x^2+k_y^2`$) and so either fall or grow exponentially as functions of $`z`$. We take the absence of the growing behaviour to be our boundary condition for $`z=z_{}`$.

It is instructive to examine the qualitative properties of the eigenvalue problem we have obtained. In the limit $`B_00`$ (and so $`v_A0`$), Eq. (5) reduces to the standard evolution equation for ‘pure’ helioseismic $`g`$-modes, which is usually expressed in terms of the variable $`v_z`$ (which is related to $`b_z`$ through the relation $`v_z=\omega b_z/k_x`$). Since more detailed analyses – see Fig. 7 – show that $`N`$ rises from zero at the solar center, remains approximately constant $`NN_0`$, through the radiative zone, and then falls to zero again at the bottom of the convective zone, these $`g`$-modes can be thought of as the eigen-modes of oscillations inside the cavity in between the two regions where $`N`$ goes through zero. For a given wave frequency, $`\omega <N_0`$, the turning points of this cavity are given by the condition $`\omega N`$. For smaller $`\omega `$ the lower turning point gets closer to the solar center, $`z0`$, and the upper one gets slightly closer to the bottom of the convective zone (CZ).

Conversely, if gravity is turned off ($`N0`$), then Eq. (5) describes Alfvén waves, which oscillate with frequency $`\omega =k_xv_A`$ and propagate along the magnetic field lines. Notice that since $`v_A\rho ^{1/2}`$ this frequency grows with $`z`$, since the density of the medium falls.

Keeping both magnetic and gravitational fields introduces qualitatively new behavior, as may be seen mathematically because equation (5) acquires a new singular point which occurs when the coefficient of the second-derivative term vanishes. Since this singularity appears at the Alfvén frequency,

$$\omega =k_xv_A$$
(8)
it can be interpreted as being due to a resonance between the $`g`$-modes and Alfvén waves. This resonance turns out to occur at a particular radius because the Alfvén frequency varies with radius while the $`g`$-mode frequency is independent of radius (see Fig. 7). The resonance occurs where the growing Alfvén mode frequency crosses the frequency of one of the $`g`$-modes, and the resulting waveforms would be expected to vary strongly at these points. The resonance can occur inside the radiative zone if the Alfvén frequency climbs high enough to cross a $`g`$-mode frequency before reaching the top of the radiative zone. Since $`v_AB_0`$, whether this occurs or not depends on the field value, $`B_0`$.

The qualitative behaviour of the solutions can be seen from equation (5), and are oscillatory (or exponentially growing or damped) if the coefficient of the last term of this equation is positive (or negative). In the absence of magnetic fields, this leads to the onset of damping when $`\omega >N(z)`$. Since $`N(z)`$ vanishes at the top of the radiative zone, all modes necessarily become unstable there, and the corresponding instability gives rise to the convection which defines the convective zone. This is the usual picture of helioseismic $`g`$ modes. The presence of magnetic fields makes the coefficient of $`b_z`$ in eq. (5) go negative (for some modes) at smaller values of $`z`$, indicating the onset of an instability inside of the radiative zone.

The magnetic field leads to several new effects. First, the shortening of the cavity due to the presence of magnetic field causes the eigen-frequencies to depend on the magnetic field value. Second, there can be energy transfer between $`g`$-modes and Alfvén waves, within the narrow singular resonance layer, leading to the corresponding eigen-frequencies acquiring imaginary parts. Third, the nearer the upper boundary of MHD cavity is to the solar center, the stronger the $`g`$-modes are confined to the solar core.

More details of these waveforms can be obtained from a WKB-type analysis of the master Eq. (5). Near the solar center where $`N^20`$, if $`(k_xv_A)^2/\omega ^21`$ one obtains the exponential solution $`b_ze^{\pm k_{}z}`$. (Which combination of these solutions appears is fixed by boundary conditions, e.g. if $`b_z(0)=0`$ at the solar center we have $`b_z\mathrm{sinh}(k_{}z)`$.) One finds similar behaviour for the solution above the singular resonant layer, $`z>z_r`$, up to the top of the radiative zone. This exponential growth happens because of the exponential decrease of density (and so exponential increase in $`v_A`$) with $`z`$. The requirement for complex frequencies arises from the demand that the solution remain regular in the narrow Alfvén resonance layer.

Figures 8 and 9 present our numerical solution of the eigenvalue spectrum. Figure 8 plots $`v_{ph}(n)/c_s=\omega _1/(c_sk_x)`$ vs mode number $`n`$ for various magnetic fields, $`B_0`$. Here $`\omega _1`$ is the real part of the eigenfrequency $`\omega `$, which in general is complex. (More about this later.) Figure 9 similarly plots the imaginary part of omega, Im $`\omega =\omega _1d`$ against mode number for the same magnetic fields. In both plots the parameter $`\alpha =k_x/k_{}`$ is taken to be unity, where as before $`k_x`$ and $`k_{}^2=k_x^2+k_y^2`$ are the mode’s wave-numbers in the transverse directions. In both figures a dashed line is used if the resonance of interest occurs above the top of the radiative zone (and so outside the domain of many of our approximations).

Both of the figures 8 and 9 show two branches of solutions up to a maximum mode number, as expected. Their dependence on $`n`$ can also be understood analytically from following approximate expressions

$`{\displaystyle \frac{\omega _1}{N}}={\displaystyle \frac{2k_{}H}{\pi n}}\mathrm{ln}\left({\displaystyle \frac{4NH}{\alpha v_{Ac}}}{\displaystyle \frac{1}{\pi n}}\right),`$
$`d={\displaystyle \frac{\mathrm{ln}\left|\mathrm{tan}(\pi /\gamma )\right|}{2\pi n}}.`$ (9)
These modes correspond to the branches of the figures 8 and 9 for which both $`v_{ph}/c_s`$ and $`d`$ fall with $`n`$. Analytic expressions are also possible for the other branch , and the spectrum in this case is

$`\omega _1={\displaystyle \frac{\alpha K_{}v_{Ac}}{H}},`$
$`d=\left({\displaystyle \frac{\alpha v_{Ac}}{4NH}}\right)^2\pi n\mathrm{ln}\left|\mathrm{tan}{\displaystyle \frac{\pi }{\gamma }}\right|.`$ (10)
Note that for this branch $`v_{ph}(n)/c_s`$ is independent of the mode number, $`n`$, and $`d`$ grows with $`n`$, as is also seen in the figure.

In addition to requiring the resonance to occur in the radiative zone (the solid line in the figures), the validity of our approximations also demand the frequency not to be smaller than $`10^5`$ s-1 due to our neglect of the 27-day period solar rotation. Inspection of Fig. 8 shows that these two conditions are consistent with one another for a reasonably large range of modes only for magnetic fields larger than a kG or so.

The resonance alluded to above appears causes two distinctive features to appear in these solutions. First, the eigenfrequencies are complex, implying both damped and exponentially-growing modes. Second, the eigenmodes are not smooth as functions of $`z`$ about the singular resonant point, $`z=z_r(n,k_x,k_y)`$, whose position depends on the quantum numbers of the mode in question.

The necessity for complex frequencies imply the mode functions grow exponentially in amplitude with time. This signals an instability in the physics which pumps energy into these resonances, and this section aims to discuss the nature of this instability, and the physical interpretation to be assigned to the imaginary part, $`d`$. Exponentially-growing instabilities within an approximate (e.g. linearized) analysis reflect the system’s propensity to leave the small-field regime, on which the validity of the approximate analysis relies. The question becomes: where does the instability lead, and what previously-neglected effects stabilize the runaway behaviour?

In the present instance the normal modes are strongly peaked near the resonant radii, and the energy flow near these radii is directed along the resonance plane, much as would be true for a pure Alfvén wave. Since helioseismic waves are likely generated by the turbulence at the bottom of the convective zone, it is natural to imagine starting the system with a regular helioseismic $`g`$-mode and asking how it evolves. In this case the resonance allows the energy in this mode to be funnelled into the Alfvén wave, and so to be channelled along the magnetic field lines away from the solar equatorial plane. The imaginary part of the frequency, $`\mathrm{\Gamma }=\text{Im}\omega =\omega _1d`$, describes the rate at which the Alfvén mode is excited in this process.

Once excited, the mode amplitudes near the resonances grow until the energy in them becomes dissipated by effects which are not captured by the approximate discussion we present here. The rate for this dissipation must grow as the mode amplitude grows, until it equals the production rate, $`\mathrm{\Gamma }`$, at which point a steady state develops and the mode growth is stabilized. Because the resonant mode grows until its damping rate equals its production rate, it suffices to know the mode’s growth rate in order to determine the overall on-resonance amplitude. If the mode stabilizes once it is large enough to require a nonlinear analysis, then the final production and dissipation rates may be very different from the linearized expressions derived above. If, on the other hand, the modes saturate at comparatively small amplitudes by dissipating energy into non-hydrodynamical modes (such as by Landau damping), then it can happen that the stabilized mode amplitude is not large enough to significantly change the linearized prediction for its production rate.

Fig. 10 plots the position, $`z_r(n)`$, predicted for the case of longitudinal wave propagation ($`k_y=0`$) and for magnetic fields in the range $`B_0=`$10G–1MG. (The same result for an obliquely-propagating wave with $`k_y0`$ is produced by a larger value for $`B_0`$.) This plot is the analog of Figs. 5a,b in . Knowing the position of the resonances also permits us to determine the distance between them. This quantity is relevant to the propagation of particles like neutrinos through the resonant waves. The spacing is:

$$z_r(n+1)z_r(n)=\left(\frac{k_x}{k_{}}\right)\frac{\pi Hv_{Ac}\gamma }{c_s\sqrt{\gamma 1}}e^{z_r(n)/2H}.$$
(11)
This dependence of this quantity on mode number, $`n`$, is shown in Fig. 11 for longitudinal wave propagation, $`k_y=0`$, and for different magnetic field values. From this figure we see, in the region $`z_r\text{ }\stackrel{>}{}\text{ }0.3R_{}`$, that the distance between layer positions grows with magnetic field and with distance from the solar center. Both of these features were already seen in preliminary WKB calculations of the resonances (see Figs. 7a, 7b of ref. ), and our present numerics also confirm that the spacing is approximately proportional to $`B_0`$ for $`z\text{ }\stackrel{>}{}\text{ }0.3R_{}`$, that was seen earlier.

An approximate expression for the differenece between adjacent resonance positions is given by

$$z_r(n+1)z_r(n)\frac{2H}{n},$$
(12)
where the inverse mode number $`n^1`$ is proportional to the background magnetic field $`B_0`$.

Numerically, it is noteworthy that the spacing between resonances can be hundreds of kilometers. This is significant because this is close to the resonant oscillation length, $`l_{osc}^{res}`$, for $`E\mathrm{MeV}`$ neutrinos, if – as now seems quite likely – resonant LMA oscillations are responsible for explaining the solar neutrino problem since

$$l_{osc}^{res}=\frac{250\mathrm{km}(\mathrm{E}/\mathrm{MeV})}{\mathrm{\Delta }m_5^2\mathrm{sin}2\theta },$$
(13)
where $`\mathrm{\Delta }m_5^2=\mathrm{\Delta }m^2/10^5\text{eV}^2`$. Repeatedly perturbing neutrinos over distance scales comparable to $`l_{osc}^{res}`$ has long been known to be a prerequisite for disturbing the standard MSW picture of oscillations in the solar medium. This raises the possibility – recently explored in more detail in ref.  – that $`g`$-mode/Alfvén resonances can alter neutrino propagation. If so, the observation of resonant oscillations of solar neutrinos may provide some direct information about the properties of the MG waves we discuss here.

## 4 Summary

Within the approximations given it appears that sufficiently large magnetic fields can cause significant changes to the profiles of helioseismic $`g`$-waves, while not appreciably perturbing helioseismic $`p`$-waves. The comparatively large $`g`$-wave effects arise because of a resonance which occurs between the $`g`$-modes and magnetic Alfvén waves in the solar radiative zone. Although the radiative-zone magnetic fields required to produce observable effects are larger than are often considered – more than a few kG – they are not directly ruled out by any observations.

Although the density profiles at their maxima could have amplitudes as large as a few percent or more on resonance, we do not believe the corresponding radiative-zone magnetic fields can yet be ruled out by comparison with helioseismic data, since the density excursions are sufficiently narrow (hundreds of kilometers) as to evade the assumptions which underlie standard helioseismic analyses.

For magnetic fields in the 10 kG range, the best hopes for detection of the resonant waves may be through their influence on neutrino propagation. This influence essentially arises because the presence of strong density variations affects the solar neutrino survival probability, with a corresponding change in the resulting solar neutrino fluxes. As described above, and shown in ref. , the measurement of neutrino properties at KamLAND provides new information about fluctuations in the solar environment on correlation length scales close to 100 km, to which standard helioseismic constraints are largely insensitive. We have already seen how the determination of neutrino oscillation parameters from a combined fit of KamLAND and solar data depends strongly on the magnitude of solar density fluctuations.

Since the resonances rely on the condition that $`\stackrel{}{B}\rho `$, there are several magnetic-field geometries to which our analysis might apply, and it is instructive to consider two illustrative examples to see what kinds of observable effects might be possible. Suppose first that, in spherical coordinates $`(r,\theta ,\varphi )`$, we imagine $`B_r0`$ but $`B_\theta 0`$. Then the field is always perpendicular to a radial density gradient and the resonance we find might be expected to arise in all directions as one comes away from the solar center. In this case the solar $`g`$-modes would tend to be trapped behind the resonance, and so are kept away from the solar surface even more strongly than is normally expected. This would make the prospects for their eventual detection very poor.

Alternatively, if the magnetic field has more of a dipole form it might be imagined that the condition $`\stackrel{}{B}\rho `$ only holds near the solar equatorial plane, and not near the solar poles. In this case a more detailed calculation is needed in order to determine the resulting wave form. This kind of geometry could have interesting consequences for the solar neutrino signal, because in this case the deviation from standard MSW analyses only arises for neutrinos which travel sufficiently close to the solar equatorial plane. Given the roughly 7-degree inclination of the Earth’s orbit relative to the plane of the solar equator, there is a possibility of observing a seasonal dependence in the observed solar neutrino flux. Since the presence of the MG resonance tends to decrease the MSW effect, the prediction would be that the observed rate of solar electron-neutrino events is maximized when the Earth is closest to the solar equatorial plane (December and June) and is minimized when furthest from this plane (March and September).

## Acknowledgements

This work was supported by Spanish grants BFM2002-00345, by the European Commission RTN network HPRN-CT-2000-00148, by the European Science Foundation network grant N. 86, by Iberdrola Foundation (VBS) and by INTAS grant YSF 2001/2-148 and CSIC-RAS agreement (TIR). C.B.’s research is supported by grants from NSERC (Canada), FCAR (Quebec) and McGill University. M.M. is supported by contract HPMF-CT-2000-01008. VBS, NSD and TIR were partially supported by the RFBR grant 00-02-16271. C.B. would like to thank the University of Valencia for its hospitality during part of this work.


text:
# On the C and N Abundances of 47 Tucanae Main Sequence StarsBased on observations made with ESO telescopes at the Paranal Observatory under program ID 67.D-0153 and the MPG/ESO 2.2 m telescope at La Silla Observatory

## 1 Introduction

It has been more than three decades since the first reports of star-to-star differences in light element abundances among members of Galactic globular clusters (Osborn, 1971). That these intra-cluster inhomogeneities are confined to C, N, O, and often Na, Al, and Mg, suggests an origin involving element production via the CNO-cycle and associated proton-addition reactions. Two likely sites for this nucleosynthesis exist: One, within the presently evolving cluster stars during red giant branch (RGB) ascent, where conditions above the H-burning shell may allow material exposed to partial CN(O)-processing to be mixed to the stellar surface (i.e., a deep dredge-up mechanism). A second possible site lies in the asymptotic giant branch (AGB) phase of more massive cluster members (now white dwarfs), whose ejecta were incorporated into the present stars.

Unfortunately, the signatures of both sources appear remarkably similar and observing the presence of CN(O)-cycle material in the photosphere of an evolved cluster RGB star is not a sufficient discriminant as to its origin. This has been made abundantly clear by the detection of O-Na and Mg-Al anti-correlations among NGC 6752 main-sequence turn-off (MSTO) and subgiant branch (SGB) stars by Gratton et al. (2001). The two possibilities may however be disentangled by exploring a span in luminosity (evolutionary state) from the RGB tip to the main-sequence (MS). Abundance changes with position on the RGB signal the operation of internal processes, while composition patterns present among MS stars are in place before RGB ascent and the onset of deep mixing.

As telescope apertures and instrumentation have improved, such tests have become increasingly feasible and data now exist on abundance variations among low luminosity stars of a number of clusters (47 Tuc, M71, M5, NGC 6752, and M13, see Cannon et al., 1998; Cohen, 1999b; Gratton et al., 2001; Briley, Cohen, & Stetson, 2002; Cohen, Briley, & Stetson, 2002; Cannon, Croke, da Costa, & Norris, 2002; Harbeck, Smith, & Grebel, 2003, for recent studies). Perhaps not so surprisingly, we find ourselves now facing a wealth of evidence that suggests not one origin or the other but rather both. The progressive depletions in \[C/Fe\] with decreasing MV seen among the stars of the more metal-poor clusters (e.g. M92, M15, NGC 6397, see Bellman et al., 2001, and references therein), the drop in Li abundances at the luminosity function “bump” observed by Grundahl et al. (2002) in NGC 6752, and the change in isotopic C abundances with luminosity in NGC 6528 and M4 found by Shetrone (2003) provide dramatic evidence of a deep dredge-up mechanism. Yet also, significant star-to-star differences in C and N have been found among the faintest stars of every cluster studied to the level of the SGB or below. Most recently, Briley et al. (2002) documented substantial \[C/Fe\] differences in a sample of M13 SGB stars which appeared to be subsequently modified during RGB ascent. Thus it would seem that in order to fully quantify the effects of deep mixing one must know the compositions of the stars before they reach the RGB.

The process(es) which set the compositions of the cluster pre-RGB stars presumably operated early in the cluster histories. That large fractions of the members of every cluster studied to date show significant differences of at least C and N imply this was not the result of some unique condition among one or two clusters and it involved a non-trivial portion of each cluster’s mass. As has been discussed by Briley, Smith, & Claver (2001) – the process appears not to be a simple accretion event. It is evident then that a complete understanding of the origin of abundance inhomogeneities in globular clusters must take account of the pattern of light element abundances among the main sequence stars in these objects. With this in mind, we present a C and N abundance analysis of MS stars in 47 Tucanae, one of the most favorable globular clusters in the Milky Way for the study of main sequence abundances.

We base our analysis on the data set obtained by Harbeck et al. (2003), who detected significant CN and CH variations as faint as 2.5 mag below the MSTO of 47 Tuc. The occurrence of CN variations among these low-mass stars favors accretion or primordial scenarios, since the capability of such stars to maintain interior CNO-cycle burning is reduced by a factor of at least ten compared to the main-sequence turn-off; thus such stars are not expected to be able to modify their surface C and N abundances.

## 2 Observations, Indices, and Membership

The sample selection, observations, and reduction of the 47 Tuc MS and MSTO spectra used in this analysis are described in Harbeck et al. (2003). To summarize, 115 program stars were selected from MPG/ESO 2.2 m images of 47 Tuc to span a range in luminosity from the MSTO to $`2.5`$ mag fainter. Photometry of the stars was also extracted from these images. The spectra were obtained at the ESO VLT using multi-slit masks with the Focal Reducer and Low Dispersion Spectrograph 2 and cover from 3800 to 5500 Å at R = 815. The locations of the stars on the 47 Tuc color-magnitude diagram (CMD) are shown in Figure 1.

Membership was assessed in Harbeck et al. (2003) based on radial velocity (unfortunately only $`18.7`$ km s-1 for 47 Tuc, Harris, 1996) and position in the CMD (within 0.025 mag of the cluster fiducial). As a further test of membership, we have measured an MgH sensitive index (MgH, with a feature bandpass from 5160 to 5195 Å, and continuum bandpasses from 5130 to 5160 and 5195 to 5225 Å) which is plotted in Figure 2 (see also Table 1). Of the stars in the sample, an additional 11 have MgH band strengths inconsistent with the other members of the cluster. We will presume these stars to also be non-members of 47 Tuc and they will not be included further.

Indices sensitive to absorption by CN at 3883 Å and CH at 4300 Å were measured from these spectra and are presented in Harbeck et al. (2003). However, here we use the I(CH) index as defined in Cohen (1999a) rather than CH(4300) to aid in comparison with the low luminosity stars of M71 (see Cohen, 1999b; Briley & Cohen, 2001). The CN (S(3839)) and CH (I(CH)) indices of the stars in the sample which are presumed members and further have measurable CN, CH indices are given in Table 1 and plotted in Figure 3, where the CN-bimodality and CN-CH anti-correlation noted by Harbeck et al. (2003) is quite evident. The error bars shown are at the one-sigma level as determined from Poisson statistics in the index bandpasses.

As noted in Cohen et al. (2002) single-sided indices such as S(3839) can suffer from zero-point offsets due to differences in slopes caused by flat fielding, instrumental response, etc. The observed spectra were therefore roughly flux calibrated before the indices were measured by dividing by a sensitivity function derived from a spectrum of LTT 377 taken during the same run. Thus the observed and model (see below) indices are assumed to share essentially the same zero point which should put the abundances on a scale consistent with that used in previous analyses of other stars in this cluster.

## 3 Model Atmospheres, Synthetic Spectra, and Resulting Abundances

The approach we have taken is identical to that of Cohen et al. (2002): establishing a $`T_{eff}`$/log($`g`$) versus luminosity relation for a grid of model atmospheres based on an isochrone fit, then adjusting the C and N abundances of corresponding synthetic spectra to simultaneously fit the observed CN and CH band strengths. The isochrone used was that of Bergbusch & VandenBerg (2001) with an alpha-element enhancement of 0.3 dex and an overall metallicity of $`0.83`$. An observed distance modulus of 13.33 and a reddening of 0.04 were assumed following Grundahl, Stetson, & Andersen (2002). The fit of the isochrone to the region of the CMD of interest is shown in Figure 1.

For each of the program stars, a $`T_{eff}`$ and log($`g`$) was derived by transposing its V magnitude onto the isochrone in Figure 1 and interpolating between grid points. These values were then used as inputs to MARCS (Gustafsson et al., 1975) for the generation of model atmospheres. The same overall metallicity (\[Fe/H\] = $`0.8`$) as was used in previous studies of more luminous 47 Tuc stars (see Briley, 1997) and of M71 MS/MSTO stars (Briley & Cohen, 2001) has been employed. Synthetic spectra were computed from each model atmosphere for various C and N abundances using the SSG program (Bell, Paltoglou, & Tripicco, 1994, and references therein) and the line-list of Tripicco & Bell (1995) (see Briley & Cohen, 2001, for details). A microturbulent velocity of 2 km s-1 was assumed throughout the sample, and an \[O/Fe\] of +0.45 was used (see Cannon et al., 1998, for a discussion of the minor role O plays in setting the CN/CH band strengths of 47 Tuc MS/MSTO stars). The synthetic spectra, calculated at steps of 0.05 Å, were convolved with a Gaussian of width 4.5 Å to match the resolution of the observed spectra and S(3839) and I(CH) indices measured from the result. Sample synthetic indices appropriate to two sets of C and N abundances over the luminosity range of the present sample are shown in Figure 3. Errors in the resulting abundances were calculated in a similar fashion using indices adjusted by the $`\pm 1\sigma `$ uncertainties from Poisson statistics.

With two exceptions (stars 4554 and 6622 – the brightest and faintest star in the sample), all of the program star models yielded C and N abundances which simultaneously fit the observed indices to 0.005 or less (i.e., to within 0.01 dex in \[C/Fe\] and \[N/Fe\], which is the limiting precision of the abundance inputs to the synthesis program). These values are shown in Table 1. Under the present method of extracting abundances, it is the V photometry of a star which is used to assign its $`T_{eff}`$ and log($`g`$) through fitting to the isochrone. Thus, errors in the photometry as well as the assumed distance modulus will lead to differences in the derived $`T_{eff}`$, log($`g`$), and resulting abundances. Fortunately, with the exception of the MSTO region, the change in $`T_{eff}`$, log($`g`$) with V is quite gradual. As a test of this, we have recomputed the results for three pairs of CN-weak/strong stars (one pair near the MSTO, one pair towards the luminous end of the MS sample, and one pair near the faint end) assuming a distance modulus of 13.23 rather than 13.33 (this is also equivalent to the change expected from a photometric error of 0.1 in the brightness of an individual star). These results are presented in Table 2 where it can be seen that indeed, the MS results are largely insensitive to errors in photometry/distance.

The resulting abundances (sans stars 4554 and 6622) are plotted in Figures 3 and 4, where a clear C versus N anti-correlation may be seen, especially in the case of the CN-strong stars. As is cautioned by Cohen et al. (2002), the method employed here can naturally lead to such anti-correlations in the results (e.g., underestimating \[C/Fe\] will naturally require larger \[N/Fe\] abundances to reproduce a necessary CN band strength). However, the indices of Figure 3 lead us to conclude that the anti-correlation of Figure 4 is likely quite real.

## 4 Discussion

We begin with a very brief summary of what is currently known about abundance variations among low luminosity (MSTO/MS) stars in the Galactic globular clusters, and in particular for 47 Tuc. We consider here only such stars as they presumably reflect the compositions of cluster members before the action of any deep mixing mechanism:

1) Star-to-star abundance differences have been observed among the low luminosity stars of 47 Tuc, NGC 6752, M5, M71, and M13 (see above references). In this sense, the intra-cluster light element abundance differences between low luminosity stars appear ubiquitous among the Galactic globular clusters. That similar patterns are not observed among metal-poor field stars (e.g., Carretta, Gratton, & Sneden, 2000) implies that some aspect of the cluster environment plays a role.

2) The pattern of the MS/MSTO variations exhibits the same correlations and anti-correlations expected from proton-capture nucleosynthesis, e.g., the C-N anti-correlation (see Figure 4), the N-Na correlation of 47 Tuc (Briley et al., 1996), and the O-Na, Mg-Al anti-correlations of NGC 6752 and M71 (Gratton et al., 2001; Ramírez & Cohen, 2002). It is widely accepted that the present-day cluster MS stars are incapable of reaching the $`10^8`$ K temperatures required for these reactions, which necessitates a source external to the present-day stars – this consideration has lead to several so-called “primordial enrichment scenarios.” The most often suggested site is intermediate-mass AGB stars undergoing hot bottom burning and third dredge-up (see Ventura et al., 2001), although difficulties such as the establishment of an O-Na anti-correlation remain (see for example Denissenkov & Herwig, 2003). Primordial scenarios include a number of possibilities as discussed in Cannon et al. (1998), among them: a pollution/accretion process where already formed low-mass cluster stars acquire AGB ejecta at a later date, or an extended period of star formation in which some stars formed from gas already mixed with AGB ejecta. Of course any viable scenario of this sort must not only explain the star-to-star variations in light elements but also the remarkable consistency of heavier (e.g. Fe-peak and $`\alpha `$-capture) elements within the clusters, which can be a problem for some primordial theories (such as supernovae of massive zero metallicity stars).

3) In the case of 47 Tuc and M71, little change in the range of C and N abundances is observed from the tip of the RGB to the MSTO and below. This is significant as the depth of the convective envelope varies significantly over this interval (see Figure 1) – from a mass fraction of 0.068 at 0.65 $`M_{}`$ (the lower limit of the present sample) and 0.012 at 0.88 $`M_{}`$ (our upper limit) to a maximum fraction of 0.708 on the RGB at MV = +1.8 (VandenBerg, 2003). In a pollution/accretion origin, any surface contamination would therefore be diluted as the convective envelope deepens during RGB ascent, resulting in increasing C and decreasing N abundances on RGB ascent. That this is not observed then requires that a considerable fraction ($`70\%`$ or more) of the mass of a CN-strong star include captured AGB ejecta. As pointed out in Briley et al. (2001), a reservoir of sufficient C-poor/N-rich material may have been produced by more massive AGB stars if the cluster initial mass function is less steep than a Salpeter power law ($`dN=m^\alpha dm,`$ i.e., $`\alpha =1.5`$ rather than 2.35), although without taking into account several other factors such as efficiency with which the low-mass stars incorporate this ejecta, $`\alpha =1.5`$ is at best an upper limit.

Note also that if a substantial fraction of cluster stars were to acquire significant quantities of mass via the accretion of AGB ejecta, one might expect a measurable change in the present day cluster mass function. Moreover, if the accretion process is indeed tied to the dense environment of GCs as suggested by the lack of similar abundance variations among field main sequence stars, a difference in the mass functions in low density environments such as in the Local Group dwarf spheroidal galaxies should be present. Yet the mass functions in the low-density environments of the dwarf spheroidal galaxies Draco (Grillmair et al., 1998) and Ursa Minor (Wyse et al., 2002) appear largely indistinguishable from those of at least the low-metallicity globular clusters, in conflict with the expectations from the AGB ejecta scenario. However, if the accretion of CNO processed material happened during the epoch of star formation, the mass function might remain unchanged.

The present results further highlight two issues in regard to such scenarios. First, the pattern of C and N abundances is the same for M71 and 47 Tuc - see Figure 4 where the C and N abundances for a sample of M71 MSTO/MS stars from Briley & Cohen (2001); Briley, Cohen, & Stetson (2002) are plotted. Yet the abundances found among similar stars in M5 by Cohen et al. (2002) show greater C depletions with similar if not smaller N enhancements among the CN-strong stars (see also Figure 4). 47 Tuc and M71 have nearly identical metallicities (\[Fe/H\] $`0.8`$), and M5 is lower (\[Fe/H\] $`1.2`$). The question is whether we are seeing differences in the initial C and N abundances of the clusters, or a metallicity dependence in the composition of the AGB-star ejecta presumed responsible for the variations among the present-day MS stars. Theory certainly predicts such a dependency (see for example Ventura et al., 2001), although one should also therefore expect larger N over-abundances in M5 relative to 47 Tuc. Clearly, without O, Na, Al, and/or Mg abundances among such stars, and with only three clusters (two of the same metallicity) in the sample, little more can be said.

Also note that despite the similarities in MSTO/MS compositions, the clusters M71 and 47 Tuc have very different structural properties. In the case of an origin for the CN-strong stars in material ejected from AGB stars, one might expect a dependence on some aspect of the cluster environment. Thoul et al. (2002) discuss an accretion model in which centrally concentrated clusters with high escape velocities can trap large quantities of AGB-star ejecta, leading to the accretion of significant mass by low mass members passing near the cluster center before the gas is subsequently stripped during disk crossing (as is apparently taking place in NGC 6779, Hopwood et al., 2000). In the case of 47 Tuc, an additional 0.8$`M_{}`$ could be accumulated by a 1$`M_{}`$ star (0.4$`M_{}`$ for M5) (see their Table 3). Yet both the central concentration and escape velocity for M71 are considerably lower ($`c`$ = 1.15 versus 2.03, and $`v_{esc,core}`$ = 16.7 versus 68.8 km s-1, Harris, 1996; Gnedin et al., 2002, respectively), and this process should be far less efficient in M71 - to the point of being unable to explain the $`70\%`$ or more accreted material required for the Briley & Cohen (2001) results.

Smith (2002) recently revisited a correlation between the fraction of CN-strong stars in a cluster and cluster ellipticity originally noted by Norris (1987). Again, the ellipticity of M71 is 0.000 while that of 47 Tuc is 0.100 (see Smith’s Table 1). However, also note the fraction of CN-strong stars is more than two times greater in 47 Tuc compared to M71 ($`r`$ = 1.8 versus 0.7). While this fraction is based on the distribution seen in larger surveys of RGB stars, neither cluster exhibits evidence of significant composition changes during RGB ascent (see above), thus this ratio also likely reflects the fraction of contaminated stars on the MS. Considering the similarities of such stars in 47 Tuc and M71, it would then appear that we are seeing the action of the same process, resulting in the same composition modifications - just a smaller fraction of the cluster stars were involved in M71. As pointed out by Smith, the CN-strong/ellipticity correlation may be a selection effect – the higher ellipticity clusters in the sample tend also to have higher masses (according to the compilation of Pryor & Meylan, 1993, 47 Tuc has both a considerably higher velocity dispersion and a higher inferred mass than M71). Whether the difference between M71 and 47 Tuc is the result of a smaller central reservoir of AGB ejecta or less mass involved in the tail end of star formation (and how the necessary gas could survive being blown out by supernovae), or some aspect of the cluster dynamics, remains to be seen.

We wish to thank Roger Bell for the use of the SSG program, Don VandenBerg for providing the convective envelope mass fractions, and the anonymous referee for their helpful comments. MMB acknowledges support from the National Science Foundation (under grant AST-0098489), the F. John Barlow professorship, and the UW Oshkosh Faculty Development Program. DH is supported in part through the McKinney postdoctoral fellowship at the University of Wisconsin-Madison. GHS gratefully acknowledges the support of NSF grant AST-0098453.