text:
# CO Isotopes in Planetary Nebulae

## 1 Introduction

Measuring isotopic abundance ratios is one method of probing how stars of various masses process material and how the Galaxy redistributes this material over time—that is, stellar and Galactic chemical evolution. Some of the light elements H, 2H, 3He, 4He, and 7Li, were produced during the era of primordial nucleosynthesis, so these ratios provide additional constraints for Big Bang nucleosynthesis. H $`\mathrm{II}`$ regions or molecular clouds are examples of zero-age objects which are young relative to the age of the Galaxy. Therefore these abundance ratios chronicle the results of billions of years of Galactic chemical evolution. Observations of isotopic ratios in stars and planetary nebulae (PNe) measure locally how stars influence chemical evolution. In particular, PNe probe material which has been ejected from low-mass ($`M2`$$`M_{}`$) to intermediate-mass ($`M25`$$`M_{}`$) stars to be further processed by future generations.

Measurements of $`{}_{}{}^{3}\mathrm{He}`$ in the Galaxy have led to what has been called “The $`{}_{}{}^{3}\mathrm{He}`$ Problem” (Galli et al., 1997). Standard stellar evolution theory predicts that $`{}_{}{}^{3}\mathrm{He}`$ is produced in significant quantities by stars of $`12`$$`M_{}`$ and that PNe with progenitors in that mass range should have $`{}_{}{}^{3}\mathrm{He}/\mathrm{H}`$ $``$ several $`\times 10^4`$ (Iben, 1967; Rood, 1972; Vassiliadis & Wood, 1993; Weiss et al., 1996). (All abundance ratios will be by number unless otherwise noted.) Galactic chemical evolution models which integrate over a range of initial masses predict that $`{}_{}{}^{3}\mathrm{He}`$ is produced in prodigious amounts (Rood et al., 1976; Dearborn et al., 1986). A survey of 6 PNe produced a detection of $`{}_{}{}^{3}\mathrm{He}`$ in NGC 3242 (Rood, Bania, & Wilson, 1992; Balser, Rood, & Bania, 1999) and probable detections in two other sources (Balser et al., 1997) with $`{}_{}{}^{3}\mathrm{He}/\mathrm{H}`$ $``$ $`10^4`$$`10^3`$ as predicted by standard stellar evolution models. This PNe sample was highly biased, however, and was selected to maximize the likelihood of a detection. Observations of $`{}_{}{}^{3}\mathrm{He}`$ in protosolar material (Geiss, 1993), the local interstellar medium (Gloeckler & Geiss, 1996), and Galactic H $`\mathrm{II}`$ regions (Balser et al., 1999) all indicate a $`{}_{}{}^{3}\mathrm{He}/\mathrm{H}`$ abundance ratio of $``$ $`2\times 10^5`$. Although H $`\mathrm{II}`$ regions are expected to be zero-age objects there is no evidence for any stellar $`{}_{}{}^{3}\mathrm{He}`$ enrichment during the last 4.5$`\mathrm{Gyr}`$ (Bania et al., 2002). The observed abundances of $`{}_{}{}^{3}\mathrm{He}`$ throughout the Galaxy when combined with the observations that some stars do indeed produce $`{}_{}{}^{3}\mathrm{He}`$ strongly disagree with standard chemical evolution models (Steigman & Tosi, 1992; Vangioni-Flam et al., 1994; Steigman & Tosi, 1995; Dearborn et al., 1996; Prantzos, 1996; Fields, 1996; Scully et al., 1996; Galli et al., 1997).

Anomalies in the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ abundance ratio in red giant stars has stimulated studies of non-standard stellar evolution models (Charbonnel, 1994, 1995; Wasserburg et al., 1995; Weiss et al., 1996; Sackmann & Boothroyd, 1999; Boothroyd & Sackmann, 1999). All of these models include some form of extra mixing which transports material in the convective envelope to lower regions which are hot enough for further nuclear processing and then transports it back up into the convective envelope. The extra processing that occurs is called cool bottom processing (CBP).

During the first dredge-up on the red giant branch (RGB), the convective envelope reaches regions where $`{}_{}{}^{12}\mathrm{C}`$ had previously been processed into $`{}_{}{}^{13}\mathrm{C}`$ and 14N; thus we expect $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ to decrease from $``$ 90 to 25 for both low-mass and intermediate-mass stars (Iben, 1964; Iben & Renzini, 1984). Observations of carbon in the surfaces of red giant stars reveal lower $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios (5–20) for stars below 2$`M_{}`$ than theoretically expected (e.g., Gilroy 1989). Studies which measure the evolution of the surface $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ as a star ascends the RGB in M67 (Gilroy & Brown, 1991) and in field stars (Charbonnel et al., 1998) indicate that additional processing occurs between the base and tip of the RGB. This has been attributed to extra mixing below the convective envelope.

During the asymptotic giant branch (AGB) the theoretical abundances depend on the interplay between dredge-up events, mass loss, hot bottom burning, and non-standard mixing (CBP). During the third dredge-up on the AGB significant pollution of $`{}_{}{}^{12}\mathrm{C}`$ will increase the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratio if hot bottom burning does not take place (Forestini & Charbonnel, 1997). In low-mass stars, $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ may decrease to $`4`$ if AGB cool bottom processing occurs (Wasserburg et al., 1995). For AGB stars more massive than $`4`$$`M_{}`$, hot bottom burning drives the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ towards its equilibrium value of 3.5 (Frost et al., 1998).

Observations of abundances in PNe measure the cumulative effects of different mixing and nuclear processing events during the entire evolution of their progenitors. In general the abundance of $`{}_{}{}^{3}\mathrm{He}`$ should follow the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratio; that is, if $`{}_{}{}^{3}\mathrm{He}`$ is destroyed then the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratio should decrease. Therefore given the current observational challenges of detecting $`{}_{}{}^{3}\mathrm{He}`$ in PNe, an alternative approach is to measure $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$.

There have been several studies of $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ in the molecular gas in PNe using millimeter transitions of CO (Bachiller et al., 1989; Jaminet et al., 1991; Cox et al., 1992; Bachiller et al., 1997; Palla et al., 2000). In these studies it is assumed that the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios measured in molecular clouds near the PN refer to the abundances in the PN itself. Alternatively, Clegg et al. (1997) and Miskey et al. (2000) have used the $`{}_{}{}^{13}\mathrm{C}`$ $`{}_{1/2}{}^{3}\mathrm{P}_{0}^{o}_{1/2}^1\mathrm{S}_0`$ transition at 1909.597Å to measure the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratio in several PNe. From these studies the measured $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ abundance ratios are between $``$ 10 to 40. Palla et al. (2000; hereafter PBSTG) concluded that in the majority of PNe the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratio is lower than expected from standard stellar models and are consistent with non-standard mixing models which incorporate CBP.

In this paper we present new observations of $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ which enlarge the PN sample. Usually isotope ratios are estimated using optically thin lines, however, since the CO intensities are small, and since the 16O/18O ratios may vary, the $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ lines are used to deduce the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratio. Because of the resulting uncertainty in opacity, it is important that at least two transitions of CO be observed to constrain radiative models used to calculate the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ abundance ratio. The observations are discussed in §2 and the general results are in §3. In §4 the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ abundance ratios are calculated using the assumption that the optical depths of both $`{}_{}{}^{13}\mathrm{CO}`$ and $`{}_{}{}^{12}\mathrm{CO}`$ are small and that the transitions are themalized. More sophisticated methods are explored in §5. In §6 the implications with respect to stellar and Galactic evolution are discussed. A summary is given in §7.

## 2 Observations

The observed sample of PNe are listed in Table 1. The PNe were selected from the CO survey of Huggins & Healy (1989) based on the estimated intensity of $`{}_{}{}^{13}\mathrm{CO}`$ for an expected $`{}_{}{}^{12}\mathrm{CO}/^{13}\mathrm{CO}`$ ratio of 25. Listed are the source name, the Galactic coordinates rounded to the nearest degree, the equatorial coordinates, the observed offset position, the LSR velocity, the angular size, the progenitor main sequence mass, and the morphology. The progenitor masses, $`M_{\mathrm{ms}}`$, are taken from PBSTG for objects in common or determined by first using the central star masses determined by Górny, Stasińska, & Tylenda (1997) or Stasińska, Górny, & Tylenda (1997) and the initial-final mass relation of Herwig (1996). The PNe morphology is defined as round (R), elliptical (E), or bipolar (B) (Balick, 1987; Manchado et al., 1996; Huggins et al., 1996).

An observing log is summarized in Table 2 which includes the spectral transitions, the telescope’s half-power beam-width, and the velocity resolution. Observations of the $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ (J=2$``$1) rotational transitions were made with the National Radio Astronomy Observatory111The National Radio Astronomy Observatory is a facility of the National Science Foundation operated under cooperative agreement by Associated Universities, Inc. 12$`\mathrm{m}`$ telescope on 1997 May and 1997 December. The 200–300$`\mathrm{GHz}`$ dual-polarization SIS receiver was employed. Both the filter-bank and hybrid spectrometers were used simultaneously and were configured to provide spectral resolutions of 1$`\mathrm{MHz}`$ and 2$`\mathrm{MHz}`$ per channel which produce a velocity resolution of 1.4$`\mathrm{km}\mathrm{sec}^1`$ and 2.7$`\mathrm{km}\mathrm{sec}^1`$ per channel, respectively. The spectra were obtained by beam switching plus position switching (BSP). The reference position (OFF) was typically 4 arcmin from the source position (ON). The pointing should be accurate to within $`5`$ arcsec. The data were calibrated on-line to the radiation temperature (T$`{}_{\mathrm{R}}{}^{}{}_{}{}^{}`$) scale which corrects for atmospheric attenuation, radiative loss, and rearward and forward scattering and spillover (Kutner & Ulich, 1981). In this paper we report all intensities in units of the main beam brightness temperature, Tmb, using an efficiency of $`\eta _\mathrm{m}^{}=0.5`$ which converts the radiation temperature to the main beam brightness temperature (Mangum 1999).

Observations of the $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ (J=3$``$2) transitions were made with the Henrich Hertz Telescope (HHT)222The Heinrich Hertz Telescope is operated by the Submillimeter Telescope Observatory on behalf of Steward Observatory and the Max-Planck-Institut für Radioastronomie. during February 1998 and May 1998. The facility 345$`\mathrm{GHz}`$ dual-polarization SIS receiver was employed, although because of technical difficulties only one channel was used. The spectrometers consisted of a 1$`\mathrm{GHz}`$, 1024-channel AOS and a 250$`\mathrm{MHz}`$, 1024-channel AOS with spectral resolutions of 1.0$`\mathrm{MHz}`$ (0.87$`\mathrm{km}\mathrm{sec}^1`$) and 0.40$`\mathrm{MHz}`$ (0.35$`\mathrm{km}\mathrm{sec}^1`$) per channel, respectively. Beam switching was used to produce the spectra using a procedure called WSWITCH with a reference offset from the source by 4 arcmin. (This procedure is effectively identical to BSP at the 12$`\mathrm{m}`$.) The pointing model was determined to be accurate to within 1–2 arcsec for most locations in the sky. The data were calibrated on-line to the T$`{}_{\mathrm{A}}{}^{}{}_{}{}^{}`$ temperature scale which corrects for atmospheric attenuation. The line intensities were scaled to main beam brightness temperatures using several standard calibration sources ($`\chi `$ Cyg, R Aql, and IRC+10216). These intensities are consistent with an assumed efficiency of $`\eta _{\mathrm{mb}}=0.55`$, if the mixer sideband ratios are equal. We estimate that the uncertainties in the intensity scale are between 15–20% (Wilson et al., 2001).

## 3 Results

The results are summarized in Table 3. Listed are the source name, the transition, the peak intensity ($`T_{\mathrm{peak}}`$), the root-mean-squared ($`RMS`$) variations in the spectral baseline, the velocity of the centroid ($`V_{\mathrm{centroid}}`$), the full-width half-maximum of the linewidth based on a single-component Gaussian fit ($`\mathrm{\Delta }V`$), and the integrated intensity ($`I`$). Only $`\mathrm{\Delta }V`$ is based on a Gaussian model. The data were analyzed independently using the spectral line analysis software CLASS, UniPops, and COMB. The majority of the sources were detected in both the $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ J=2$``$1 transitions. Four objects were also detected in the J=3$``$2 transitions which require higher excitation.

Figure 1 summarizes the spectra taken with both telescopes. The first two columns correspond to the NRAO 12$`\mathrm{m}`$ J=2$``$1 transitions while the last two columns correspond to the HHT 10$`\mathrm{m}`$ J=3$``$2 transitions. The spectra have been shifted in velocity to align the profiles relative to the LSR velocity in Table 1. For several sources there is evidence for complicated velocity structure (e.g. multiple peaks, asymmetries, etc.). In these cases this structure is being averaged over the telescope’s beam and thus our analysis refers to the bulk properties of the molecular emitting gas.

For PNe observed at both transitions we measure the J=3$``$2 transition to be significantly more intense than the J=2$``$1 transition. This effect is expected based on solely the radiative transfer in gas with physical conditions characteristic of PNe.

In Figure 2, the radiation temperature is plotted versus the upper J level rotational transition for CO over a range of densities. Across a broad range of physical conditions, the J=3$``$2 transition produces the brightest intensity; only in the lowest volume density-column density models is the J=2$``$1 intensity higher. The generally better coupling between smaller beams (obtained at higher frequencies) and PN structure also suggests that CO searches in PNe may be better carried out in the sub-millimeter regime; atmospheric and hardware limitations may mitigate this effect.

## 4 $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ Abundance Ratios

The derived $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios for our PN sample are listed in Table 4. In columns 3 and 4 the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios are shown for the J=2$``$1 and J=3$``$2 transitions, respectively. Values based on the integrated intensities are expected to be more accurate, but values based on the peak intensities are also given (in parentheses). The $`{}_{}{}^{12}\mathrm{CO}/^{13}\mathrm{CO}`$ abundance ratios are believed to be a direct measure of the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratio; therefore we have assumed that $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ = $`{}_{}{}^{12}\mathrm{CO}/^{13}\mathrm{CO}`$. Chemical fractionation will lower the observed $`{}_{}{}^{12}\mathrm{CO}/^{13}\mathrm{CO}`$ ratios, but this effect should be reduced by the observed high temperatures in the molecular gas ($`2560`$ K, Bachiller et al. 1997). Selective photo-dissociation of $`{}_{}{}^{13}\mathrm{CO}`$ causes an increase in the $`{}_{}{}^{12}\mathrm{CO}/^{13}\mathrm{CO}`$ ratio but this will be counterbalanced by the charge exchange reaction 13C+ \+ 12CO $``$ 12C+ \+ 13CO (Likkel et al., 1988).

The range of $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios are lower than the values predicted by standard stellar models. About half of the objects observed in this study were also included in the sample of PBSTG. In most cases our ratios agree with the PBSTG values shown in the last column in Table 4. The peak main brightness temperatures for our observations are lower in all cases, as expected for the larger beam in our study. However the intensity is not proportional to the beam areas.

These $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios were calculated under the assumptions that the filling fraction for the gas is the same for both isotopomers, both rotational transitions are thermalized to the same temperature, and that the emission lines are optically thin. The first assumption seems reasonable based on the spectra shown in Figure 1. A comparison of the different transitions reveals nearly identical $`V_{\mathrm{LSR}}`$ and line shape characteristics, indicating that similar volumes of gas are sampled by the two transitions. The excitation characteristics of the low-J transitions of CO also supports the assumption of a thermalized, uniform temperature. The small dipole moment of CO results in a low critical density for rotational excitation; thus most PNe have densities high enough for significant CO excitation. The effect of line opacity in CO may be accounted for in multi-transition models of radiation transport.

## 5 Large Velocity Gradient Models

In order to better interpret our results, the radiative transfer in the PNe is examined through the use of the Sobelev or Large Velocity Gradient (LVG) approximation (Sobolev, 1958, 1960). This model provides a numerical solution to the coupled equations of radiative transfer and statistical equilibrium for the rotational transitions of CO (see e.g. Leung & Liszt 1976). In particular, we used the RAD package which is a part of the MIRIAD333Multichannel Image Reconstruction, Image Analysis, and Display software developed by Berkeley-Illinois-Maryland Association. software package (Mundy, 1992). Only the four PNe for which we have measured multiple transitions of $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ emission are considered. For each source a two-dimensional grid composed of volume density ($`n`$) and column density ($`N`$), over a range of physical conditions relevant to molecular gas in PNe, is generated. We adopt a kinetic temperature $`T_\mathrm{k}=60`$ K. This value is towards the high end of temperatures suggested by Bachiller et al. (1997); however, the higher J transitions of CO typically trace warmer regions and the LVG grids are consistent with $`T_\mathrm{k}=60`$ K. Data taken with the HHT (J=3$``$2) have been spatially smoothed to match the resolution of the 12$`\mathrm{m}`$ (J=2$``$1) data. We assume that the emission fills the beam for both pairs of transitions, or has a similar beam filling fraction. We obtained $`\chi ^2`$ values from a comparison of the data to the model grids as a measure of the quality of the fits.

The derived $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios are listed in column 5 of Table 4. The quoted uncertainties were determined by propagating the measured error in the integrated intensities from Table 3. These statistical uncertainties may be much smaller than the systematic uncertainties. For the PNe for which there only exists J=2$``$1 detections the integrated intensity ratios are adopted. Where the line shape is parabolic the line is optically thick so the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios are cited as lower limits. For PNe where we have detected both the J=2$``$1 and J=3$``$2 transitions an LVG analysis is performed and a $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratio determined. Below we discuss each of these sources in more detail.

AFGL 618—Using $`{}_{}{}^{13}\mathrm{CO}`$ we derive $`n`$ = $`6.8\times 10^4`$$`\mathrm{cm}^3`$and $`N`$ = $`6.8\times 10^{15}`$$`\mathrm{cm}^2`$. Using these values, we obtain $`N`$ = $`3.2\times 10^{16}`$$`\mathrm{cm}^2`$ for $`{}_{}{}^{12}\mathrm{CO}`$, yielding a ratio of $`{}_{}{}^{12}\mathrm{CO}/^{13}\mathrm{CO}`$ $`5`$. Overall this value is consistent with the observations. This model overproduces the J=2$``$1 emission by about 20% and the best fit density to the J=2$``$1 is a factor of $``$2 lower than that derived from the $`{}_{}{}^{13}\mathrm{CO}`$. This may suggest that the observed $`{}_{}{}^{12}\mathrm{CO}/^{13}\mathrm{CO}`$ ratio is a lower limit, since the $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ emission regions are not identical, or that density, temperature, and abundance gradients exist which are not isolated in velocity space. Although the derived opacity in the J=3$``$2 and J=2$``$1 lines from the LVG grids is low, the line shape is strongly suggestive of high opacity. We therefore cautiously cite the adopted ratio as a lower limit only.

M1-16—The model for M1-16 produces $`n`$ = $`6.8\times 10^4`$$`\mathrm{cm}^3`$ and $`N`$ = $`6.8\times 10^{15}`$$`\mathrm{cm}^2`$ for $`{}_{}{}^{13}\mathrm{CO}`$, and $`N`$ = $`1.5\times 10^{16}`$$`\mathrm{cm}^2`$ for $`{}_{}{}^{12}\mathrm{CO}`$, yielding $`{}_{}{}^{12}\mathrm{CO}/^{13}\mathrm{CO}`$ = 2.2, very similar to that derived from both the peak and integrated intensities. Although this value is very low, the model suggests that the opacity remains low for both transitions (e.g., $`\tau <0.05`$ for $`{}_{}{}^{12}\mathrm{CO}`$). PBSTG consider this source to be a lower limit, however, presumably because of high opacity. Examining the shape of the spectrum, the profile seems most like the classical flat topped structure of an optically thin, unresolved shell. Therefore we adopt our model-derived value for the ratio in this source.

IRAS 21282+5050—Using $`{}_{}{}^{13}\mathrm{CO}`$ we obtain $`n`$ = $`3.2\times 10^5`$$`\mathrm{cm}^3`$ and $`N`$ = $`1.0\times 10^{15}`$$`\mathrm{cm}^2`$. This produces a column density of $`3.2\times 10^{16}`$$`\mathrm{cm}^2`$ and a ratio of $`{}_{}{}^{12}\mathrm{CO}/^{13}\mathrm{CO}`$ = 32. The volume density for IRAS 21282+5050 is larger than that obtained from the other sources, perhaps indicative of its youth. The shape of its spectral profile is difficult to discern given the signal-noise in our data. For $`{}_{}{}^{12}\mathrm{CO}`$ in the J=3$``$2 line the opacity should be higher but this line seems more parabolic than flat topped. Moreover, the values obtained from the ratios of the integrated intensities and peak intensities are not consistent, most likely due to the low signal-to-noise ratio in the $`{}_{}{}^{13}\mathrm{CO}`$ data. Due to these concerns, we adopt $`{}_{}{}^{12}\mathrm{CO}/^{13}\mathrm{CO}`$ = 32 as a lower limit.

NGC 7293— Numerous models were run for NGC 7293 without obtaining any reasonable fits to the individual $`{}_{}{}^{12}\mathrm{CO}`$ or $`{}_{}{}^{13}\mathrm{CO}`$ data. In addition, the derived physical characteristics for $`{}_{}{}^{13}\mathrm{CO}`$ were highly inconsistent with those from $`{}_{}{}^{12}\mathrm{CO}`$. Since NGC 7293 is an evolved PN, the simplest explanation for these results is that structure in the CO line emitting gas couples differently in the two beams, so the J=2$``$1 and J=3$``$2 regions are distinct and different. The emission profiles from the lines do not suggest a high opacity, however. Since the J=3$``$2 $`{}_{}{}^{13}\mathrm{CO}`$ transition only has a signal-to-noise ratio of $`4`$ we adopt the J=2$``$1 ratio.

## 6 Discussion

The carbon isotopic ratio, $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$, along with other elements are processed during different stages of stellar evolution such as the red giant branch (RGB), the asymptotic giant branch (AGB), and the planetary nebula (PN). Observations have been made during these different stages of evolution which can be used to constrain stellar evolution theory.

Observations of anomalous $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ abundance ratios in RGB field stars (Tomkin, Lambert, & Luck, 1975; Tomkin, Luck, & Lambert, 1976) and open clusters (Gilroy, 1989) have stimulated theoretical studies which introduce an additional mixing process to reduce $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ in stars (Sweigart & Mengel, 1979; Charbonnel, 1994, 1995; Wasserburg et al., 1995). A detailed study of M67 revealed that the anomalously low $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios occurred only for stars at the tip of the RGB (Gilroy & Brown, 1991). Reanalysis of 191 field and cluster giants stars using HIPPARCOS parallaxes has determined that 96% of evolved stars, objects which have passed the luminosity bump, have $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios lower than standard models predict (Charbonnel & Do Nascimento, 1998). Gratton et al. (2000) have measured $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ in 62 field metal-poor stars, including 43 additional stars from the literature. They determine two distinct mixing episodes. The first dredge-up which produces $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios consistent with standard models ($`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ = 20–30); and a second mixing episode occurring after the luminosity bump, when the molecular weight barrier is removed, where $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ is reduced to 6–10. These results are consistent with the theoretical models which predict extra mixing (Charbonnel, 1994, 1995).

Additional processing of carbon and other elements occurs on the AGB. Observations have been made of $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ in AGB stars using mm-wave, optical, and IR spectroscopy. CO mm-wave observations of the circumstellar envelops of carbon and oxygen-rich AGB stars yield $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ values between $`1050`$ (Wannier & Sahai, 1987; Kahane & Jura, 1994; Greaves & Holland, 1997). Evidence for some deep mixing process is suggested from an optical study of carbon stars where Li-rich objects typically have low ($`<15`$) $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios (Abia & Isern, 1997). IR spectroscopy of cool carbon stars produces $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios between $`3070`$ with a few very low ratios, $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ $`4`$ (Lambert et al., 1986). The high $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ values are expected from third dredge-up of $`{}_{}{}^{12}\mathrm{C}`$, while depending on the mass of the star the very low $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios could be explained by hot bottom burning or by deep mixing. In contrast Ohnaka & Tsuji (1996) determine $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios that are 2–3 times smaller than Lambert et al. even though there is a significant overlap in the stars observed. There have been heated discussions about the analysis of spectra using the iso-intensity method (de Laverny & Gustafsson, 1998; Ohnaka & Tsuji, 1998; de Laverny & Gustafsson, 1999; Ohnaka, Tsuji & Aoki, 2000). Nevertheless, it appears that at least some objects do have reduced $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios compared to that expected from standard models and that some extra mixing mechanism is at work.

Observations of $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ in PNe measure how carbon has been processed through both the RGB and AGB phases. Because this processed material escapes into the interstellar medium these $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios can be used to constrain Galactic chemical evolution models. Figure 3 summarizes the current status of $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ abundance ratios in PNe. Plotted are the determined $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios versus the estimated main sequence progenitor mass ($`M_{\mathrm{ms}}`$). Only PNe for which progenitor masses could be found are shown. The results from other CO observations are taken from PBSTG and displayed as triangles. Measurements using the CIII\] transition are shown as stars. Notice that there are five PNe that we have in common with PBSTG. These objects are indicated by open triangles.

For clarity only the (statistical) uncertainties in the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratio are shown in Figure 3. The uncertainties in the derived progenitor masses are significant and are mostly due to estimates of the PN distance and calculations of the progenitor mass from the current central star mass using some initial mass–final mass relationship. These uncertainties can be as large as 50% (e.g, Galli et al. 1997; PBSTG). Therefore care must be used when interpreting these results. The error bars in Figure 3 for $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios are based on the uncertainties in the measured integrated intensity for each transition. Note that the systematic errors may be much larger than these, in some cases.

Our results are consistent with PBSTG with the exception of NGC 6720 where our $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratio is a factor of two smaller. However, the signal-to-noise ratio for the $`{}_{}{}^{13}\mathrm{CO}`$ J=2$``$1 transition is only 5 and we do not have multiple transitions for this PN. Therefore the $`{}_{}{}^{12}\mathrm{CO}`$ transition may be optically thick and our $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ value a lower limit. However from the line shapes these data are consistent with optically thin emission.

The $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios range between 2.2–31 with progenitor masses from $`14`$$`M_{}`$. This indicates that at least some stars do undergo some non-standard processing such as cool bottom processing (CBP). Because of the large uncertainties in determining the progenitor mass it is difficult to make any strong statements about the relationship of $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ with $`M_{\mathrm{ms}}`$. Since CBP is expected to operate only for stars with $`M_{\mathrm{ms}}<2`$$`M_{}`$ the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios should be larger for stars above this mass. Above 4$`M_{}`$ hot bottom burning can drive the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratio towards it equilibrium value of 3.5 (Frost et al., 1998). The low $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ values between 2–4$`M_{}`$ might indicate further processing on the AGB or that we have underestimated the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios. It is interesting to note that the lowest $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ values are classified as bipolar PNe. Bipolar PNe tend to lie closer to the Galactic plane with a wider distribution of central star masses and have higher N/O ratios (Górny, Stasińska, & Tylenda, 1997)

Since most of the PNe in our sample have only been detected in one transition we do not have sufficient constraints to determine the CO line opacity. If the $`{}_{}{}^{12}\mathrm{CO}`$ line is optically thick then the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ values will be underestimated and the points in Figure 3 must be increased. Clearly more $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratios are required to improve the statistics, together with multiple CO transitions to constrain radiative transfer models. Because of the higher kinetic temperatures in PNe, observations in the sub-mm regime may be more useful for this purpose.

A significant effort has been made to detect CO toward the six PNe in which $`{}_{}{}^{3}\mathrm{He}`$ have been measured to directly establish the relationship between non-standard mixing in low-mass stars and the $`{}_{}{}^{3}\mathrm{He}`$ and $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ abundances (PBSTG). The stellar models predict that non-standard mixing (CBP) will lower both $`{}_{}{}^{3}\mathrm{He}/\mathrm{H}`$ and $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$. Based on chemical evolution models more than 90% of low-mass stars must undergo CBP to reconcile the observations of $`{}_{}{}^{3}\mathrm{He}`$ in H $`\mathrm{II}`$ regions (Charbonnel et al., 1998). Only one object with a measured $`{}_{}{}^{3}\mathrm{He}/\mathrm{H}`$ ratio, NGC 6720, has been detected in CO. Since the measured $`{}_{}{}^{3}\mathrm{He}`$ upper limit of $`{}_{}{}^{3}\mathrm{He}/\mathrm{H}`$ = $`0.54\pm 0.29\times 10^3`$ is consistent with both the standard model and with models which include CBP, the current data do not provide any interesting constraints. Note that the progenitor mass of NGC 6720 is slightly above 2$`M_{}`$ which is where CBP should turn off.

CO is typically detected in young PNe where the nebula has not evolved sufficiently to dissociate the molecular material. In contrast, $`{}_{}{}^{3}\mathrm{He}`$ is best detected in evolved, older PNe which tend to fill the telescope’s beam. Therefore it is not a coincidence that there is little overlap between the two samples. As pointed out by PBSTG, observations using the CIII\] transition toward PNe will have a better opportunity to match the $`{}_{}{}^{3}\mathrm{He}`$ samples. Currently, the number of PNe with either CIII\] or $`{}_{}{}^{3}\mathrm{He}`$ detections is sparse. The Green Bank Telescope will offer a significant improvement in our ability to detect $`{}_{}{}^{3}\mathrm{He}_{}^{+}`$ in PNe and thus increase the current sample of one well determined ratio (NGC 3242) and two probable detections. Observations of CO in PNe will complement any CIII\] data. Because of the large uncertainties in estimating the progenitor masses we feel that the number of detections must be significantly increased to build a strong observational case for non-standard mixing in PNe. It is interesting to note that even if CBP occurs in most low-mass stars it will not be a significant factor in the chemical evolution of $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ in the Galaxy since the production of carbon is dominated by intermediate-mass stars (PBSTG). Therefore, measurements of $`{}_{}{}^{3}\mathrm{He}`$ is one of the few methods to probe how low-mass stars influence the chemical evolution of the Galaxy.

## 7  Summary

The following are the main results of the J=2$``$1 and J=3$``$2 millimeter transitions of $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ in 11 Galactic planetary nebulae (PNe).

1. We have detected both $`{}_{}{}^{12}\mathrm{CO}`$ and $`{}_{}{}^{13}\mathrm{CO}`$ emission in 9 PNe for the J=2$``$1 transition using the NRAO 12$`\mathrm{m}`$ telescope and in 4 PNe for the J=3$``$2 transition using the HHT 10$`\mathrm{m}`$ telescope. Assuming $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ = $`{}_{}{}^{12}\mathrm{CO}/^{13}\mathrm{CO}`$ the ratios range from 2.2–31. These results support an extra mixing process in at least some stars.

2. Large velocity gradient (LVG) models have been applied to the 4 PNe where both transitions have been detected. These models were primarily used to determine if the $`{}_{}{}^{12}\mathrm{CO}`$ line is optically thick, an effect which would overestimate the $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratio. A very low $`{}_{}{}^{12}\mathrm{C}/^{13}\mathrm{C}`$ ratio of 2.2 was determined for M1-16 which does not appear to be a lower limit.

3. The J=3$``$2 transition is significantly more intense than the J=2$``$1 transition. This result has been confirmed by models using typical densities measured in PNe. This suggests that CO observations in PN may be better carried out in the sub-millimeter regime.

We thank the technical support staff at both the NRAO 12$`\mathrm{m}`$ telescope and the SMTO. In particular, we thank Harold Butner for help during our observations at the HHT. We acknowledge fruitful discussions with Corinne Charbonnel and Bob Rood about stellar evolution. Lastly, we thank Daniele Galli for communicating his results before publication.


format: text

text:
# The polynomial algorithm for graphs’ isomorphism testing

## 1 Lemma

Let $`A,B`$ be permutation groups on $`V`$, and $`Orb(A),Orb(B)`$ be corresponding systems of orbits on $`V`$ then a partition $`Q=Orb(A)Orb(B)`$ is the system of orbits $`Orb(C)`$ of group $`C`$ generated by subgroups $`A,B`$ (where the union of partitions means the union of its classes containing non-empty intersection).

#### Proof

1. As any permutation $`aA`$ keeps partition $`Orb(A)`$ and any permutation $`bB`$ keeps partition $`Orb(B)`$ then any product of permutations from $`A,B`$ keeps partition $`Q`$.
2. Let $`U`$ be a class of $`Q`$ and points $`x,yU`$ then there exists a sequence of classes $`V_1,W_1,V_2,W_2,V_3,\mathrm{}`$ so that $`V_iOrb(A)`$, $`W_iOrb(B)`$, any two of nearby classes have non-empty intersection and points $`x,y`$ are located in the utmost classes. It follows that there exists a product of elements from $`A,B`$ that transforms $`x`$ in $`y`$.

## 2 Algorithm

Let $`A`$ be a polynomial algorithm that colors vertices of the graph $`G(V)`$ in ordered colors and thus determines an ordered partition on the set $`V`$. Let a graph $`AG`$ obtained from graph $`G`$ by action of algorithm $`A`$ has the same automorphism group as graph $`G`$. Let also this algorithm recognize graphs with trivial automorphism group so that it colors every vertex of such graph in its unique color.

Now we can construct a polynomial algorithm that finds the orbits of automorphism group of any graph $`G(V)`$.

1. Obtain graphs $`R=AG`$, $`R(x_1)=AG(x_1)`$, $`R(x_1,x_2)=AG(x_1,x_2)`$, …, $`R(x_1,x_2,\mathrm{},x_k)=AG(x_1,x_2,\mathrm{},x_k)`$ consecutively, by a consecutive coloring of vertices $`x_1,\mathrm{},x_k`$ from the set $`V`$, until the graph $`R(x_1,\mathrm{},x_k)`$ with regular automorphism group from the initial graph $`G`$ is obtained.

It means that the graph $`R(x_1,\mathrm{},x_k)`$ has a non-trivial automorphism group, but for any other vertex $`x_{k+1}`$ a graph $`G(x_1,\mathrm{},x_{k+1})`$ obtained from the graph $`R(x_1,\mathrm{},x_k)`$ by coloring a vertex $`x_{k+1}`$ has the trivial automorphism group.

Of course, if $`G`$ has trivial automorphism group, the algorithm stops immediately.
2. By testing the graphs $`R(x_1,\mathrm{},x_{k+1})=AG(x_1,\mathrm{},x_{k+1})`$ on isomorphism, construct the system of orbits of the group $`Aut(G(x_1,\mathrm{},x_k))`$.

It is possible to iterate steps 1 and 2 taking every time a new sequence of fix-vertices so as to obtain a new automorphic partition (a system of orbits of some subgroup of $`Aut(G)`$) on $`V`$. By union of automorphic partitions we will every time (accordingly to Lemma) obtain a larger and larger automorphic partition until we will have obtained a system of orbits of automorphism group of initial graph $`G`$.

We can notice that the larger automorphic partition is the more restricted is the possibility for the choice of fix-vertices, because all vertices in the same orbit are identical for such choice.

By finding the isomorphic graphs we find not only automorphic partitions, but also automorphisms. So by constructing the system of orbits of $`Aut(G)`$ we also obtain a defining system of this group.

## 3 How to select fix-vertices

At first we can see that it is sufficient to find at most $`n1`$ automorphic partitions, having non-trivial intersection, in order to obtain the system of orbits of $`Aut(G)`$ as union of these partitions. More precise, if a system of orbits $`Q`$ of some subgroup of $`Aut(G)`$ is already found then it is sufficient to discover further at most $`|Q|1`$ corresponding automorphic partitions.

It can be used different strategies for separation of fix-vertices. One of these can be based on a partial order of set $`V`$ that orders vertices by its number of previous fixations.

The main problem is to verify whether an obtained system of orbits $`Q`$ is $`Orb(Aut(G))`$. Let $`O_1,O_2Q`$. Now we will give the verifying algorithm that verifies whether there exists an automorphism of graph $`G`$ connecting these two orbits (suborbits). It is sufficient to search a corresponding automorphism for any pair of vertices $`o_1O_1`$ and $`o_2O_2`$.

1. By fixation of vertices, obtain the graph $`R(x_1,\mathrm{},x_k)=AG(x_1,\mathrm{},x_k)`$ relative to that the vertices $`o_1,o_2`$ are equivalent (have the same color), but for any other vertex $`x_{k+1}`$ the graph $`R(x_1,\mathrm{},x_{k+1})=AG(x_1,\mathrm{},x_{k+1})`$ distinguishes the vertices $`o_1`$ and $`o_2`$.
2. Obtain the graphs $`T(o_1)R(x_1,\mathrm{},x_k,o_1)`$ and $`T(o_2)R(x_1,\mathrm{},x_k,o_2)`$.

If connecting automorphism for $`o_1`$ and $`o_2`$ exists then the graphs $`T(o_1)`$ and $`T(o_2)`$ are isomorphic and, accordingly to the result of the previous step, have the empty intersection of the same color classes.
3. Using the whole algorithm from beginning, find a system of orbits of automorphism groups of graphs $`T(o_1)`$ and $`T(o_2)`$ and examine these graphs on isomorphism.

It is clear that, in order to finish the verifying algorithm, we need to examine some tree of graphs.

It is clear that the number of levels of consecutive pairs of vertices (the same as $`o_1,o_2`$) of such tree is smaller than $`\mathrm{log}_2n`$, because the graphs $`T(o_1)`$, $`T(o_2)`$ halve the set $`V\backslash \{x_1,\mathrm{},x_k\}`$. And hence the total number of graphs of this tree is smaller than $`n`$.

This consideration is sufficient for polynomiality of whole algorithm.

There is one more property of the base tree of whole algorithm that ascertains the polynomiality of the described algorithm: the total number of graphs of the base tree, that need to be investigated, has a non-principal difference for various cases and is by order not greater than $`n^2`$, because

* the greater the number of levels of the tree, the smaller the number of branches of every level (as for complete graph in the limit, that requires algorithm $`A`$ for about $`n^2`$ graphs )

and vice versa,

* the greater the number of branches that levels of tree have, the smaller the number of levels of the tree (as for graph with regular automorphism group in the limit, that gives $`n`$ graphs for testing of their pairs on isomorphism).

The complexity of algorithm $`A`$ is obviously greater than the complexity of graphs isomorphism testing for graphs with trivial automorphism group. Hence the case “more levels” has greater complexity than the case “more branches”. It follows that the whole algorithm has its greatest complexity for complete graph and this complexity is equal $`n^2|A|`$, where $`|A|`$ is complexity of algorithm $`A`$.

## 4 Conclusion

### 4.1 Algorithm $`A`$

The sequence of graph stabilization algorithms can be constructed in the way B. Weisfeiler described in . We will denote such algorithms as $`A_k`$. The algorithm $`A_k`$ stabilizes the graph $`G`$ by stabilizing the structure on $`V^k`$, generated by this graph. B. Weisfeiler examined the algorithm $`A_2`$ and began to examine the algorithm $`A_3`$.

The given algorithm puts up a question whether there exists a natural $`l`$ so that for all $`kl`$ algorithm $`A_k`$ recognizes the graphs with trivial automorphism group. To author an existence of counter-example for $`l=3`$ is unknown, i.e. an existence of graph with trivial automorphism group in that any two isomorphic triangles have the same number of spanning isomorphic quadrangles.

### 4.2 Perspectives

The algorithms $`A_k`$ present a combinatorial direction of investigation of graphs’ symmetries. This direction for $`k>2`$ cannot have a two-place algebraic interpretation. It is obviously more difficult but substantially stronger for considered problem than a generalization of $`A_2`$ on the algorithms acting on $`V^2`$ that leads to conception of distance-regular graph systematically described in .

On the way to solve the graphs isomorphism problem Author has discovered the original combinatorial objects, that was, most likely, not investigated earlier. These combinatorial objects present the transformation of principal properties of graph symmetries on the set $`C_k`$ of $`k+1`$ classes from $`V^k`$ that leads to conception of being $`C_k`$ assembled or non-assembled in the subset of $`V^{k+1}`$. The simplest example of this representation on $`V^2`$ we obtain in the next matrix notation:

1. Assembled case:

$$C_2^{}=\{\left[\begin{array}{c}12\\ 45\end{array}\right],\left[\begin{array}{c}23\\ 56\end{array}\right],\left[\begin{array}{c}31\\ 64\end{array}\right]\}=\left[\begin{array}{c}123\\ 456\end{array}\right]$$
2. Non-assembled case:

$$C_2^{\prime \prime }=\{\left[\begin{array}{c}12\\ 45\end{array}\right],\left[\begin{array}{c}23\\ 56\end{array}\right],\left[\begin{array}{c}34\\ 61\end{array}\right]\}\left[\begin{array}{c}123\\ 456\end{array}\right]$$

Both sets $`C^{}`$ and $`C^{\prime \prime }`$ have the same projection on the set $`V`$ that consists of subsets

$$\left[\begin{array}{c}1\\ 4\end{array}\right],\left[\begin{array}{c}2\\ 5\end{array}\right],\left[\begin{array}{c}3\\ 6\end{array}\right].$$

The investigation of these objects for greater dimensions shows the way of the problem solution.

### 4.3 Related problems

The existence of the polynomial algorithm of graphs’ isomorphism testing allows to suggest the existence of a full polynomial invariant of a graph and then the existence of a full polynomial invariant of any finite group, because any finite group $`G`$ has a special for this group graph representation as a partition of the set $`G^2`$ generated by the left (right) action of group $`G`$ on the set $`G^2`$. Indeed, this graph can be simplified to a $`d+1`$ color digraph (where one color is empty) with geometrical structure of many-dimensional tore by dimension $`d`$ being equal to a number of group’s generators. The canonization of group’s generating system can allow to identify a group by full invariant of its graph.

On the other hand, the existence of the polynomial algorithm of graphs’ isomorphism testing shows a possible existence of the full polynomial invariant of a partition of $`V^l`$. It gives a different approach for group identification as a permutation representation of minimal degree, acting on $`V`$, and then as automorphism group of its system of orbits on $`V^l`$ for minimal possible $`l`$; where for symmetric group $`l=1`$, for graph automorphism group $`l=2`$, for alternating group $`l=n1`$.

It can be seen that described symmetries are equivalent to symmetry of the set of left (right) cosets of cyclic subgroups of group $`G`$ generated by the group’s generators. This shows that considered approach is different from the group symmetry investigation through the structure of the set of its normal subgroups.

It should be said further if above considerations are correct then it follows that an independent natural full invariant of an abstract group is needed to exist. An investigation of invariants of described symmetrical objects generated by these groups gives the way for obtaining this full invariant.

## Acknowledgements

I would like to express many thanks to Dr. M.Tabachnikov for the proposal (in 1983) to solve this problem, to Dr.V.Grinberg for productive contacts, to Dr. M.Klin for acquainting with contemporary achievements in this field and to my sons Dr. S. Golubchik and R. Golubchik for priceless computer assistance.


text:
# Light Hadron Masses in QCD with Valence Wilson Quarks at 𝛽=6.25 from a Parallel PC Cluster 11footnote 1Supported by the National Science Fund for Distinguished Young Scholars (19825117), National Science Foundation, Guangdong Provincial Natural Science Foundation (990212), Ministry of Education, and Foundation of Zhongshan University Advanced Research Center.

## 1 Introduction

The computational physics group at the Zhongshan University is in a period of rapid development. The group’s interests cover such topics as lattice study of hadron spectroscopy, glueball decay and mixing, QCD at finite density, quantum instontons and quantum chaos using the quantum action. Most of these topics can be studied through Monte Carlo simulation, but can be quite costly in terms of computing power. In order to do large scale numerical investigations of these topics, we built a high performance parallel computer using PC components.

QCD has been accepted to be the most successful gauge theory of strongly interacting particles. Calculation of hadron spectroscopy remains to be an important task of non-perturbative studies of QCD using lattice methods. This paper is one of our first steps in the direction of studying hadron properties. Our reasons for performing simulations with Wilson valence quarks are twofold: First, we are interested in the analysis of hadronic spectroscopy for quenched Wilson fermions; Second, we are interested in exploring the performance of our cluster in actual lattice QCD simulations. We decided to do a spectrum calculation on the $`8^3\times 32`$ lattice at $`\beta =5.7`$ and 5.85, and on the $`12^3\times 36`$ and $`16^3\times 32`$ lattices at $`\beta =6.25`$. The results at smaller $`\beta `$ values are used to compare with those in the literature, while those at a larger $`\beta `$ value, corresponding to smaller lattice spacing $`a`$, are useful for extrapolating the lattice results to the continuum limit. The smearing technique is employed to improve the spectrum data. We hope our data will be an important addition to the lattice study of QCD spectrum.

The remainder of this paper is organized as follows. In Sect. 2, the basic ideas of lattice QCD with Wilson quarks are given. In Sect. 3, we describe the simulation parameters as well as some basic features of our cluster. In Sect. 4, we present new calculation of some light hadron masses. Finally, conclusions and outlook are given in Sect. 5.

## 2 Basic ideas of lattice QCD

Our starting point is the Wilson action

$`S={\displaystyle \frac{\beta }{6}}{\displaystyle \underset{p}{}}\mathrm{Tr}(U_p+U_p^{}2)+{\displaystyle \underset{x,y}{}}\overline{\psi }(x)M_{xy}\psi (y),`$ (1)
where $`\beta =6/g^2`$, and $`U_p`$ is the ordered product of gauge link variables $`U_\mu (x)=e^{ig_{xa}^{xa+\widehat{\mu }a}𝑑x^{}A_\mu (x^{})}`$ around an elementary plaquette. $`a`$ is the lattice spacing, $`\widehat{\mu }`$ is the unit vector and $`M`$ is the fermionic matrix:

$`M_{xy}=\delta _{xy}\kappa {\displaystyle \underset{\mu =1}{\overset{4}{}}}\left[(1\gamma _\mu )U_\mu (x)\delta _{x,y\widehat{\mu }}+(1+\gamma _\mu )U_\mu ^{}(x\widehat{\mu })\delta _{x,y+\widehat{\mu }}\right].`$ (2)
The fermion field $`\psi `$ on the lattice is related to the continuum one $`\psi ^{cont}`$ by $`\psi =\psi ^{cont}\sqrt{a^3/(2\kappa )}`$ with $`\kappa =1/(2ma+8)`$. Then all the physical quantities are calculable through Monte Carlo (MC) simulations with importance sampling. Fermion fields must be integrated out before the simulations, leading to the vacuum expectation value for an operator $`F`$

$`F={\displaystyle \frac{[dU]\overline{F}([U])e^{S_g([U])}\left(detM\right)^{N_f}}{[dU]e^{S_g([U])}\left(detM\right)^{N_f}}}{\displaystyle \frac{1}{N_{\mathrm{config}}}}{\displaystyle \underset{𝒞}{}}\overline{F}[𝒞].`$ (3)
Here $`N_f`$ is the number of flavors, and $`\overline{F}`$ is the corresponding operator after Wick contraction of the fermion fields and the summation is over the gluonic configurations, $`𝒞`$, drawn from the Boltzmann distribution. In quenched approximation, $`det(M)=1`$.

The correlation functions of a hadron is:

$`C_h(t)={\displaystyle \underset{\stackrel{}{x}}{}}O_h^{}(\stackrel{}{x},t)O_h(\stackrel{}{0},0),`$ (4)
where $`O_h(\stackrel{}{x},t)`$ is a hadron operator. For sufficiently large values of $`t`$ and the lattice time period $`T`$, the correlation function is expected to approach the asymptotic form:

$`C_h(t)Z_h[\mathrm{exp}(m_hat)+\mathrm{exp}(m_hatm_haT)].`$ (5)
Fitting the large $`t`$ behavior of the correlation function according to Eq. (5), the hadron mass $`m_h`$ is obtained.

## 3 Simulations

The quenched gauge field configurations and quark propagators were obtained using the recently built ZSU cluster. The full machine has 10 nodes, and each node consists of dual Pentium CPUs. The total internal memory is 1.28GB and the sustained speed is around 2 Gflops. An upgrade is planned. The cluster runs on the Linux operating system. To operate the cluster as a parallel computer, the programmer must design the algorithm so that it appropriately divides the task among the individual processors. We use MPI (Message Passing Interface), one of the most popular message passing standards.

We updated the pure SU(3) gauge fields with Cabibbo-Marinari quasi-heat bath algorithm, each iteration followed by 4 over-relaxation sweeps. The simulation parameters are listed in Tab. 1. The distance between two nearest stored configurations is 100. The auto-correlation time was computed to make sure that these configurations are independent.

The $`u`$ quark and $`d`$ quark are assumed to be degenerate. The quark propagators are calculated in the Coulomb gauge using the independent configurations mentioned above, and conjugate gradient for inversion of the Dirac matrix with preconditioning via ILU decomposition by checkerboards.

At $`\beta =6.25`$, we computed the $`\pi `$, $`\rho `$ and proton and $`\mathrm{\Delta }`$ masses on the $`12^3\times 36`$ and $`16^3\times 32`$ lattices at four values of hopping parameter: $`\kappa =0.1480`$, 0.1486, 0.1492, and 0.1498 with point and smeared sources. In order to compare our results with those by MILC and GF11, the quenched simulations were also performed at $`\beta =5.7`$ and $`\beta =5.85`$ on the $`8^3\times 32`$ lattice. We repeated the quenched simulations of the MILC and GF11 collaborations, using the same set of $`\kappa `$ values  but on a $`8^3\times 32`$ lattice. Detailed data will be presented elsewhere .

## 4 Light hadron masses

To extract masses from the hadron propagators, we must average the correlation function of the hadron over the ensemble of gauge configurations, and use a fitting routine to evaluate the hadron masses $`m_h`$. We determined hadron masses by fitting our data under the assumption that there is a single particle in each channel .

Point source means a delta function, and smeared source means a spread-out distribution (an approximation to the actual wave-function of the quantum state). For example, the simplest operator for a meson is just $`O_h(\stackrel{}{x})=\overline{q}(\stackrel{}{x})q(\stackrel{}{x})`$, i.e. the product of quark and anti-quark fields at a single point. A disadvantage of this point source, is that this operator creates not only the lightest meson, but all possible excited states too. To write down an operator which creates more of the single state, one must “smear” the operator out, e.g. $`_\stackrel{}{y}\overline{q}(\stackrel{}{x})f(\stackrel{}{x}\stackrel{}{y})q(\stackrel{}{y})`$ where $`f(\stackrel{}{x})`$ is some smooth function. Here we choose

$`f(\stackrel{}{x})=N\mathrm{exp}(|\stackrel{}{\mathrm{x}}|^2/\mathrm{r}_0^2),`$ (6)
with $`N`$ a normalization factor. The size of the smeared operator should generally be comparable to the size of the hadron created. There is no automatic procedure  for tuning the smearing parameter $`r_0`$. One simply has to experiment with a couple of choices.

The effective meson mass for $`\pi `$ and $`\rho `$ as a function of time $`t`$ at $`\beta =6.25`$ and on the $`16^3\times 32`$ lattice is depicted in Figs. 1 and 2 respectively. As one sees, the plateau from which one can estimate the mass, is very narrow for point source, due to the reason mentioned above. When the smearing source is used, the width of the plateau changes with the smearing parameter $`r_0`$. We tried many values of $`r_0`$ and found that when $`r_0>16`$, the effective mass is almost independent on $`r_0`$ and the optimal value for $`r_0`$ is 26, where we observe the widest plateau.

The $`\pi `$, $`\rho `$, proton and $`\mathrm{\Delta }`$ masses  at $`\beta =5.7`$ and on $`8^3\times 32`$ are consistent with the previous GF11 data  and those for the mesons at $`\beta =5.85`$ and on $`8^3\times 32`$ are comparable with those in  on the $`16^3\times 32`$ lattice. This means the finite size effects are very small. At $`\beta =5.85`$, we also see that our results for the proton and $`\mathrm{\Delta }`$ masses are bigger than the MILC data and consistent with the finite size behavior analysis in . More detailed results will be given in . In Fig. 3, we show the Edinburgh plot, ($`m_N/m_p`$) vs. ($`m_\pi /m_\rho `$) mass ratios.

In Figs. 4, and 5, we show respectively $`(m_\pi a)^2`$, $`m_\rho a`$, $`m_pa`$ and $`m_\mathrm{\Delta }a`$ as a function of $`1/\kappa `$ at $`\beta =6.25`$ and on the $`16^3\times 32`$ lattice. Assuming that $`(m_\pi a)^2`$ is linear in $`1/\kappa `$, we can compute the critical coupling $`\kappa _c`$ at which the pion becomes massless. We extrapolate the data using:

$$(m_\pi a)^2=A\left(\frac{1}{\kappa }\frac{1}{\kappa _c}\right).$$
(7)
The results for $`\kappa _c`$, and $`m_\rho a`$, $`m_Pa`$ and $`m_\mathrm{\Delta }a`$ at $`\kappa _c`$ on the $`12^3\times 36`$ and $`16^3\times 32`$ are given in Tab.2, where the experimental value for the $`\rho `$ mass has been used as an input. The finite size effects are consistent with . Also, the mass ratios $`m_P/m_\rho `$ and $`m_\mathrm{\Delta }/m_\rho `$ seem to be closer the experimental values, 1.222 and 1.604 on the larger lattice. To compare $`m_P`$ and $`m_\mathrm{\Delta }`$ with experiment, we need to do simulation on a larger lattice volume and carefully study the lattice spacing errors. In this aspect, it might be more efficient to use the improved action and some progress has been reported by some lattice groups in China.

## 5 Conclusions and Outlook

In this paper, we have presented new data on hadron masses in QCD with Wilson valence quarks at $`\beta =6.25`$ on the $`12^3\times 36`$ and $`16^3\times 32`$ lattices. Our results at $`\beta =5.7`$ are consistent and those at $`\beta =5.85`$ are comparable with the results in the literature. We have also made a more careful and systematic study of the smearing method. Such large scale simulations had usually required supercomputing resources , but now they were all done on our recently built cluster for parallel computer. In the task of lattice QCD simulations, we are confident that ZSU’s Pentium cluster can provide a very flexible and extremely economical computing solution, which fits the demands and budget of a developing lattice field theory group. We are going to use the machine to produce more useful results of non-perturbative physics.

Acknowledgements

This work was in part based on the MILC collaboration’s public lattice gauge theory code. (See reference .) We are grateful to C. DeTar, C. McNeile, and D. Toussaint, for helpful discussions.


text:
# 1 Introduction

## 1 Introduction

One of the primary goals of the current generation of colliders is the discovery of the mechanism by which the mass scale of the weak interactions is generated and stabilized. Whether that mechanism is supersymmetry (SUSY), technicolor, extra dimensions or something not yet considered, it will generically require a number of new particle states. One or more of these particles are often stable, at least on experimental time scales, thanks to discrete symmetries in the model or suppressed couplings. The archetypes for such particles are the “LSP” (lightest SUSY particle) in the most SUSY models, and the NLSP (next-to-lightest SUSY particle) in gauge-mediated models.

If a new particle produced at a collider lives long enough to escape the detector, there is no simple way to measure its lifetime. Of utmost importance will be determining whether the lifetime is indeed finite or if the particle is stable on cosmological timescales. It has become standard lore that any stable remnant of new weak-scale physics must be, at most, weakly interacting. The possibility that the stable particle is electrically charged or strongly interacting has long been dismissed based on a combination of astrophysical and cosmological arguments. If such a particle, often called a CHAMP , exists in abundance in our galaxy, then a steady flux of CHAMPs would have been captured by the earth over the course of our planet’s lifetime. Such stable, charged particles would form hydrogen-like atoms which could in turn form superheavy water molecules. These superheavy atoms and molecules would collect in our oceans in significant abundances. Searches for superheavy isotopes of water in ocean water have all yielded negative results . These searches have effectively ruled out charged dark matter at the weak scale.

This result is obviously dependent on the flux of CHAMPs onto the earth. There are two theoretical avenues for determining this flux. The first assumes that these relics represent the majority component of the galactic dark matter halo; then their local abundance can be estimated directly from astrophysical data. From this we easily deduce the rate at which the CHAMPs would be captured and accumulate in our oceans.

The second approach is to use standard cosmological assumptions and techniques to calculate the relic abundance of the CHAMPs from first principles. Then using calculations of dark matter clumping, a local density can be estimated. These calculations also predict far greater abundances of superheavy water than are observed experimentally, even when the abundance is too small to account for the dark matter.

There has been recent interest in whether charged stable relics are ruled out in more cosmologies. Motivated in part by theories with large extra dimensions , studies have been done to calculate relic abundances for the case of extremely low reheat temperatures (as low as an MeV). The surprising result of these studies is that particles with TeV masses can still be created during the reheat process even if the final temperature is around the MeV scale . When applied to electrically charged relics, a bound of roughly 1 TeV on the relic’s mass can be deduced using experimental search bounds .

But can we bound, exclude or search for stable, charged relics without recourse to cosmology? In some models, physics near the TeV scale is so non-canonical that it is unclear whether existing cosmological bounds apply at all (see for example the model of Ref. ). In this paper, we will approach the same problem but from an original point of view. We will consider charged relics created by cosmic rays in the earth’s upper atmosphere. In order to be specific, we will consider SUSY as our template model; however extension of these bounds to any other model of weak-scale physics with a stable charged relic would be trivial.

We will not place a bound on the mass of the SUSY relic itself. Instead we will place a bound on the SUSY scale, or more specifically, the masses of the squarks and gluinos. Direct production of the charged relic in cosmic ray collisions may be very small (for example, the relic could be a stau with only small electromagnetic cross-sections). However, production of any SUSY state would, after a cascade of decays, result in one or more of the charged LSPs. Thus the production cross-section will generally be dominated by production of squarks and gluinos, not the LSP itself.

None of these results depend in any way on new models which attempt to explain cosmic rays at the very highest energies. These results are generated using the usual spectrum of incident protons and heavy nuclei. Our only assumption is that this spectrum has remained essentially unchanged over the last 3 billion years of the earth’s history.

## 2 Cosmic Rays as Particle Colliders

While the energy frontier in particle accelerators is slowing moving higher and higher, collisions with center-of-mass energies in the TeV range have been common throughout the history of the earth in the form of cosmic rays scattering off atmospheric nuclei. While the vast majority of these cosmic rays lose their energy through one or more hard QCD interactions, a small fraction can undergo new-physics interactions. With energies ranging up to $`10^{11}\mathrm{GeV}`$ (in the earth’s frame), the incident cosmic rays are capable of generating center-of-mass energies above $`10^5\mathrm{GeV}`$ when scattering off a proton. And with the exception of the very highest energies, the incident cosmic ray energy spectrum is well-measured.

While these weak interactions in the upper atmosphere are useless for studying short-lived SUSY states, such states will be produced nonetheless. And if R-parity is conserved, those SUSY states will eventually decay down to the LSP. Thus LSPs are produced in our atmosphere in standard particle physics processes, independent of cosmological or astrophysical assumptions

### 2.1 Spectrum and Composition

The all-particle spectrum for cosmic rays is well-measured up to energies of about $`10^{11}\mathrm{GeV}`$. At energies up to about $`10^6\mathrm{GeV}`$ the spectrum follows a simple $`E^{2.7}`$ power law. Around $`10^6\mathrm{GeV}`$, the so-called “knee” is hit at which point the spectrum begins falling more dramatically, as $`E^{3.2}`$. Finally, at energies around $`10^9\mathrm{GeV}`$ (the “ankle”) the spectrum flattens out slightly to $`E^3`$ .

The composition of the incident cosmic rays is not as well known. The primary spectrum is composed of a number of elements including protons, helium, iron, etc. For energies below the “knee” ($`10^6\mathrm{GeV}`$ per nucleus), protons are the most abundant constituents, representing roughly 50% to 80% of the spectrum. At the knee, iron nuclei begin to dominate the spectrum; however, only qualitative details of the composition can be inferred using, for example, the depth of shower maximum. The data seems to indicate that there is a relative rise in heavier elements and then a gradual decrease again as the “ankle” is approached . Above $`10^9\mathrm{GeV}`$ per nucleus, lighter elements appear to be more abundant again, consistent with fragmentation of the heavier elements by the CMBR. Extractions from the data of the average particle number as a function of primary energy vary from 5 to 15 for different experiments . Given this level of ignorance, we assume a fixed fraction of protons (50%) at energies below $`10^6\mathrm{GeV}`$, exponentially decreasing to a constant 1% at energies above $`10^8\mathrm{GeV}`$.

We will also assume that the remainder of the primary composition (the part that is not protons) has on average $`A`$ nucleons/incident nucleus; we will choose $`A=13`$ as a typical experimental value . We will model these nuclei as a collection of $`A`$ loosely bound nucleons, each carrying $`E/A`$ of the incident energy $`E`$ of the entire nucleus; in the center-of-mass frame, the total energy of the collision, $`\sqrt{s}`$, is therefore reduced by a factor $`A^{1/2}`$. To be conservative, we assume that any one nuclei can only participate in a single hard interaction and we do not consider the interactions of the daughter nuclei, if any are formed, at all. We find that nuclei heavier than hydrogen are responsible for only about 10% of the SUSY interactions and therefore our result is almost independent of how we model their interactions. We believe the above choices to be very conservative. Looser constraints with more protons in the spectrum could allow our calculated rates to more than double.

### 2.2 Atmospheric Cascades

Protons 111This discussion is easily generalized to neutrons or heavy nuclei. incident on our atmosphere will usually lose their energy either by QCD or QED processes. Typical QED energy loss rates through brehmstrahlung and ionization are about 2 MeV/(g/cm2 so that charged particles lose only a few $`\mathrm{GeV}`$ traversing 15 km of atmosphere. This is negligible compared to losses in hard QCD interactions.

Energetic primary protons have a nuclear interaction length of $`80`$g/cm2, corresponding to 12 QCD interactions over the depth of the atmosphere . For simplicity, we only consider the first such interaction and neglect secondaries produced in the resulting air showers. That is, we assume that once an incident proton has had a single hard QCD interaction, its energy is degraded beyond the point at which it is kinematically possible to create SUSY particles. In effect, we are discarding the possibility that a secondary could participate in a SUSY interaction, a choice which again results in a conservative bound on the number of SUSY states produced. A complete analysis, unneccesary here, involves solving a set of coupled cascade equations with appropriate boundary conditions, taking into account energy loss processes and energy-dependent cross sections. Our method reduces to analyzing a fixed-target $`pp`$ collision with one of the protons at very high energies.

## 3 The Calculation

Now we present the details of our calculation. After passing through a distance $`dx`$ of the atmosphere with local density $`\rho _{atm}`$, the flux of protons with energies between $`E`$ and $`E+\delta E`$ which have undergone hard QCD interactions and are therefore “lost” for SUSY interactions is simply given by

$$\frac{1}{n_{\mathrm{atm}}(x)}\frac{d\mathrm{\Phi }_p(E)}{dx}=\mathrm{\Phi }_p(E)\sigma _{pp}(E)$$

where $`\mathrm{\Phi }_p`$ is the flux of protons measured in GeV$`^1`$cm$`^2`$s$`^1`$sr-1, $`n_{\mathrm{atm}}`$ is the number density of protons in the atmosphere at a depth $`x`$, and $`E`$ is the energy of the incident proton in the earth’s frame. At the energies under consideration here, the QCD $`pp`$ cross-section, $`\sigma _{pp}`$, is roughly 100 mb, independent of energy. Over that same slice of atmosphere, the total number of SUSY interactions (integrated over all energies), corresponding to the flux of LSPs, is

$$\frac{1}{n_{\mathrm{atm}}(x)}\frac{d\mathrm{\Phi }_{\mathrm{SUSY}}}{dx}=𝑑E\mathrm{\Phi }_p(E)\sigma _{\mathrm{SUSY}}(E).$$

Here $`\sigma _{\mathrm{SUSY}}`$ represents the inclusive cross-section for $`ppX_{\mathrm{SUSY}}`$, independent of whether $`X`$ is squarks, sleptons, etc. ($`\sigma _{\mathrm{SUSY}}`$ is a function of $`s=2m_pE`$.) This system of equations is simple to solve. For a spectrum of only protons, the total number of SUSY interactions in cm$`^2`$s$`^1`$sr-1 is then:

$$\mathrm{\Phi }_{\mathrm{SUSY}}=𝑑E\left(\frac{d\mathrm{\Phi }_p}{dE}\right)\left(\frac{\sigma _{\mathrm{SUSY}}}{\sigma _{pp}}\right)$$

where $`\mathrm{\Phi }_p`$ and $`\sigma _{\mathrm{SUSY}}`$ are implicitly functions of $`E`$. Accounting for the observation that only a fraction $`f_p`$ of the incident cosmic rays are protons, and using our stated assumption that the remaining cosmic rays contain on average $`A`$ nucleons, then the above discussion generalizes to:

$$\mathrm{\Phi }_{\mathrm{SUSY}}=𝑑E\left(\frac{d\mathrm{\Phi }_i}{dE}\right)\frac{1}{\sigma _{pp}}\left[f_p\sigma _{\mathrm{SUSY}}(E)+(1f_p)\sigma _{\mathrm{SUSY}}\left(\frac{E}{A}\right)\right]$$

where $`d\mathrm{\Phi }_i/dE`$ is the all-particle incident flux of cosmic rays.

There is another, simpler way to understand the above calculation. Since we only allow each incident cosmic ray one hard interaction (and only 1 LSP per SUSY interaction), we can simply calculate the probability that that one interaction will be either QCD or SUSY. Since $`\sigma _{\mathrm{SUSY}}\sigma _{pp}`$, the probability of a SUSY interaction is then roughly $`\sigma _{\mathrm{SUSY}}/\sigma _{pp}`$. That probability is then integrated over the entire flux of incident cosmic rays to find the flux of LSPs produced.

The only unknown quantity that remains is $`\sigma _{\mathrm{SUSY}}`$. Since this represents a total cross-section, summed over all possible SUSY final states, there is considerable room for model-dependence. However, we can make several simplifying observations/assumptions. First, because the cosmic rays and the atmospheric nuclei are baryonic, strong interactions should dominate the SUSY production processes. Second, in most “realistic” SUSY models, there is an approximate degeneracy among the strongly-interacting sparticles, that is, among the squarks and the gluinos. These two statements allow us to greatly simplify the calculation by only considering superQCD interactions and by assuming that all squarks and gluinos have a common mass $`M_{\mathrm{SUSY}}`$.

Independent of the details, the parton-level cross-section for SUSY production obviously goes like $`\alpha _s^2/s`$. To be slightly more realistic, we will use actual calculated cross-sections for the range of superQCD processes. The processes that we consider are:

$$\begin{array}{ccccc}\text{(a) }q\overline{q}\stackrel{~}{q}\overline{\stackrel{~}{q}}\hfill & & \text{(c) }qq\stackrel{~}{q}\stackrel{~}{q}\hfill & & \text{(e) }gg\stackrel{~}{g}\stackrel{~}{g}\hfill \\ \text{(b) }gg\stackrel{~}{q}\overline{\stackrel{~}{q}}\hfill & & \text{(d) }q\overline{q}\stackrel{~}{g}\stackrel{~}{g}\hfill & & \text{(f) }qg\stackrel{~}{q}\stackrel{~}{g}\hfill \end{array}$$

The relevant cross-sections have been calculated and tabulated in Ref. ; we confine our calculation to tree level. To go to the proton-proton cross-section we use the CTEQ5M parton distribution functions . The total $`pp`$SUSY cross-section is shown in Figure 1 as a function of $`\sqrt{s}`$ for several choices of $`M_{\mathrm{SUSY}}`$. Over the energy range of interest, the SUSY cross-section is dominated by $`q\overline{q}\stackrel{~}{q}\overline{\stackrel{~}{q}},\stackrel{~}{q}\stackrel{~}{q}`$ and $`qg\stackrel{~}{q}\stackrel{~}{g}`$. There is a subtlety associated with exact versus approximate degeneracy of the squarks and gluinos and so we consider both the case in which $`m_{\stackrel{~}{q}}m_{\stackrel{~}{g}}0`$ and $`m_{\stackrel{~}{q}}m_{\stackrel{~}{g}}M_{\mathrm{SUSY}}`$. We find a difference of about 50% in our calculated number densities as we vary over this range, and so we show a range of limits.

Once a squark or gluino has been produced, it will decay in a cascade down to the the LSP. Presumably, one LSP of positive charge will be produced for each of negative charge. Searches for the latter are more difficult and we will concentrate only on the former. Further, we will assume conservatively that only one positively charged LSP is produced per SUSY interaction, though the number can be significantly higher as the cascades of decays progress. One expects that these positively charged LSPs form superheavy hydrogen by attracting a nearby electron and that this superheavy hydrogen eventually (over the lifetime of the earth) bonds into a superheavy water molecule in the earth’s oceans.

So as a final step we calculate the concentration of superheavy water in the oceans. Considering the age of the oceans, $`t_O`$, to be roughly 3 billions years, and assuming that the flux of cosmic rays has remained essentially unchanged over that time period, the number of superheavy water molecules per usual water molecule in the oceans (the “anomalous concentration”, $`\chi `$) is

$$\chi =7.6\times 10^{13}\left(\frac{t_O}{3\text{Gyr}}\right)\left(\frac{2.6\text{km}}{d_O}\right)\times \frac{\mathrm{\Phi }_{\mathrm{SUSY}}}{\text{cm}^2\text{s}^1\text{sr}^1}$$

where $`d_O`$ is the average depth of the oceans.

The constraints on $`\chi `$ come from searches for superheavy molecules in large samples of water. Experiments then place bounds on $`\chi `$ as a function of the mass of the stable, charged particle. In our approach, this mass in unknown, though it is bounded from above by $`M_{\mathrm{SUSY}}`$. The strongest bound comes the experiment of Smith et al. , who find limits ranging from $`\chi <3.5\times 10^{30}`$ for $`m_{LSP}=100\mathrm{GeV}`$ to $`\chi <8\times 10^{28}`$ for $`m_{LSP}=1.2\mathrm{TeV}`$. Because their bound monotonically weakens as $`m_{LSP}`$ increases, we can place a very conservative bound on $`M_{\mathrm{SUSY}}`$ by setting $`m_{LSP}=M_{\mathrm{SUSY}}`$. For that case, we find $`M_{\mathrm{SUSY}}>230\mathrm{GeV}`$ if the LSP is stable and charged. (Our bound has an uncertainty of roughly $`\pm 10\mathrm{GeV}`$ due to a $`50\%`$ uncertainty that comes from the subtleties for defining $`M_{\mathrm{SUSY}}`$ discussed above.) These results are summarized in Figure 2 where we have shown the predicted anomalous concentration as a function of $`M_{\mathrm{SUSY}}`$ and the experimental bound of Ref. .

## 4 Conclusions

We have shown that there exist bounds from cosmic ray production in the upper atmosphere on charged stable relics (a charged SUSY LSP in particular) which are independent of cosmological constraints. Under the assumption that the incident cosmic ray flux has remained constant over the last 3 billion years, we have calculated a conservative lower bound on the scale of new physics ($`M_{\mathrm{SUSY}}`$), using nothing more than standard particle physics. If SUSY has a stable and charged LSP, then we can place a lower bound on the mass scale of the squarks and gluinos at 230 GeV.

This procedure can easily be extended to other models of weak-scale physics. In such a case, limits similar to those found here could be placed on the masses of new strongly-interacting particles.

Importantly, these bounds will not change if the basic paradigms of cosmology at and below the weak scale are questioned, such as happens in models with large extra dimensions. For example, it is unclear how standard cosmological bounds can be used to constrain models such as that of Ref.  which predicts a light, stable top squark but becomes non-perturbative and higher-dimensional at the TeV scale. The bound presented here should hold even in these highly non-standard cases.

Finally, if a charged particle is discovered which is stable on collider timescales, ruling out or verifying that it is stable on cosmological timescales will require that searches for superheavy water be examined again, though with larger initial samples. Unfortunately, the steeply falling cosmic ray spectrum requires us to go to exponentially larger samples in order to significantly increase our sensitivity. For example, the method used in Ref.  would require an initial heavy water sample equal to that contained in the SNO experiment in order to probe squark masses up to 700 GeV. While such a large-scale search seems unnecessary at present, any future discovery of a new stable, charged particle might require just such an effort. Otherwise there may be no other way to study the stability of that state on long timescales.

## Acknowledgements

CK would like to thank N. Arkani-Hamed, G. Domokos, L. Hall, H. Murayama and J. Poirier for enlightening conversations. This research was supported by the National Science Foundation under grant PHY00-98791.


text:
# Information transfer and nontrivial collective behavior in chaotic coupled map networks

## Acknowledgments

M.G.C and A.P acknowledge support from Consejo de Desarrollo Científico, Humanístico y Tecnológico of Universidad de Los Andes, Mérida, Venezuela. L. Cisneros thanks Consejo Nacional de Investigaciones Científicas y Tecnológicas, Venezuela, for support.


text:
# Magic Structures of Helical Multi-shell Zirconium Nanowires

## Abstract

The structures of free-standing zirconium nanowires with 0.6$``$2.8 nm in diameter are systematically studied by using genetic algorithm simulations with a tight-binding many body potential. Several multi-shell growth sequences with cylindrical structures are obtained. These multi-shell structures are composed of coaxial atomic shells with the three- and four-strands helical, centered pentagonal and hexagonal, and parallel double-chain-core curved surface epitaxy. Under the same growth sequence, the numbers of atomic strands in inner- and outer-shell show even-odd coupling and usually differ by five. The size and structure dependence of angular correlation functions and vibrational properties of zirconium nanowire are also discussed.

In the past decade, there has been tremendous interests on ultrathin metal nanowires from both fundamental low-dimensional physics and technological applications such as nanoelectronic devices . Most of the previous studies are based on the tip-surface contact  or mechanically controllable break junctions , which can be seen as very short metal nanowire. In recent experiments, Takayanagi’s group has successfully fabricated stable gold wires with various diameter of sufficient length suspended between two stands , and the helical multi-shell structures are observed in those ultrathin wires . On the theoretical side, the noncrystalline structures, melting behavior, and electronic properties of ultrathin free standing Pb, Al, Ag nanowires have been investigated by Tosatti and co-workers . By using molecular dynamics-based (MD) genetic algorithm (GA) simulations, our group has studied structural, vibrational, electronic, magnetic properties of gold, titanium, and rhodium nanowires . However, the current knowledge on the detailed structural characters and growth sequences of the metal nanowires is still quite limited. In this paper, we report systematical multi-shell structural growth sequence of zirconium nanowires and their vibrational properties.

In our simulations, the zirconium nanowires with diameters from 0.6 to 2.8 nm are modeled by a supercell with one-dimensional (1D) periodical boundary condition along the wire axis direction. For the most zirconium nanowires studied, the length of 1-D supercell is chosen to be 1.413 nm, which is a reasonable compromise between discovering the helicity in 1D direction and avoiding the nanowires breaking into clusters upon relaxation. To match the periodicity of the multi-shell structures with pentagonal symmetry, several different supercell lengths, i.e., 1.177 nm for S5-1, S5-2 wires (see Fig.1), 1.648 nm for S5-3, S5-4 wires, have been used. The “effective diameter”, which can be controlled by adjusting the number of atoms in the supercell, is used to denote the cross section area of the nanowires . To search the most stable structure of nanowire, we adopt the genetic algorithm (GA) based on MD relaxation . The vibrational densities of states can be then calculated by diagonalizing the dynamical matrix for the optimized configurations.

The interaction between zirconium atoms is described by a tight binding (TB) many body potential , which has successfully reproduced the high-temperature hcp-bcc transition in bulk zirconium . To further check the validity of TB potential in low-dimensional systems with reduced coordination number (CN), we have performed density functional calculations on small zirconium clusters (Zr7, Zr13, Zr15, Zr19) and bulk zirconium solid by using a DMol package  within local density approximation (LDA). In both TB and LDA calculations, the cluster structures are optimized and the similar equilibrium configurations are obtained. Typically, the discrepancy for interatomic distance obtained from tight-binding potential and LDA calculation are within 0.1 Å. The comparison of LDA and TB results on the binding energies of small clusters with equilibrium structures are given in Table I. The satisfactory agreement between LDA and TB results on Zr clusters in Table I demonstrates the validity of the TB potential in the other low-dimensional systems like nanowire.

From the GA simulation, we have obtained numerous isomers for the zirconium nanowires. Among these structural candidates, we found that the ones with closed-shell structures are usually the lowest-energy configurations at a given wire diameter. The energy difference between the lowest-energy closed-shell structures and other open-shell structures is about 0.005eV/atom. For example, the binding energy of helical closed-shell S6-4 wire (see Fig.1) is -5.829 eV/atom, while a low-energy non-helical open-shell configuration at the same size has a binding energy -5.824 eV/atom. In our simulations, we have tested different supercell length and found that the multi-shell arrangements and the helicity of magic closed-shell nanowires do not sensitively depend on the chosen supercell length. On the other hand, those non-magic open-shell wires are either crystalline-like or irregular non-crystalline and have higher binding energies. These results indicate the relative stability of magic closed-shell structure and the possibility of observing these metal nanowires experimentally.

Shown in Fig.1 are five kinds of growth sequences of zirconium nanowires obtained from GA optimization. These multi-shell nanowires are composed of coaxial cylindrical atomic shells. Each shell is usually formed by a number of atomic strands winding up helically side-by-side. The pitches of the helix for the outer and inner shells are different. The lateral surface of the cylindrical shell forms a triangular atomic network. Such helical multi-shell structures have been theoretically predicted for Al, Pb , Au , Ti, and Rh nanowires , and experimentally observed in Au nanowires . In Fig.1, we introduce the index S3-x, S4-x, S5-x, S6-x, and S2-x (x=1-4) to characterize the multi-shell (x: number of shells) growth patterns based on trigonal (S3-x), tetragonal (S4-x), centered pentagonal (S5-x), centered hexagonal (S6-x), and double-chain-core (S2-x) packing, respectively. These twenty structures shown in Fig.1 represent five kind of continuous structural growth sequences of zirconium nanowires. Such systematic multi-shell structural growth patterns in metal nanowires have never been reported before.

With the notation of $`n`$-$`n1`$-$`n2`$-$`n3`$-$`n4`$, we can describe the closed-shell magic nanowires consisting of coaxial shells, each of which is composed of $`n,n1,n2,n3`$, and $`n4`$ helical atomic strands from outer to inner ($`n>n1>n2>n3>n4`$) . As shown in Fig.1, wires S3-1(3), S3-2(8-3), S3-3(13-8-3), and S3-4(18-13-8-3) constitute the sequence of trigonal multi-shell packing. The structure of the thinnest wire S3-1 (diameter D=0.574 nm) is formed by three-strands packing helically (see Fig.1). Similarly, wires S4-1(4), S4-2(9-4), S4-3(14-9-4), and S4-4(19-14-9-4) constitute growth patterns with one-, two-, three-, and four-shell tetragonal packing. This sequence starts from the single-shell wire S4-1 with four-strands helical packing. In both the trigonal and tetragonal multi-shell growth sequence, the outer and inner shells differ by five helical atomic strands.

In contrast to the helical structures found in S3-x, S4-x wires, the wires S5-1(5-1), S5-2(10-5-1), and S5-3(15-10-5-1) under the centered pentagonal multi-shell sequence are not helical, while S5-4(20-15-10-5-1) come back to helical. These structures can be viewed as deformed icosahedral packing. Wires S6-1(6-1), S6-2(11-6-1), S6-3(16-11-6-1), and S6-4(21-16-11-6-1) constitute the sequence of centered hexagonal packing. They are the most ideal cylindrical configurations in all the structures studied. In the double-chain-core sequence of wires S2-1(8-2), S2-2(13-8-2), S2-3(18-13-8-2), and S2-4(23-18-13-8-2), the central core consists a nearly parallel dimer-chain. The surrounding shells are made up of 8, 13, 18, and 23 identical atomic strands respectively. Each strand helically wound with a linear pitch around the cylinder. The differences between the number of atomic strands in outer and inner shells are $`5,5,5,`$ and $`6`$, respectively. The thickest wire in Fig.1 is double-chain-core S2-4 wire with D=2.787 nm. It is worthy noted that, in most growth sequence of zirconium nanowires, the number ($`n,n^{}`$) of atomic stands that forming the inner and outer shells typically have even-odd coupling with $`nn^{}=5`$, with exception at S2-1 and S5-1 wires. These behaviors are similar to those obtained from previous experiment on suspended gold nanowires . In the experiment, the number $`n,n^{}`$ of atomic strands in inner and outer shells also exhibit even-odd coupling with $`nn^{}=7`$, which can be related to the string tension.

The structures, diameters, coordination number, and binding energies of the zirconium nanowires in equilibrium structures are summarized in Table II. The average binding energies $`E`$ and coordination numbers of zirconium nanowires generally increase as the diameter increase. The largest binding energy of the atoms in the thickest wires (S2-4) reach about $`95\%`$ of the bulk value (-6.167 eV/atom)  and that of the thinnest one (S3-1) is about $`75\%`$. If we plot the average binding energies $`E`$ as function of $`1/D`$, all the nanowires from different growth patterns fit well to the same linear dependence as: $`E=E_b\lambda /D`$, where coefficient $`\lambda `$=0.85 eV and the wire diameter D is in unit of nm. The linear dependence with 1/D can be understood by the surface effect . The binding energies of all these multi-shell wires fitting to the same line indicate that they are almost equally stable.

To illustrate the structural characters of zirconium nanowires, the angular correlation functions (ACFs) of these nanowires are calculated and presented in Fig.2. For most structures, two major peaks in the ACFs located around $`60^o`$ and $`120^o`$ are found, which are related to the hexagonal packing. The broad distributions of ACFs in addition to the two major peaks demonstrate that noncrystalline multi-shell structures do not have any definite bond angle. From the ACFs shown in Fig.2, it is clearly seen that discrete feature of the ACF peaks gradually disappear as the wire diameter increases and the bulk character gradually come out. On the other hand, there are considerable differences between the angular correlation functions for the different growth sequences in Fig.2. For instance, the features related to deformed icosahedron can be observed in S5-x (x=1-4) sequence, while perfect icosahedral bond angles are 63.4o, 116.6o, and 180o .

Based on the optimized structures of zirconium nanowires, we now their vibrational properties. The vibrational densities of states of zirconium nanowires in Fig.3 exhibit remarkable dependence on the wire diameter and the corresponding growth classes. Except wires S3-1 and S4-1, the position and shape of the major vibrational peak located around $`2.12.7`$ THz do not sensitively change with the size and growth sequence of wires. This behavior is similar to that found for gold nanowires in our previous work . The first major peak around $`2.12.7`$ THz can be related to the normal vibrational model of hcp bulk zirconium. In an earlier study , the frequency distribution of the bulk zirconium evaluated at room temperature range between 1.2 $``$ 6.1 THz. The first bulk frequency peak were found at 2.1$``$2.3 THz, which is comparable to the major peak ($`2.12.7`$ THz) in zirconium nanowires. Different from the first major peak, although the positions of second peak ($``$ 5 THz) found for the nanowires do not substantially change, its shape and height sensitively vary with the size and growth sequence. Such differences might be understood by the different vibrational couplings between the atomic shells for the nanowires with various sizes and growth classes. In another words, the vibrational modes of zirconium nanowires sensitively depend on the specific growth patterns of the constituent atomic strands and diameters of nanowires. The vibrational bands of the thinner wires (S3-1, S4-1) are rather discrete and molecule-like. Some high frequency modes around $`7.5THz`$ are found in S2-x, S3-x and S4-x sequences. These high frequency vibration modes can be attributed to the contraction of interatomic distance in those wires, which will be further illustrated in the following discussion on zirconium clusters.

In order to further analyze the relations between the vibrational spectra and atomic structures, we have calculated the vibrational densities of states for the zirconium clusters (Zr7, Zr13, Zr15 , Zr19, Zr38, and Zr55) (Fig.4). These vibrational densities of states of clusters are comparable to those of the bulk and nanowires. Because of the structural similarity, the shapes of vibrational densities of states for Zr13 and Zr15 are very closed to those of the pentagonal S5-1 and hexagonal S6-1 wires, respectively. The vibrational models of the clusters can be also understood from their atomic structures. The high frequencies are due to the large force constants, which are the consequence of the contraction of the interatomic distance in the clusters. For instance, the distances between surface atoms and central atom in an icosahedral Zr13 cluster (about 2.96$``$2.98 Å) is smaller than these between surface atoms (3.10 $``$ 3.14 Å), which leads to the high frequency about 6.5 THz. Similar effect are observed in Zr15, Zr19, Zr38, and Zr55. The appearance of some low-frequency modes for zirconium clusters may be attributed to the atoms weakly bonded to the cluster.

From the above discussions, we can make the following conclusions. Multi-shell cylindrical structural growth sequences are found for the zirconium nanowires up to about 3 nm. The number of atomic strands in outer and inner shells exhibit even-odd and the differences between them is usually five for most growth sequence of zirconium nanowires. The angular correlation functions and vibrational properties are related to the wire diameter and growth pattern. The helical multi-shell metal nanowires represent a novel structural form of matter. They are similar but not identical to clusters, which is yet to be fully understood . In principle, the mechanical, thermal, electronic, magnetic, and transport properties of such nanowires may sensitively depend on their atomic structures and helicity. The future studies in these directions are under way. We expect future experiments on transition metal nanowires to validate our arguments.

B. L. Wang and G. H. Wang would like to thank financial support from National Nature Science Foundation of China (No:$`29890210,10023001`$).